#!/usr/common/usg/languages/perl/5.16.0/bin/perl -w
#
# blastMap.pl by Jarrod Chapman <jchapman@lbl.gov> Fri Jun 12 08:49:31 PDT 2009
# Copyright 2009 Jarrod Chapman. All rights reserved.
#


use Getopt::Std;
my %opts = ();
my $validLine = getopts('s:f:m:q:MSFGT:NR', \%opts);
my @required = ("s","f","m",);
my $nRequired = 0;
map {$nRequired += exists($opts{$_})} @required;
$validLine &= ($nRequired == @required);
if ($validLine != 1) {
    print "Usage: ./blastMap.pl <-m merSize> <-f fastaFile> <-s seqTxtFileGlob> <<-q queue>> <<-M(emtime?)>> <<-S(hortReadMode)>> <<-F(astqFormat?)>> <<-G(zipped?)>> <<-T truncateReadLength>>  <<-N(oFormatdb?)>> <<-R(eUseTmpFiles?)>>\n";
    exit;
}

#########
my $meraculousPath = undef;  #  DEFINE THIS FOR YOUR BUILD!
unless (defined($meraculousPath)) {
    die "Edit blastMap.pl to define the path to meraculous source code.\n";
}
#########

my $merSize = $opts{"m"};
my $fastaFile = $opts{"f"};
my $seqFileGlob = $opts{"s"};
my $queue = "normal.c";
if (exists($opts{"q"})) {
    $queue = $opts{"q"};
}
my $gzipped = 0;
if (exists($opts{"G"})) {
    $gzipped = 1;
}
my $memtime = 0;
if (exists($opts{"M"})) {
    $memtime = 1;
}
my $shortReadMode = 0;
if (exists($opts{"S"})) {
    $shortReadMode = 1;
}
my $fastQ = 0;
if (exists($opts{"F"})) {
    $fastQ = 1;
}
my $truncateReads = 0;
if (exists($opts{"T"})) {
    $truncateReads = $opts{"T"};
}

my $formatDatabase = 1;
if (exists($opts{"N"})) {
    $formatDatabase = 0;
}

my $reuseTmpFiles = 0;
if (exists($opts{"R"})) {
    $reuseTmpFiles = 1;
}

my @seqFiles = glob($seqFileGlob);

my $maxReads = 500000;
my $readIndex = 0;
my @filesCreated = ();
my $openFile = undef;
my $fileIndex = 0;
my $outIsOpen = 0;


if ($reuseTmpFiles) {
    @filesCreated = glob("blastMap.tmp.*");
} else {

    foreach my $seqFile (@seqFiles) {
	
	if ($gzipped) {
	    open (S,"gunzip -c $seqFile |") || die "Couldn't open pipe to $seqFile\n";
	} else {
	    open (S,$seqFile) || die "Couldn't open $seqFile\n";
	}
	
	while (my $line = <S>) {
	    if ($readIndex%$maxReads == 0) {
		if (defined($openFile)) {
		    close OUT;
		    $outIsOpen = 0;
		}
		$openFile = "blastMap.tmp.fa.$fileIndex";
		$fileIndex++;
		open (OUT, ">$openFile") || die "Couldn't open $openFile\n";
		print STDERR "Opening temporary fasta file : $openFile\n";
		$outIsOpen = 1;
		push(@filesCreated,$openFile);
	    }
	    if ($fastQ) {
		chomp $line;
		unless ($line =~ /^\@/) {
		    die "Invalid fastQ $line\n";
		}
		
		my $name = undef;
		my $passing = 1;
		# Latest Illumina header format
		if ($line =~ /^\@\S+\s+[12]:[YN]:\d+:/) {
		    my ($id,$end,$fail) = $line =~ /^\@(\S+)\s+([12]):([YN]):/;
		    if ($fail eq "Y") {
			$passing = 0;
		    }
		    $name = "$id/$end";
		} else {
		    # older format
		    ($name) = $line =~ /^\@(\S+)/;
		}
		$line = <S>;
		chomp $line;
		my $seq = $line;
		$line = <S>;
		unless ($line =~ /^\+/) {
		    die "Invalid fastQ $line\n";
		}
		$line = <S>;
		if ($passing) {
		    
		    if ($truncateReads) {
			$seq = substr($seq,0,$truncateReads);
		    }
		    
		    print OUT ">$name\n$seq\n";
		}
	    } else {
		if ($line =~ /.+\/\d+\:.+\:/) { 
		    my ($name,$seq) = $line =~ /(.+\/\d+)\:(.+)\:/; 
		    
		    if ($truncateReads) {
			$seq = substr($seq,0,$truncateReads);
		    }
		    
		    print OUT ">$name\n$seq\n";
		}
	    }
	    $readIndex++;
	}
	close S;
    }
    if ($outIsOpen) {
	close OUT;
    }
}


my $cwd = `pwd`;
chomp $cwd;
my $formatdb = "/usr/common/jgi/aligners/blast/2.2.26/bin/formatdb";
my $formatdbOptions = " -i $fastaFile -p F";

my $blastOut = "blastn.tmp.out";
my $blastn = "/usr/common/jgi/aligners/blast/2.2.26/bin/blastall";
my $blastnOptions = " -b 100 -v 100 -K 100 -p blastn -S 3 -d $cwd/$fastaFile -e 1e-10 -F F -W $merSize -o $blastOut -i "; 

if ($shortReadMode) {
    $blastnOptions = " -b 10 -v 10 -K 10 -p blastn -S 3 -d $cwd/$fastaFile -e 1e-1 -F F -W $merSize -o $blastOut -i "; 
}

my $bv3 = "$meraculousPath/blastView3.pl";
my $bv3Options = " -b $blastOut -p hits=all,HSPs=best,sequence=0 -S $merSize";

my $jobFile = "blastMap.jobs";
my $clusterSubmit = "$meraculousPath/clusterSubmitJobsA.pl";
my $clusterSubmitOptions = " -j $jobFile -q $queue -t 6:0:0 -m 11 ";

if ($memtime) {
    $clusterSubmitOptions .= "-M ";
}

if ($formatDatabase) {
    doIt ("$formatdb $formatdbOptions");
}

open (J, ">$jobFile") || die "Couldn't open $jobFile\n";
foreach my $file (@filesCreated) {
    print J "$blastn $blastnOptions $cwd/$file && $bv3 $bv3Options && rm $blastOut && rm $cwd/$file\n";
}
close J;

doIt ("$clusterSubmit $clusterSubmitOptions");

sub doIt {
    my ($cmd) = @_; 
    print STDERR "$cmd\n";
    system $cmd;
}
