#perl! -w

my $infile=shift(@ARGV); chomp $infile;

my $species=shift(@ARGV); chomp $species;

my $outfile="$infile"."$species"."_extracted_transcripts.fa";
print "$outfile\n";

if(-e "transcripts.fa.nhr"){
  #do nothing
} else{
  system("ncbi-blast-2.2.30+/bin/makeblastdb -in transcripts.fa -dbtype nucl");
}

my $blast_out="blast_out_"."$infile";
system("ncbi-blast-2.2.30+/bin/tblastn -db transcripts.fa -query $infile -out $blast_out -evalue 1e-5 -num_alignments 10 -outfmt 6");
#system("ncbi-blast-2.2.30+/bin/blastx -db human_proteins_ncbi.fa -query temporary_blast_file.fa -evalue 1e-5 -num_alignments 20 -outfmt 6");

my $infile_mod="$blast_out"."_trim";

system("cut -f 2 $blast_out | uniq > $infile_mod");
open IN, $infile_mod or die "cannot open trimmed infile";

my $counter=0;
while (my $line =<IN>){

    chomp $line;
    print "$line\n";
    my $blast_out="$infile"."_"."$species"."_"."$line"."_reciprocal_blast_results";;
    $blast_out =~ s/\//-/g;
    print "$blast_out\n";
    open OUT, ">temporary_blast_results";

    if ($counter==0){
	system("perl getScaffold_samtools.pl transcripts.fa $line > $outfile");
	system("perl getScaffold_samtools.pl transcripts.fa $line > temporary_blast_file.fa");
	my $results=qx(ncbi-blast-2.2.30+/bin/blastx -db human_proteins_ncbi.fa -query temporary_blast_file.fa -evalue 1e-5 -num_alignments 20 -outfmt 6); 

#	my $results = qx(perl web_blast.pl blastx nr temporary_blast_file.fa);
	print OUT "$results\n";
    }
else{
    system("perl getScaffold_samtools.pl transcripts.fa $line >> $outfile");
    system("perl getScaffold_samtools.pl transcripts.fa $line > temporary_blast_file.fa");

 my $results=qx(ncbi-blast-2.2.30+/bin/blastx -db human_proteins_ncbi.fa -query temporary_blast_file.fa -evalue 1e-5 -num_alignments 20 -outfmt 6); 
   # my $results = qx(perl web_blast.pl blastx nr temporary_blast_file.fa);
    
    print OUT "$results\n";

}
    $counter++;

    system("mv temporary_blast_results $blast_out");
}
