##################################
# Finding Polymorphisms          #
#  to Differentiate Species      #
#                                #
# This program will take two     #
# very large fasta flatfiles &   #
# search one against the other   #
# sequence by sequence looking   #
# for local alignment hits that  #
# are top scoring but not        #
# identical. For each alignment  #
# priming sites are identified   #
# to generate an (optionally     #
# nested) primer design.         #
# This is meant for trawling     #
# through genomic/rna data sets  #
# of pairwise taxa in order to   #
# find points of diagnostic      #
# divergence. The default        #
# settings for target primed     #
# sequence length were set with  #
# ancient DNA in mind but can be #
# set to any value you require.  #
##################################

# Version
# 1.0.1
# 09/06/2008
# Robin Allaby, Sarah Palmer, Jay Moore, University of Warwick

# Changes since previous version:
#
# 1.0.1 09/06/2008
# Added fipods.ini functionality for optional supply of default values for parameters
# Improved usefulness of comments
# 


# Known Bugs:
#
# None.
#

# Outstanding requirements:
#
# Make similarity threshold for BLAST into a paramter
# Stop reporting nested primer pairs which don't contain any polymorphisms (this happens occasionally)
# Make program dispose of its own trash
#


# Parameters:
# The script will read default values for all parameters from the defaults file fipods.ini, if it exists. 
# Each line in the defaults file should take the following form:
#
#   parametername=parametervalue
#
# An example would be:
#
#   FASTA_file_1=sequences.fas
#
# Lines beginning with # are ignored
# See below for a list of all the parameters, or look at the sample fipods.ini file that came with this script.
#
# You can alternatively override the defaults with command line parameters in the following strict order.
# Command line parameter values take precedence over default values.
#
#   FASTA_file_1  (mandatory - query file name)
#   FASTA_file_2  (mandatory - search file name)
#   reciprocate   (optional - 1/0 - whether to carry out reciprocal BLAST.  1 is recommended, i.e. yes, and is the default)
#   target_length (optional - goal for length in bp of nested primed sequences, defaults to 150)
#   min_length    (optional - minimum allowed length of nested primed sequences, defaults to 50)
#   max_length    (optional - maximum allowed length of nested primed sequences, defaults to 200)
#   min_target_gc (optional - minimum % GC content of target primed sequences, defaults to 40)
#   max_target_gc (optional - maximum % GC content of target primed sequences, defaults to 65)
#
# If you do specify command line params you must at least include the 2 file names.  If you want to specify any of the optional
# parameters, then you must specify ALL the optional ones up to and including the last one you want. Typical command lines might look like this:
#
#   perl fipods_X.X.X.pl file1.fas file2.fas
#   perl fipods_X.X.X.pl file1.fas file2.fas 1 120 100 200
#   perl fipods_X.X.X.pl file1.fas file2.fas 1 120 100 200 30 70
#
# If fipods.ini doesn't exist, the script expects you to at least supply the filename parameters on the command line (plus any others you want to override).


my $param_file = 'fipods.ini';
my $no_params = 8;
my $qf;
my $sf;
my $reciprocate = 1;
my $target_length = 150;
my $min_length = 50;
my $max_length = 200;
my $min_target_gc = 40;
my $max_target_gc = 65;

my $using_ini = 0;
if (-e $param_file)
{
  $using_ini = 1;  
  open INI,$param_file;
  my @lines = <INI>;
  close INI;
  my @params;
  foreach my $line (@lines)
  {
    #print $line;
    if ($line =~/^\s*#/)
    {
      #print "hello $line\n";
      $line = '';
    }
    $line =~s/\s//g;
    if (length $line > 0)
    {
      #print "$line\n";
      push @params, $line;
    }
  }
  
  if (($#params + 1) != $no_params)
  {
    print "Parameters found in $param_file:\n";
    foreach my $param (@params)
    {
      print "$param\n";
    }
    die "\nWrong number of parameters, $no_params expected.\n";
  }
 
  my %params;
  for (my $i=0; $i<=$no_params - 1; $i++)
  {
    #print "$params[$i]\n";
    my @param = split /=/, $params[$i];
    $params{$param[0]} = $param[1];
    #print "$param[0] $param[1]\n";
  }
  
  if (defined $params{FASTA_file_1})
  {
    $qf=$params{FASTA_file_1};
  }
  else
  {
    undef $qf;
  }
  if (defined $params{FASTA_file_2})
  {
    $sf=$params{FASTA_file_2};
  }
  else
  {
    undef $sf;
  }
  if (defined $params{reciprocate})
  {
    $reciprocate=$params{reciprocate};
  }
  if (defined $params{target_length})
  {
    $target_length=$params{target_length};
  }
  if (defined $params{min_length})
  {
    $min_length=$params{min_length};
  }
  if (defined $params{max_length})
  {
    $max_length=$params{max_length};
  }
  if (defined $params{min_target_gc})
  {
    $min_target_gc=$params{min_target_gc};
  }
  if (defined $params{max_target_gc})
  {
    $max_target_gc=$params{max_target_gc};
  }
}

if ($ARGV[0])
{
  $qf = $ARGV[0]; 
}
elsif (($using_ini == 0) || (!(defined $qf)))
{
  die "Please supply the name of the query FASTA file.\n";
}

if ($ARGV[1])
{
  $sf = $ARGV[1];
}
elsif (($using_ini == 0) || (!(defined $sf)))
{
  die "Please supply the name of the other FASTA file to search.\n";
}

if (defined $ARGV[2])
{
  $reciprocate = $ARGV[2];
}
if (defined $ARGV[3])
{
  $target_length = $ARGV[3];
}
if (defined $ARGV[4])
{
  $min_length = $ARGV[4];
}
if (defined $ARGV[5])
{
  $max_length = $ARGV[5];
}
if (defined $ARGV[6])
{
  $min_target_gc = $ARGV[6];
}
if (defined $ARGV[7])
{
  $max_target_gc = $ARGV[7];
}

#print "$qf\n$sf\n$reciprocate\n$target_length\n$min_length\n$max_length\n$min_target_gc\n$max_target_gc\n";
#exit;


$r = 0; $idhits = 0; $up = 0;
$gameon = 0; $firstid = 0;
$similarity = 85;
@snpscores = ();
$collectseq = 0; $norms = 0; $indels = 0; $biases = 0;
print "making database...\n";
make_blast_databases ($sf, "DNA");
make_blast_databases ($qf, "DNA");
#open (LOG, ">log") or die "Cannot lay my log\n";
open (Q, $qf) or die "Cannot open query file\n";

open (OUT, ">$qf.out") or die "Cannot make output file 1\n";
#open (OUTIN, ">$qf.indel") or die "Cannot make output file 2\n";
#open (OUTBI, ">$qf.bias") or die "Cannot make output file 3\n";
print OUT "F_primer\tR_primer\tSNPs_caught\tdels_caught\tns_caught\tsubs_caught\tindel_caught\tprimed_length\tstart_position\tQuery_name\tQtarget_seq\tSubject_name\tStarget_seq\n";

$counter = 1;
print "searching database .....\n";
my $finished_file = 0;
while (!($finished_file)) {
        if (eof Q) {
          $finished_file = -1;
        }
        else
        {
          $line = <Q>;
        }
	if (($line =~/\>/) || ($finished_file)) { 
		#$queryname = $line;
		if ($gameon == 1) {
		# must have finished loading a sequence in to the file $$,
		# so we must close and search
		
		close O;
		blast_off ($$, $sf, "DNA");
		unlink ("formatdb.log");

		#read blast output
		@output = ();
		open (STUFF, "$$.out") or die "Cannot open output file\n";
		@output = <STUFF>;
		close STUFF;
		$matchno = 1;
@sigseqs = ();
$sigseqsection = 0;
$matchcount = $matchno;
$anothermatchcount = $matchno;
$g = $matchno + 1;
$live = 0;
@taxa = ();
@ids = ();
@sequences = ();
@q = ();

        open (WHO, "$$") or die "Cannot open $$\n";
        @q = <WHO>;
        close WHO;
        $qu_name = $q[0]; $qu_name =~s/\>//g; chomp $qu_name;
        shift @q;
        $qu_seq = join ('', @q);
        $qu_seq =~s/\n//g;
        push (@sequences, $qu_seq);
        #unshift (@taxa, $qu_name);
        @qle = (); @qle = split ('', $qu_seq);
        $query_length = scalar @qle;

foreach $line (@output) {
  #print $line;
                if ($sigseqsection == 1) {
                        if ($matchcount > 0) {
                                        push (@sigseqs, $line);
                                        $matchcount = $matchcount - 1;
                        }
                }

                if ($line =~/Sequences producing significant alignments/) {
                                $sigseqsection = 1;
                }
                if ($line =~/\>/) {
			$firstid = 0;
                        unless ($g > $matchno) {
                        if ($g > 0) {
                                $identities = join (" ", @identities);
                                $ident = process_identities ($identities);
                                push (@ids, $ident);
                        }
                        }
                        if ($anothermatchcount > 0) {
                                $name = capture_name ($line);
				$nameline = $line;
                                push (@taxa, $name);
                                $sigseqsection = 0;
                                @identities = ();
                                $live = 1; $anothermatchcount = $anothermatchcount - 1; $g = $g - 1;
                        } else {
                                $live = 0; 
                        }
                }
                if ($live == 1) {
                if ($line =~/Identities/) {
                  #print "$line";
                 	unless ($firstid == 1) {
			    $firstid = 1;
                            $len = find_length ($line, $query_length);
			    $length = really_find_length ($line); 
                            $id = find_identity ($line);
			    $snpcount = count_snps ($line); $collectseq = 0;
                            if ($id == 100) { $idhits++; push (@identities, "$len\texact");} else { if ($id > $similarity) {push (@identities, "$len\tsimilar");} else { push (@identities, "$len\tdistant");}}
			    if ($id > $similarity) { unless ($id == 100) { unless ($length < 40) { 

				if ($reciprocate == 1) {
				# before we celebrate, we must check the reciprocal hit is indeed to the same query.
				print "checking reciprocal\n";
				open (RS, ">rec") or die "Cannot make reciprocal search\n";
				print RS "$nameline";
				$info = return_blast ($$, 1);
					#print LOG "Info from $$.out:$info\n";
					@info = (); @info = split (/\t/, $info);
					print RS $info[6];
					close RS;
					blast_off ("rec", $qf, "DNA");
					$rechit = return_name_blast ("rec", 1);
                                        chomp $rechit;
                                        $rechit =~s/ //g;

				#	$first = return_identifier_gb ($rechit);
			#		$second = return_identifier_gb ($queryname);
			#	 	print "$rechit against $queryname\n";
			#		print "$first against $second\n";
					if ($queryname =~/$rechit/) { if ($rechit =~/$queryname/) {  
					print "reciprocated\n"; $r++;



					#### categorize interesting target by GC content, indels etc.
					$qseq = (); $sseq = ();
					$qseq = $info[5]; $sseq = $info[6]; @qseq = (); @sseq = ();
					@qseq = split ('', $qseq); @sseq = split ('', $sseq);
				#	print LOG "trial $counter\n";
			#		print LOG "indel: $ind\n";					
					my $gc = gc_content ($qseq);
				#	print LOG "GC: $gc\nqseq: $qseq\n";
					# could also suggest primers/find out if primers can be made easily.	
	
					my $all_primers = primer_design ($qseq, $sseq, $snpcount, $sf, $qf, -1, $target_length, $max_length, $max_target_gc, $min_target_gc);
                                        #print "hello $all_primers hello";
                                        my @all_primers = split (/\n/, $all_primers);
                                        
                                        for my $i(0..$#all_primers)
                                        {
                                                $primers = $all_primers[$i];
        					@primers = (); @primers = split (/\t/, $primers);
                				$primerinfo = scalar @primers; 
                        			if ($primerinfo == 1) {
                                		print "Unsuitable for primer design\n"; $up++;
                                        	} else {
               					#$ind = indel_check ($qseq, $sseq);
                                                #$dels = deletion_check ($qseq, $sseq);
                                                #$ncount = n_check ($qseq, $sseq);

                                                #if ($primers[2] == 0) {
                                                #print "Unsuitable target, primers cannot capture snps\n"; $up++
                                                #} else {
                                                #$type = "norm";  				
					        #if ($ind > 3) {$type = "indel";}
					        #if ($gc > 75) { $type = "bias"; }
					        #if ($gc < 25) {$type = "bias";}
						
						my $targetqseq = substr ($qseq, $primers[8], $primers[7]);
                                                my $targetsseq = substr ($sseq, $primers[8], $primers[7]);
                                                my $subjectname = $nameline;
                                                chomp $subjectname;
                                                print OUT "$primers[0]\t$primers[1]\t$primers[2]\t$primers[3]\t$primers[4]\t$primers[5]\t$primers[6]\t$primers[7]\t$primers[8]\t$queryname\t$targetqseq\t$subjectname\t$targetsseq\n";
					
                                                #$collectseq = 1; 
                                                #if ($type =~/norm/) {
                                                #$norms++;
                                                #print "Normal type\n";
                                                #print OUT "$nameline\n SNP count: $snpcount id: $id length: $len\ntrial $counter\n"; push (@snpscores, $snpcount);
                                                #print OUT "QUERY: $queryname\n";
                                                #print OUT "F primer: $primers[0] R primer: $primers[1] snps in target: $primers[2] target length $primers[3]\n";
                                                #print OUT "Indel: $ind\n";
                                                #print OUT "GC content: $gc\n";
                                                #}
                                                #if ($type =~/indel/) {
                                                #$indels++;
                                                #print "Indel type\n";
                                                #print OUTIN "$nameline\n SNP count: $snpcount id: $id length: $len\ntrial $counter\n"; push (@snpscores, $snpcount);
                                                #print OUTIN "QUERY: $queryname\n";
                                                #print OUTIN "F primer: $primers[0] R primer: $primers[1] snps in target: $primers[2] target length $primers[3]\n";
                                                #print OUTIN "Indel: $ind\n";
                                                #print OUTIN "GC content: $gc\n";
                                                #}
                                                #if ($type =~/bias/) {
                                                #$biases++;
                                                #print "Bias type\n";
                                                #print OUTBI "$nameline\n SNP count: $snpcount id: $id length: $len\ntrial $counter\n"; push (@snpscores, $snpcount);
                                                #print OUTBI "QUERY: $queryname\n";
                                                #print OUTBI "F primer: $primers[0] R primer: $primers[1] snps in target: $primers[2] target length $primers[3]\n";
                                                #print OUTBI "Indel: $ind\n";
                                                #print OUTBI "GC content: $gc\n";
                                                #}
					} }}  }}
					else {
					#$collectseq = 1; print OUT "$nameline\n SNP count: $snpcount id: $id length: $len\ntrial $counter\n";
                                        push (@snpscores, $snpcount);
					#print OUT "QUERY: $queryname\n";

					}}}  }
		}	}
			#if ($collectseq == 1) { unless ($line =~/Identities/) {unless ($line =~/Strand/) { unless ($line =~/Score/) { 
					#if ($type =~/norm/) {
					#print OUT $line;
					#}
					#if ($type =~/indel/) {
					#print OUTIN $line;
					#}
					#if ($type =~/bias/) {
					#print OUTBI $line;
					#}
				#			}}}}
		  	#if ($line =~/Score/) { $collectseq = 0;}
                }

}

print "$counter\n"; $counter++;
#if ($counter == 21) { exit;} 
		open (O, ">$$") or die "Cannot write search file\n";
		print O "$line\n";

		} 
		chomp $line; open (O, ">$$") or die "Cannot write search file\n";
		$queryname = $line; 
		print O "$line\n"; $gameon = 1;
	} else {
		#then it must be sequence, so we must be loading a file
		print O "$line";
	} 



}

close OUT;
#close OUTIN; close OUTBI;
close Q; #close LOG;
print "$r possible targets found with SNPs\n$idhits cases of loci identity between databases\n";
print "$norms normal snp hits\n$indels indel targets\n$biases gc biased targets\n";
print "$up targets were unsuitable for primer design\n";

####################################
# SUBROUTINES                      #
####################################

sub indel_check {
        my $seqone = $_[0]; my $seqtwo = $_[1];
		my $gapone = (); my $gaptwo = (); 
		$gapone = indel_find ($seqone);
		$gaptwo = indel_find ($seqtwo);
		
		if ($gapone > $gaptwo) { return $gapone;} else {return $gaptwo;}
}

sub indel_find {
		my $seqone = $_[0]; 
		my @seqone = (); 
		my $res = (); my $gap = 0; my $gapsize = 0;
		my @gaps = ();  


		@seqone = split ('', $seqone); 
		foreach $res (@seqone) {
				if ($res =~/\-/) {$gap++;} else { 
					#gap ended
						$gapsize = $gap; if ($gapsize > 0) { push (@gaps, $gapsize); $gapsize = 0; $gap = 0;}
				}
		}
		$gapcount = scalar @gaps;
		if ($gapcount > 0) {
			$highgap = $gaps[0];
			foreach $gap (@gaps) { if ($gap > $highgap) { $highgap = $gap;}}
			return $highgap;
		} else { return $gapcount;} 
} 

sub gc_content {
		my $seq = $_[0];
		my @seq = (); my $length = ();
		my $gc = 0; my $res = ();

		@seq = split ('', $seq);
		$length = scalar @seq;
		foreach $res (@seq) { if ($res =~/[gc]/i) {$gc++;}}
		$content = ($gc/$length)*100;
	return $content;
}

sub primer_design {
		my ($seqone, $seqtwo, $snpcount, $dbone, $dbtwo, $nested, $target_length, $max_length, $max_target_gc, $min_target_gc) = @_;
		my @seqone = (); my @seqtwo = (); my $rseqone = (); my $rseqtwo = ();
		my $c = 0; my @fprim = (); my @rprim = (); my $fprim = (); my $rprim = ();		
		my $fprimer = (); my $rprimer = (); my $snipads = ();
		my @answer = (); my $answer = (); my $snps = (); my $primed_length = (); my $length = ();

                @seqone = split ('', $seqone); @seqtwo = split ('', $seqtwo);		
                $length = scalar @seqone;
                $snipads = find_characterise_snps ($seqone, $seqtwo);
		@snipads = split ("\n", $snipads);

                @snp_d = split ("\t", $snipads [0]);
                @snp_n = split ("\t", $snipads [1]);
                @snp_s = split ("\t", $snipads [2]);
                #print "received dels $snipads[0]  Ns $snipads[1]  subs $snipads[2]\n";
                if ($nested)
                {
                        my @nprimer_pairs = ();
                        my $start_position = 0;
                        my $finished = 0;
                        my @primers = ();
                        while (!($finished))
                        {
                                # identifies priming sites (3' to 5')
                                #print "hello2 Start position: $start_position \n";
                                #print "sequences:\n$seqone\n$seqtwo\n";
                                $fprim = make_primer (substr ($seqone, $start_position), substr ($seqtwo, $start_position), $dbone, $dbtwo);
                                @fprim = (); @fprim = split (/\t/, $fprim);
                                $fprim[1] = $fprim[1] + $start_position;
                                $fprim[2] = $fprim[2] + $start_position;
                                $fprimer = $fprim[0]; my $fb = $fprim[1]; my $flb = $fprim[2];
                                
                                if ($fprimer =~/0/)
                                {
                                        $finished = -1;
                                }
                        
                                else
                                {
                                        push @primers, [ @fprim ];
                                        #print "$fprim[0] $fprim[1] $fprim[2]\n";
                                        $start_position = $fb + 1;
                                        if ($start_position > $length - $min_length)
                                        {
                                                $finished = -1;
                                        }
                                }
                        }
                        if ($#primers > 1)
                        {
                          #for my $i(0..$#primers -1)
                          #{
                            #print "hello3 $primers[$i][0] $primers[$i][1] $primers[$i][2]\n";
                          #}
                                $finished = 0;
                                my $nprimerS = 0;
                                my $prev_nprimerE = -1;
                                
                                while (!($finished))
                                {
                                        #this loop is executed once per primer pair found
                                        my $found_End = 0;
                                        my $end = 0;
                                        my $nprimerE = $nprimerS + 1;
                                        #print "hello4 nprimerE: $nprimerE nprimerS: $nprimerS\n";
                                        while (!($found_End))
                                        {
                                                 #this loop is executed once per possible end primer per pair
                                          #print "hello5 $primers[$nprimerE][1]  $primers[$nprimerE][2] $primers[$nprimerS][1]\n";
                                                if ($target_length > ($primers[$nprimerE][2] - $primers[$nprimerS][1]))
                                                {
                                                        #still space for another step?
                                                        $nprimerE++;
                                                        #print "hello6 nprimerE: $nprimerE nprimerS: $nprimerS\n";
                                                        if ($nprimerE > ($#primers - 1))
                                                        {
                                                                #step drops off end of array
                                                                $nprimerE = $nprimerE - 1;
                                                                $found_End = -1;
                                                        }
                                                }
                                                else
                                                {
                                                        if (($nprimerE - $nprimerS) > 1)
                                                        {
                                                                #print "hello7\n";
                                                                #this means we have stepped - good end found within target
                                                                $nprimerE = $nprimerE - 1;
                                                                $found_End = -1;
                                                        }
                                                        else
                                                        {
                                                                #we haven't stepped - need to look at max or move on
                                                                if ($max_length > ($primers[$nprimerE][2] - $primers[$nprimerS][1]))
                                                                {
                                                                        #still space for another step?
                                                                        $nprimerE++;
                                                                        if ($nprimerE > ($#primers - 1))
                                                                        {
                                                                                #step drops off end of array
                                                                                $found_End = -1;
                                                                                $nprimerE = $nprimerE - 1;
                                                                        }
                                                                }
                                                                else
                                                                {
                                                                        #good end found within max
                                                                        #$nprimerE = $nprimerE - 1;
                                                                        $found_End = -1;
                                                                }
                                                                
                                                        }
                                                }
                                        }
                                        #found end or passed end
                                        if (($nprimerE > $nprimerS)  && ($min_length <= ($primers[$nprimerE][2] - $primers[$nprimerS][1])) && ($max_length >= ($primers[$nprimerE][2] - $primers[$nprimerS][1])) && ($nprimerE > $prev_nprimerE))
                                        {
                                                #print "hello8 nprimerS: $nprimerS nprimerE: $nprimerE \n";
                                                #end is good
                                                push @nprimer_pairs, "$nprimerS $nprimerE";
                                                $prev_nprimerE = $nprimerE;
                                        }
                                        else
                                        {
                                                #end is no good (ignore)
                                                $nprimerS = $nprimerS + 1;
                                        }

                                        #get ready for the next pair
                                        $nprimerS = $nprimerS + 1; 
                                        if (($nprimerE == ($#primers - 1)) || ($nprimerS == ($#primers - 1)))
                                        {
                                                #print "hello9\n";
                                                $finished = -1;
                                        }
                                }
                                
                        }
                        else
                        {
                                return 0;
                        }
                        
                        #print "hello10 $#nprimer_pairs\n";
                        for my $i(0..($#nprimer_pairs - 1))
                        {
                                my @answer_line = ();
                                my @this_primer_pair = split ' ',$nprimer_pairs[$i];
                                my $fnprimer = $primers[$this_primer_pair[0]][0];
                                
                                my $rnprimer = revcomp ($primers[$this_primer_pair[1]][0]);
                                my $snps = 0;
                                my $dels = 0;
                                my $ncount = 0;
                                my $subs = 0;
                                my $nstart = $primers[$this_primer_pair[0]][1];
                                my $nend = $primers[$this_primer_pair[1]][2];
                                my $primed_length = $nend - $nstart + 1;
                                my @n_caught = ();
				my @d_caught = ();
				my @s_caught = ();
				my @snpscaught = ();
                		my $sn = ();
                                my $ind = 0;
                                my $nseqone = substr($seqone,$nstart,$primed_length - 1);
                                my $nseqtwo = substr($seqtwo,$nstart,$primed_length - 1);
                                my $gc = gc_content($nseqone);
                                if (($gc <= $max_target_gc) && ($gc >= $min_target_gc))
                                {
                                  $gc = gc_content($nseqtwo);
                                  if (($gc <= $max_target_gc) && ($gc >= $min_target_gc))
                                  {
                                    #print "$nseqone\n$nseqtwo\n";

                                    #foreach $sn (@snipads) { if ($sn > $fb) { if ($sn < $rb) { push (@snpscaught, $sn); }}}
                                    foreach my $snp_n (@snp_n) { if ($snp_n > $nstart) { if ($snp_n < $nend) { push (@n_caught, $snp_n); }}}
                                    foreach my $snp_d (@snp_d) { if ($snp_d > $nstart) { if ($snp_d < $nend) { push (@d_caught, $snp_d); }}}
                                    foreach my $snp_s (@snp_s) { if ($snp_s > $nstart) { if ($snp_s < $nend) { push (@s_caught, $snp_s); }}}
                                    my $longest_indel = indel_check($nseqone, $nseqtwo);

                                    $dels = scalar (@d_caught);
                                    $ncount = scalar (@n_caught);
                                    $subs = scalar (@s_caught);
                                    $snps = $dels + $ncount + $subs;
                                                       

                                    #print "hello11 $fnprimer $rnprimer $fb $rb $snps\n";
 
                                
                                    push (@answer_line, ($fnprimer, $rnprimer, $snps, $dels, $ncount, $subs, $longest_indel, $primed_length, $nstart));
                                    my $answer_line = join ("\t", @answer_line);
                                    #print "hello13 answer line: $answer_line\n";
                                    $answer .= "$answer_line\n";
                                  }
                                }
                        }
                        return $answer;
                }
                else
                {
                        my $indels = 0;
         		my $sn = (); my $rb = (); my $fb = (); my $flb = (); my $rlb = (); 
                        $fprim = make_primer ($seqone, $seqtwo, $dbone, $dbtwo);
                        @fprim = (); @fprim = split (/\t/, $fprim);
                        $fprimer = $fprim[0]; $fb = $fprim[1]; $flb = $fprim[2];

                        if ($fprimer =~/0/) {return $fprimer;}
                        $rseqone = revcomp ($seqone); $rseqtwo = revcomp ($seqtwo);	  			
                        $rprim = make_primer ($rseqone, $rseqtwo, $dbone, $dbtwo);
                        @rprim = split (/\t/, $rprim);
                        $rprimer = $rprim[0]; $rb = $rprim[1]; $rlb = $rprim[2];
                        if ($rprimer =~/0/) {return $rprimer;}
                        if (($flb + $rlb) > $length) {return 0;}
                        $rb = $length - $rb;
                        my @snpscaught = ();
                        foreach $sn (@snipads) { if ($sn > $fb) { if ($sn < $rb) { push (@snpscaught, $sn); }}}

                        $snps = scalar @snpscaught;
                        $primed_length = $rb - $fb;
                        push (@answer, $fprimer, $rprimer, $snps, $indels, $primed_length, $fb);
                        $answer = join ("\t", @answer);
                }
	return $answer;
		
}

sub revcomp()
{
  my ($seq) = @_;
  $seq =~ tr/CAGTcagt/GTCAgtca/;
  return reverse $seq;
}

#sub invert_seq {
#		my $sequence = $_[0];
#		my @sequence = (); my $base = (); my $revbase = (); 
#		my @reverse = (); my $reverse = ();
#@sequence = split ('', $sequence);
#
#@reverse = ();
#
#foreach $base (@sequence) {
#        if ($base =~/[AGCT]/i) {
#                if ($base =~/G/i) {
#                        $revbase = "C";
#                }
#                if ($base =~/C/i) {
#                        $revbase = "G";
#                }
#                if ($base =~/T/i) {
#                       $revbase = "A";
#                }
#                if ($base =~/A/i) {
#                        $revbase = "T";
#                }
#        } else {
#                $revbase = $base;
#        }
#                unshift (@reverse, $revbase);
#
#}
#
#$reverse = join ('', @reverse);
#return $reverse;
#}

sub find_characterise_snps {
		#print "hello ch_snps\n";
		my $seqone = $_[0]; my $seqtwo = $_[1];
		my @seqone = (); my @seqtwo = (); my $c = 0;
		my @n_adds = (); 
                my @del_adds = (); 
                my @sub_adds = (); 
                my @adds = (); my $res = ();
                #print "hello seqone $seqone\n";
                #print "hello seqtwo $seqtwo\n";

                @seqone= split ('', $seqone); @seqtwo = split ('', $seqtwo);

foreach $res (@seqone) 
{
  unless ($res =~/$seqtwo[$c]/) 
  {
    if ((uc($res) =~/N/) || (uc($seqtwo[$c]) =~/N/))
    {
      push (@n_adds, $c);
    }
    elsif (($res =~/-/) || ($seqtwo[$c] =~/-/))
    {
      push (@del_adds, $c);
    }
    else    
    {
      push (@sub_adds, $c);
    }
  }
  $c++;
}
my $n_adds = join ("\t", @n_adds);
my $del_adds = join ("\t", @del_adds);
my $sub_adds = join ("\t", @sub_adds);
#print "sent dels $del_adds  Ns $n_adds  subs $sub_adds\n";
{push (@adds, $del_adds, $n_adds, $sub_adds)};
my $adds = join ("\n", @adds);
	return $adds;
}

sub find_snp_addresses {
		my $seqone = $_[0]; my $seqtwo = $_[1];
		my @seqone = (); my @seqtwo = (); my $c = 0;
		my @adds = (); my $res = ();
		@seqone= split ('', $seqone); @seqtwo = split ('', $seqtwo);
		foreach $res (@seqone) { unless ($res =~/$seqtwo[$c]/) { push (@adds, $c); } $c++; } 
		$adds = join ("\t", @adds);
	return $adds;
}

sub make_primer {
                my $seqone = $_[0]; my $seqtwo = $_[1]; my $dbone = $_[2]; my $dbtwo = $_[3];
                my @seqone = (); my @seqtwo = (); my $c = 0; my $next = 0; $fb = ();
		my $primerscore = 0; my $length = (); my @answer = (); my $answer = ();		
		my @primer = (); my $primer = (); my $lb = 0; my $gc = (); my $f = (); my $check = 0;
                
                @seqone= split ('', $seqone); @seqtwo = split ('', $seqtwo);
		$fb = $c;
		$length = scalar @seqone;
		#print "seqone: $seqone\n seqtwo: $seqtwo\n"; 
		do
                  {
                    if ($seqone[$c] =~/$seqtwo[$c]/)
                    {
                      $lb++;
                      push (@primer, $seqone[$c]);
                      if ($seqone[$c] =~/[cg]/i)
                      {
                        $primerscore = $primerscore + 4;
                      }
                      else
                      {
                        $primerscore = $primerscore + 2;
                      }
                      #print "$c $primerscore\n";
		      if ($primerscore == 64)
                      {
                        $f = shift @primer; $fb++;
                        if ($f =~/[gc]/i)
                        {
                          $primerscore = $primerscore - 4;
                        }
                        else
                        {
                          $primerscore = $primerscore - 2;
                        }
                      }
		      if ($primerscore == 62)
                      {
                      $primer = join ('', @primer); $gc = gc_content($primer); 
		      	#print "primer nascent: $primer  gc: $gc\n";
			if ($gc < 65)
                        {
                          if ($gc > 40)
                          {
			  $check = priming_site_count ($primer, $dbone, $dbtwo);
			  #print "$check priming sites\n";  
			  if ($check < 3)
                            {
                            $next = 1;
                            }
                            else
			    {
                            $f = shift @primer; $fb++; if ($f =~/[gc]/i)
                              {
                              $primerscore = $primerscore - 4;
                              }
                              else
                              {
                              $primerscore = $primerscore - 2;
                              }
                            }
                          }
                          else 
			  {
                          $f = shift @primer; $fb++; if ($f =~/[gc]/i)
                            {
                            $primerscore = $primerscore - 4;
                            }
                            else
                            {
                            $primerscore = $primerscore - 2;
                            }
                          }
                        } 
			else 
			{
                        $f = shift @primer; $fb++; if ($f =~/[gc]/i)
                          {
                          $primerscore = $primerscore - 4;
                          }
                          else
                          {
                          $primerscore = $primerscore - 2;
                          }
                        }
		      } 
		    }
                    else
                    {
                    @primer = (); $primerscore = 0; $fb = $c + 1;
                    if ($c == $length)
                      {
                      $primerscore = 0; $next = 1;
                      }
                    }
		    if ($c == ($length - 1))
                    {
                    $primerscore = 0; $next = 1;
                    }
		    $c++;
                  }
                  until ($next == 1);
		$primer = join ('', @primer);
		if ($primerscore == 0) { $primer = 0;}
		push (@answer, $primer, $fb, $lb); $answer = join ("\t", @answer);
		return $answer;
}

sub priming_site_count {
			my $primer = $_[0]; my $dbone = $_[1]; my $dbtwo = $_[2];
			my $directhits = (); my $more_hits = (); my $total_hits = ();

			open (P, ">primer") or die "Cannot make primer file\n";
			print P ">primer\n$primer";
			close P;

			blast_off ("primer", $dbone, "DNA");
			$directhits = return_direct_hits ("primer", $primer);
			blast_off ("primer", $dbtwo, "DNA");
			$more_hits = return_direct_hits ("primer", $primer);
		$total_hits = $directhits + $more_hits;
	return $total_hits;

}

sub return_direct_hits {
			my $name = $_[0]; my $primer = $_[1];
			my @contents = (); my @primer = (); 
			my $length = (); my $hitl = (); my $line = ();
			my $hitcount = 0;	
			@primer = split ('', $primer);
			$length = scalar @primer;
			open (BOUT, "$name.out") or die "Cannot open blast output\n";
			@contents = <BOUT>;
			close BOUT;
			foreach $line (@contents) {
					if ($line =~/Identities/) {
						if ($line =~/100\%/) {
								$hitl = really_find_length ($line);
								if ($hitl == $length) { $hitcount++;}
						}
					}			

			}
		return $hitcount;
}

sub return_identifier_gb {
			my $line = $_[0]; 
			my @line = (); my $c = 0; my $bit = (); 

			@line = split (/\|/, $line);
			foreach $bit (@line) { if ($bit =~/gb/) { return $line[$c+1]; $c++; } 
			} 
}

sub return_name_blast {
			my $query = $_[0];
			my $level = $_[1];


			my @contents = (); my $line = (); my $name = 0; my $c = 1;
			open (G, "$query.out") or die "Cannot open $query.out\n";
			@contents = ();
			@contents = <G>;
			close G;
			foreach $line (@contents) {
					if ($line =~/\>/) {
							if ($c == $level) {
								$name = $line;
								$c++;
							}
					}
			}
			return $name; 

}
sub count_snps {
		my $line = $_[0];
			my @line = (); my @id_area = (); my $snippet = (); my @snippet = ();
			my $snpcount = ();

                        @line = split (/\,/, $line);
                        @id_area = split (/ /, $line[0]);
			pop @id_area; 
			$snippet = $id_area[$#id_area];
			@snippet = split (/\//, $snippet);
			$snpcount = $snippet[1] - $snippet[0];
		return $snpcount;

}

sub undupe {
		my $taxalist = $_[0]; my @taxa = (); my $found = (); my $check = (); my @check = ();
		my $tasize = (); my $chasize = (); my $dup = 0; my @newlist = (); my $newlist = (); $d = 1;

		@taxa = split (/\t/, $taxalist);
foreach $taxon (@taxa) {
                $found = 0;
		$dup = 0;
                foreach $check (@taxa) {
                        if ($taxon =~/$check/) {
                            $found = $found + 1;
                          #  print TEST "$taxon $check $found\n";
                            @taxon = (); @check = ();
                            @taxon = split ("", $taxon); @check = split ("", $check);
                            $tasize = scalar (@taxon); $chasize = scalar (@check);
                            if ($tasize == $chasize) {
                                   if ($found > 1) {
					$dup = 1;
                                   }
                            }
                         }
                }
	if ($dup == 1) {
                                        unshift (@taxon, $d); $d = $d + 1;
                                        $taxon = join ('', @taxon); push (@newlist, $taxon);
	} else {
		push (@newlist, $taxon);
	}

}

	$newlist = join ("\t", @newlist);

	return $newlist;
}

sub tough_remove {
		my $dead_one = $_[0];
		my @argv = ();
		@argv = ("rm $dead_one");
		system(@argv) == 0 or die "Cannot remove this file\n";

}

sub move_file_up {
			my $file = $_[0];
			my $folder = $_[1];

			my @argv = ();

			@argv = ("mv $folder/$file ./");
			system(@argv) == 0 or die "Cannot move file up\n"; 

}

sub test_sequence {
			my $testee = $_[0];
			my $tester = $_[1];
			my $length_limit = $_[2];
			my $l = (); my $m = (); my @t = (); my @p = (); $a = (); $w = 0;

			@t = split ('', $testee); @p = split ('', $tester);
			$l = scalar @t; $m = scalar @p; $a = ($l/$m) * 100;
			if ($a > $length_limit) {
				$w = 1;
			}

			return $w;
			
			
			
}

sub prepare_prot {
		my $alignment = $_[0]; 
		my @contents = ();
		open (ALIGN, $alignment) or die "Cannot open alignment\n";
		@contents = <ALIGN>; close ALIGN;
		open (ALIGN, ">$alignment") or die "Cannot overwrite alignment\n";
		foreach $line (@contents) {
			$line =~s/\(/_/g; $line =~s/\)/_/g; $line =~s/\:/_/g; $line =~s/\;/_/g; 
			$line =~s/\,/_/g; $line =~s/\[/_/g; $line =~s/\]/_/g;
			print ALIGN $line;
		}
		close ALIGN;


}

sub ml_prot_tree {
		my $alignment = $_[0]; 
		my @argv = ();
		@argv = ("(echo $alignment; echo y;) | (./program/proml)");
		system(@argv) == 0 or die "Cannot make tree\n"; 

}

sub ml_tree {
		my $alignment = $_[0];
		my @argv = ();
		
		@argv = ("./program/fastDNAml<$alignment > output");
		system(@argv) == 0 or die "Cannot make tree\n";
}


sub single_space {
                   my $line = $_[0];
                   my $next = 0;
                      do {
                         if ($line =~/  /) {
                                   $line =~s/  / /g;
                         } else {
                              $next = 1;
                         }
    
                         } until ($next == 1);
                    return $line;
} 



sub prepare_for_DNAml {
			my $alignment = $_[0]; my @contents = (); my @line = (); $first_line = (); $first = ();

			open (ALIGN, "$alignment") or die "Cannot open alignment file\n";
			@contents = <ALIGN>;
			close ALIGN;

			$first = $contents[0]; chomp $first;
			shift @contents; 

			push (@line, $first, " K\nK 5\n"); $first_line = join ('', @line);
			unshift (@contents, $first_line);

			open (ALIGN, ">$alignment") or die "Cannot overwrite alignment file\n";
			print ALIGN @contents;
			close ALIGN;
}

sub align {
		$input = $_[0];
		my @argv = ();


		@argv = ("mv $input program");
		system(@argv) == 0 or die "Could not move alignment\n";
	#	@argv = ();
#	#	@argv = ("cd program"); system(@argv) == 0 or die "Could not get into program folder\n";
		@argv = ();
		@argv = ("(echo 1; echo program/$input; echo 2; echo 9; echo 4; echo 1; echo ; echo 1; echo ; echo ; echo x; echo ; echo x;) | (./program/clustalw)");
		system(@argv) == 0 or die "Could not align sequences\n";

}

sub return_blast {
                  my $taxon = $_[0];
                  my $level = $_[1];
                  my @blout = ();
                  my $hitcount = 0;
                  my $line = ();
                  my $collect = 0;
                  my @collection = ();
                  my $bagsize = ();
                  my $residue_register = 0;
                  my $residue_target = 0;
                  my @posline = ();
                  my @specpos = ();
                  my @subject_line = ();
                  my $residues = ();
                  my @residues = ();
                  my $residue_count = ();
                  my @subject_residues = ();
                  my @query_line = ();
                  my @query_residues = ();
                  my $distance = ();
                  my $data = ();
                  my @answer = ();
                  my $answer = ();
                  my $subject_residues = ();
                  my $query_residues = ();
                  my $posfound = 0;
                  my @subject_regions = ();
                  my $start = 0;
                  my $end = 0;
                  my @hitnumbers = ();
                  my $hitnumber = ();

                  open (BLASTOUT, "$taxon.out") or die "Cannot open blast output for $taxon\n";
                  @blout = <BLASTOUT>;
                  close BLASTOUT;

                  foreach $line (@blout) {
                          if ($line =~/\>/) {
                              $hitcount = $hitcount + 1;
                          }
                          if ($hitcount == $level) {
                               unless ($posfound == 1) {
                                 $collect = 1;
                               }
                          } else {
                                 $collect = 0;
                          }

                          if ($line =~/Number of sequences better than 10.0/) {
                                       @hitnumbers = split (/\:/, $line);
                                       $hitnumber = $hitnumbers[1];
                                       $hitnumber =~s/ //g;
                                       chomp $hitnumber;
                          }

                          if ($collect == 1) {
                              push (@collection, $line);
                              if ($line =~/\>/) {
                                          $subject_taxon = $line;
                                          $subject_taxon =~s/\>//g;
                                          chomp $subject_taxon;
                              }

                              if ($line =~/Identities/) {
                                     if ($posfound == 0) {
                                          @posline = split (/\//, $line);
                                          @specpos = split (/ /, $posline[1]);
                                          $residue_target = $specpos[0];
                                          $posfound = 1;
                                     } else {
                                          $collect = 0;
                                     }
                              }
                              if ($line =~/Sbjct/) {

                                         if ($residue_register < $residue_target) {
                                          @subject_line = ();
                                          $line = single_space ($line);
                                          @subject_line = split (/ /, $line);
                                          push (@subject_regions, $subject_line[1], $subject_line[3]);
                                          @subject_line = ();
                                          @subject_line = split (/\:/, $line);
                                          $residues = $subject_line[1];
                                          $residues =~s/[1234567890]//g;
                                          $residues =~s/ //g;
                                          chomp $residues;
                                          push (@subject_residues, $residues);
                                          #$residues =~s/\-//g;

                                          @residues = ();
                                          @residues = split ('', $residues);
                                          $residue_count = scalar @residues;
                                          $residue_register = $residue_register + $residue_count;
                                          }
                               }
                               if ($line =~/Query/) {
                                          if ($residue_register < $residue_target) {
                                           @query_line = ();

                                           @query_line = split (/\:/, $line);
                                           $residues = $query_line[1];
                                           $residues =~s/[1234567890]//g;
                                           $residues =~s/ //g;
                                           chomp $residues;
                                           push (@query_residues, $residues);
                                          }
                               }
                              if ($residue_register == $residue_target) {
                                  $collect = 0;
                              }
                          }
                   }


                   $bagsize = scalar @collection;
                   if ($bagsize == 0) {
                            $result = 0;
                            return $result;
                   } else {

                          if ($blout[0] =~/BLASTP/) {
                              $data = "Protein";
                          }
                          if ($blout[0] =~/BLASTN/) {
                              $data = "DNA";
                          }
                          $query_residues = join ('', @query_residues);
                          $subject_residues = join ('', @subject_residues);
                          $start = shift @subject_regions;
                          $end = pop @subject_regions;
                          chomp $end;
                          $subject_taxon =~s/ //g;
                        #  $distance = get_pairwise_distance ($taxon, $subject_taxon, $query_residues, $subject_residues, $data);
			  $distance = 0;         
                          push (@answer, 1, $subject_taxon, $distance, $start, $end, $query_residues, $subject_residues, $hitnumber);
                          $answer = join ("\t", @answer);
                          return $answer;
                   }

}



sub check_length {
			my $problem = $_[0]; my $desire = $_[1];
			my @problem = (); my $extent = ();
			my $solution = 0;

			@problem = split ('', $problem);

			$extent = scalar @problem;
			if ($extent == $desire) {
				$solution = 1;
			}
	return $solution;

}

sub process_identities {
			my $identities = $_[0];
			my @identities = (); my $exact = 0; my $similar = 0; my $distant = 0;
			my $exact_short = 0; my $extremely_short_exact = 0; my $similar_short = 0; my $extremely_short_similar = 0;
			my $distant_short = 0; my $extremely_short_distant = 0;
			my @answer = (); my $answer = (); my $clock = 0;
			@identities = split (/ /, $identities);
			foreach $ug (@identities) {
					if ($ug =~/short/) {
						if ($ug =~/extremely/) {
							if ($ug =~/exact/) { $extremely_short_exact = $extremely_short_exact + 1; }
                                        	        if ($ug =~/similar/) { $extremely_short_similar = $extremely_short_similar + 1; }
                                        	        if ($ug =~/distant/) {$extremely_short_distant = $extremely_short_distant + 1; }
						} else {
							if ($ug =~/exact/) { $exact_short = $exact_short + 1; }
        	                                        if ($ug =~/similar/) { $similar_short = $similar_short + 1; }
	                                                if ($ug =~/distant/) {$distant_short = $distant_short + 1; }
						}
					} else {
						if ($ug =~/exact/) { $exact = $exact + 1; } 
						if ($ug =~/similar/) { $similar = $similar + 1; }
						if ($ug =~/distant/) {$distant = $distant + 1; }
					}
			}	
			
			if ($exact > 0) {
				$clock = 1;
				if ($exact > 1) {
					push (@answer, "\t\t$exact exact matches\n");
				} else {
					push (@answer, "\t\texact match\n");
				}
			}
                        if ($similar > 0) {
				if ($clock == 1) {push (@answer, "\t");}
                                if ($similar > 1) {
                                        push (@answer, "\t\t$similar similar matches\n");
                                } else {
                                        push (@answer, "\t\tsimilar match\n");
                                }
                        }
                        if ($distant > 0) {
				if ($clock == 1) {push (@answer, "\t");}
                                if ($distant > 1) {
                                        push (@answer, "\t\t$distant distant matches\n");
                                } else {
                                        push (@answer, "\t\tdistant match\n");
                                }
                        }
                        if ($exact_short > 0) {
                                $clock = 1;
                                if ($exact_short > 1) {
                                        push (@answer, "\t\t$exact_short short exact matches\n");
                                } else {
                                        push (@answer, "\t\tshort exact match\n");
                                }
                        }
                        if ($similar_short > 0) {
                                if ($clock == 1) {push (@answer, "\t");}
                                if ($similar_short > 1) {
                                        push (@answer, "\t\t$similar_short short similar matches\n");
                                } else {
                                        push (@answer, "\t\tshort similar match\n");
                                }
                        }
                        if ($distant_short > 0) {
                                if ($clock == 1) {push (@answer, "\t");}
                                if ($distant_short > 1) {
                                        push (@answer, "\t\t$distant_short short distant matches\n");
                                } else {
                                        push (@answer, "\t\tshort distant match\n");
                                }
                        }
                        if ($extremely_short_exact > 0) {
                                $clock = 1;
                                if ($extremely_short_exact > 1) {
                                        push (@answer, "\t\t$extremely_short_exact extremely short exact matches\n");
                                } else {
                                        push (@answer, "\t\textremely short exact match\n");
                                }
                        }
                        if ($extremely_short_similar > 0) {
                                if ($clock == 1) {push (@answer, "\t");}
                                if ($extremely_short_similar > 1) {
                                        push (@answer, "\t\t$extremely_short_similar similar matches\n");
                                } else {
                                        push (@answer, "\t\textremely short similar match\n");
                                }
                        }
                        if ($extremely_short_distant > 0) {
                                if ($clock == 1) {push (@answer, "\t");}
                                if ($extremely_short_distant > 1) {
                                        push (@answer, "\t\t$extremely_short_distant distant matches\n");
                                } else {
                                        push (@answer, "\t\textremely short distant match\n");
                                }
                        }


			$answer = join ('', @answer); 
			return $answer;
				
}


sub find_length {
		my $line = $_[0]; my $max = $_[1];
		my @line = (); my @u = (); my $answer = "";
		my $length = (); my $value = ();
		@line = split (/\//, $line);
		@u = split (/ /, $line[1]);
		$length = $u[0];
		$value = ($length/$max) * 100;
		if ($value < 80) {
			if ($value < 20) {
				$answer = "extremely short";
			} else {
				$answer = "short";
			}
		}
	return $answer;
}
sub really_find_length {
		my $line = $_[0];

                my @line = (); my @u = (); my $answer = "";
		my $length = ();
                @line = split (/\//, $line);
                @u = split (/ /, $line[1]);
                $length = $u[0];

		return $length;
}


sub find_identity {
			$line = $_[0]; @line = ();

			@line = split (/\,/, $line);
			@id_area = split (/ /, $line[0]);
			$value = $id_area[$#id_area];

   			$value =~s/\(//g; $value =~s/\)//g; $value=~s/\%//g;

			return $value;
}




sub capture_name {
			my $line = $_[0];
			my @line = (); my $name_area = (); my @name_area = ();
			my @name = (); my $name = ();
			if ($line =~/\|/) {
			@line = split (/\|/, $line);
			$name_area = $line[2]; $name_area = remove_terminal_spaces ($name_area);
			@name_area = split (/ /, $name_area);
			} else {
			chomp $line;
			@line = split (/\//, $line);
			$name_area = $line[1]; @name_area = split (/ /, $name_area);
			shift @name_area;
			}
			push (@name, $name_area[0], $name_area[1]);
			$name = join (" ", @name);

		return $name;
}

sub remove_terminal_spaces {
				my $subject = $_[0];
				my @subject = (); my $next = 0;

				@subject = split ('', $subject);
				
				do {
					if ($subject[0] =~/ /) {
						shift @subject;
					} else {
						$next = 1;
					}
				} until ($next == 1);
				$subject = join ('', @subject);
 			return $subject;
}


sub make_blast_databases {
                          my $fasta_file = $_[0];
                          my $datatype = $_[1];
                          my @argv = ();
                          my $d = ();
  			  my $fat = "$fasta_file";
                          if ($datatype =~/DNA/) {
                                            $d = "F";
                          }
                          if ($datatype =~/PROT/) {
                                            $d = "T";
                          }

                          @argv = ("cp $fat BLAST");
                          system(@argv) == 0 or die "Cannot move fasta file to BLAST directory\n";

                          @argv = ();
                          @argv = ("BLAST/bin/formatdb -i BLAST/$fasta_file -p $d -o T");
                          system(@argv) == 0 or die "Could not make blast databases\n";

}

sub blast_off {
                    my $query_file = $_[0];
                    my $database = $_[1];
                    my $datatype = $_[2];
                    my $d = ();
                       if ($datatype =~/DNA/) {
                                       $d = "blastn";
                        }
                       if ($datatype =~/PROT/) {
                                       $d = "blastp";
                        }
                    my $output = $query_file;
                    $output =~s/\.fasta//g;
                    my @output = ();
                    @output = split ('', $output);
                    my $suffix = "\.out";
                    push (@output, $suffix);
                    $output = join ('', @output);


 #                   my @argv = ();
 #                   @argv = ("mv $query_file BLAST");
 #                   system(@argv) == 0 or die "Could not move $query_file to BLAST directory\n";

chdir 'BLAST';
                    @argv = ();
                    @argv = ("bin/blastall -p $d -d $database -i ../$query_file -o ../$output -a 2");
                    #print '$output ='. $output;
                    #print "BLAST $query_file against $database using 2 processors\n";
                    system(@argv) == 0; # or die "Could not blast $query_file\n";
                    # if system call fails probably means that primer is dust, so should feed this back

#                    @argv = ();
#                    @argv = ("mv BLAST/$query_file ./");
#                    system(@argv) == 0 or die "Could not put $query_file back\n";

chdir '..';
}


