#!/usr/bin/perl

#Written by Lucas B. Chaney
#2004 March 7 - April 13
#Modified by LBC & Adam D. Ewing
#v 1.02.04.08.24 - 2004 August 24
#	Reverted to RemoteBlast (for this version only)
#	Bug-fixes in primer design section
#v 1.01.04.07.09 - 2004 July 9
#	Began UNIX port
#	Began conversion to LocalBlast
#	Improved primer design somewhat
#	Fixed "Parsing for unique pairs.../Designing primers..." debug statements
#	Implement dual direction pair parsing
#	Revised naming convention to more standard 3'/5' as opposed to +/-
#v 1.00.04.06.21 - 2004 June 21
#	Added missing support for "r" argument
#	Modified date convention
#	Implemented version naming convention

#BOX-C Gap Closure
#(B)LAST (O)rganism (X)cross (-) (C)omparison
#############################################################################
#	Purpose:	This program is intended to facilitate genome sequencing by	#
#			suggesting contig alignment in whole genome shotgun sequencing	#
#			projects.														#
#																			#
#	Method:		Contig file will be accepted and processed, an array of		#
#			both terminal sequences of all contigs will be generated. This	#
#			array will be compared against known sequences through BLAST to	#
#			find homologous	sequences to suggest alignments. The sequence	#
#			of the contigs in these putative alignments will then be used	#
#			to design PCR primers which can be used to experimentaly		#
#			validate the alignment through the generation of a product,		#
#			which could then be sequenced in order to close the gap.		#
#																			#
#	Variables:	In variable decleration, true scalars follow a '$firstNext'	#
#			capitilization convention or '$i' (one letter) when used in a	#
#			for loop (or other simple placeholding function), arrays follow	#
#			a '@first_next' or '@plurals' naming style, and hashes vary in	#
#			convention.														#
#############################################################################

#---Invoke packages/protocols in use-----------------------------------------
use warnings;
use Bio::Perl;
use Bio::SearchIO;
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
#use Bio::Tools::Blast::Run::LocalBlast;
use Bio::Tools::Run::RemoteBlast; 
#</LBC>
use Bio::Search::Hit::BlastHit;
use Bio::Search::Result::BlastResult;
use Bio::Search::Result::GenericResult;
use Bio::Tools::Run::Primer3;
use POSIX;	#for floor();
use Switch;
#############################################################################

#---Control variables--------------------------------------------------------
#::::Initialize with default values::::::::::::::::::::::::::::::::::::::::::
$endSz = 500;				#Size of contig ends to BLAST
$exclSz = 0;				#Size of contig ends to exclude
$cutOffSmEnds = "";		#Ignore contigs not twice as long as both above
$outFname = "box-c";		#Filename for all files generated by program
$outFtype = "fasta";		#File extension and type of files <see below>
$dBase = 'nr';				#Database for BLAST search
$eVal = '1e-10';			#eVal for BLAST search
$v = "true";				#verbosity aka debug info
$gapSz = 3000;				#maximum size of putative gap
$overlapSz = 0;				#number of shared bases defining an overlap
$waitTime = 15;				#time to wait between requests for RIDs
$toWrite = "ps";				#What types of files to write:
							#	e	|	-e-nds
							#	b	|	-B-LAST results
							#	m	|	names of -M-atches
							#	p	|	-p-rimers
							#	s	|	-s-ummary of results
$numBfiles = 0;				#Number of BLAST files from previous runs
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
$p3path = "/usr/bin/primer3";
#@minProdIncl;
#@maxProdIncl;
$minProdIncl[0] = $endSz * 0.6;
$maxProdIncl[0] = $endSz * 0.8;
$maxPrimerTries = 3;
#</LBC>
#############################################################################

#---Validate input-----------------------------------------------------------
#::::Ensure arguments present::::::::::::::::::::::::::::::::::::::::::::::::
defined($ARGV[0]) 
	or die "No arguments given. Please call with \"-h\" for help\n";
#::::Process arguments:::::::::::::::::::::::::::::::::::::::::::::::::::::::
$hasharg = {};
for ($argNum = 0; $argNum < scalar(@ARGV); $argNum++) {
	$part = $ARGV[$argNum];
	substr($part, 0, 1) eq "-" 
		or die "Improper argument format. Use arg \"-h\" for help\n";
	$argKey = substr($part, 1, 1);
	if(length($part) > 2 && substr($part, 2) eq ":"){
		$argVal = $ARGV[++$argNum];
		while($argNum+1 < scalar(@ARGV) &&
			substr($ARGV[$argNum+1], 0, 1) ne "-"){
			$argVal = $argVal . " " . $ARGV[++$argNum];
		}
		if($argKey ne "i"){	$hasharg->{$argKey} = $argVal;	}
		elsif(defined($hasharg->{i})){	push(@{ $hasharg->{i} }, $argVal);	}
		else{	$hasharg->{i} = [ $argVal ];	}
	} elsif (length($part) == 2){
		$hasharg->{$argKey} = 1;
	}else {die "Improper argument format. Use arg \"-h\" for help\n";}
}
#::::Respond to arguments::::::::::::::::::::::::::::::::::::::::::::::::::::
if(defined($hasharg->{h})){
	print "The general format of arguments is \"-t: 1e-6 -v -i: file.ext\".\n\n";
	print "Explanation of rules:\n";
	print "1.\tEach argument begins with a \"-\" as an identifier.\n";
	print "2.\tA single lower-case letter follows the \"-\". (See below.)\n";
	print "3.\tArguments with values contain \": \" right after the letter.\n";
	print "4.\tValues for arguments do not begin with \"-\".\n";
	print "5.\tValues to turn on boolean arguments should be ommitted.\n";
	print "6.\tPaths and filenames contain no more than 1 consecutive space.\n";
	print "7.\tThe program will (should) ignore unknown arguments.\n";
	print "8.\tOptional arguments are identified [optional], below.\n";
	print "\n::::Argument identifiers::::::::::::::::::::::::::::::::::\n";
	print "id|values taken |default        |more info\n";
	print "--+-------------+---------------+-------------------------\n";
	print "h |boolean\t|false\t\t|displays this message\n";
	print "v |boolean\t|$v\t\t|turns on or off verbose debugging\n";
	print "i |[path]fname\t|none\t\t|input filename, MUST BE PROVIDED!\n";
	print "e |integer\t|$endSz\t\t|the size in bp of the contig ends\n";
	print "y |integer\t|$exclSz\t\t|the size of contig ends to exclude\n";
	print "g |integer\t|$gapSz\t\t|the maximum allowable gap size\n";
	print "l |integer\t|$overlapSz\t\t|the maximum allowable overlap size\n";
	print "d |nr,est,etc.\t|$dBase\t\t|the BLAST database\n";
	print "t |1e-10,etc.\t|$eVal\t\t|threshold for BLAST hits\n";
	print "f |e,b,m,p,s\t|$toWrite\t\t|files ouput:(e)nds, (b)last, (m)atches,\n";
	print "\t\t\t\t|\t\t(p)rimers, (s)ummary\n";
	print "o |[path]fname\t|$outFname.$outFtype\t|where to output files\n";
	print "w |integer >= 5\t|$waitTime\t\t|delay in seconds between requests\n";
	print "r |boolean\t|$cutOffSmEnds\t\t|remove contigs if ends overlap\n";
	if(!defined($hasharg->{i})){die "Thanks for viewing help\n";}
}
if(defined($hasharg->{i})){ @infiles	= @{$hasharg->{i}};	}
	else {die "Input file not specified. Use \"-h\" to see more info\n";}
if(defined($hasharg->{v})){ $v			= $hasharg->{v};	}
#<LBC edit_ver="1.00.04.06.21" edit_date="2004-06-21">
if(defined($hasharg->{r})){ $cutOffSmEnds	= $hasharg->{r};	}
#</LBC>
if(defined($hasharg->{e})){ $endSz		= $hasharg->{e};	}
if(defined($hasharg->{y})){ $exclSz		= $hasharg->{y};	}
if(defined($hasharg->{g})){ $gapSz		= $hasharg->{g};	}
if(defined($hasharg->{l})){ $overlapSz	= $hasharg->{l};	}
if(defined($hasharg->{d})){ $dBase		= $hasharg->{d};	}
if(defined($hasharg->{t})){ $eVal		= $hasharg->{t};	}
if(defined($hasharg->{f})){ $toWrite	= $hasharg->{f};	}
if(defined($hasharg->{w})){ $waitTime	= $hasharg->{w} > 5 ? $hasharg->{w} : 5;	}
if(defined($hasharg->{o})){
	$outFile	=	$hasharg->{o};	
	$outFile =~ /([^\.]+)\.(\S+)/;
	$outFname = $1;
	$outFtype = $2;
}
$infiles[0] =~ /[^\.]+\.(\S+)/;
$ext = lc($1);
$filetype = "unknown";
switch($ext){
	case "fasta"	{$filetype = "fasta";}
	case "fast"		{$filetype = "fasta";}
	case "seq"		{$filetype = "fasta";}
	case "fa"		{$filetype = "fasta";}
	case "fsa"		{$filetype = "fasta";}
	case "nt"		{$filetype = "fasta";}
	case "aa"		{$filetype = "fasta";}
	case "gb"		{$filetype = "genbank";}
	case "gbank"	{$filetype = "genbank";}
	case "genbank"	{$filetype = "genbank";}
	case "gbs"		{$filetype = "genbank";}
	case "gbk"		{$filetype = "genbank";}
	case "scf"		{$filetype = "scf";}
	case "pir"		{$filetype = "pir";}
	case "embl"		{$filetype = "embl";}
	case "ebl"		{$filetype = "embl";}
	case "emb"		{$filetype = "embl";}
	case "dat"		{$filetype = "embl";}
	case "txt"		{$filetype = "raw";}
	case "gcg"		{$filetype = "ace";}
	case "bsm"		{$filetype = "bsml";}
	case "bsml"		{$filetype = "bsml";}
	case "swiss"	{$filetype = "swiss";}
	case "sp"		{$filetype = "swiss";}
	case "phd"		{$filetype = "phd";}
	case "phred"	{$filetype = "phd";}
	case "fastq"	{$filetype = "fastq";}
	case "exp"		{$filetype = "exp";}
}
if($filetype eq "unknown"){
	die "File of unknown type provided.\n";
}
#############################################################################

#---Acknowledge previous BLAST results (if needed)---------------------------
if(scalar(@infiles) > 1){
	if($v){print "Previously output BLAST results detected...\n";}
	$blastFile = $infiles[1];
	$blastFile =~ /(\d+)(.*)/;
	$numBfiles = $1;
	$blastFile = $2;
}
#############################################################################

#---Read in contig file-------------<thank you Bioperl!!>--------------------
if($v){print "Reading file $infiles[0]...\n";}
@all_contigs = read_all_sequences($infiles[0],$filetype);
if($v){print "File read sucessfully.\n";}
#############################################################################

#---Prepare ends-------------------------------------------------------------
if($v){print "Preparing ends of ", scalar(@all_contigs), " contigs...\n";}
$tooSmall = 0;
foreach $contig (@all_contigs){
	if(length($contig->seq()) < ($endSz + $exclSz)){$tooSmall++;}
	else{
		$seq = uc($contig->seq());
		$lEnd = Bio::Seq->new(
			-seq              => substr($seq, $exclSz, $endSz),
			-display_id       => $contig->display_id() . "-L",
			-accession_number => $contig->accession_number(),
			-alphabet         => $contig->alphabet() );
		$rEnd = Bio::Seq->new(
			-seq              => substr($seq, -($endSz + $exclSz), $endSz),
			-display_id       => $contig->display_id() . "-R",
			-accession_number => $contig->accession_number(),
			-alphabet         => $contig->alphabet() );
		if($cutOffSmEnds){
			if(length($contig->seq()) < (2 * ($endSz + $exclSz))){$tooSmall++}
			else{
				push(@L_ends, $lEnd);
				push(@R_ends, $rEnd);
			}
		}
		else{
			push(@L_ends, $lEnd);
			push(@R_ends, $rEnd);
		}
	}
}
push(@all_ends, @R_ends);
push(@all_ends, @L_ends);
if($v){
	$goodEnds = (scalar(@all_contigs) - $tooSmall) * 2;
	print "$tooSmall contigs were too small, $goodEnds ends prepared.\n";
}
#::::Prepare file of ends::::::::::::::::::::::::::::::::::::::::::::::::::::
if($toWrite =~ /e/i){
	if(!$outFile){$outFile = $outFname . '.' . $outFtype;}
	if($v){print "Writing to ends-$outFile...\n";}
	write_sequence(">ends-$outFile",$outFtype,@all_ends);
}
#############################################################################

#---BLAST ends---------------------------------------------------------------
#::::Configure BLAST:::::::::::::::::::::::::::::::::::::::::::::::::::::::::
if($v){print "Preparing to blast...\n";}
$remote_blast = Bio::Tools::Run::RemoteBlast->new(
	'-prog'			=> 'blastn',
	'-data'			=> $dBase,
	'-expect'		=> $eVal );
if($v){print "Submitting blast...\n";}
#::::Run BLAST on ends:::::::::::::::::::::::::::::::::::::::::::::::::::::::
#Starting at 0 and going until you have gotten good results on all ends
for($lastGood = $numBfiles; $lastGood < $#all_ends;){
	$beenError = "false";
	if($v){
		$resNum = $lastGood;
		for($t = 1; $t * 10 <= $resNum; $t *= 10){}
		#Format the output to be neat. i.e. 001 of 100, 01 of 10, 1 of 1
		for($numEnds = scalar(@all_ends); 
			$numEnds > 10 * $t; 
			$numEnds = $numEnds / 10){$resNum = "0" . $resNum;}
	}
	#Start with the element after the last good one and go until
	#	you've done all ends or encountered an error
	for($cur = $lastGood; 
		(($cur < scalar(@all_ends)) && ($beenError ne "true")); 
		$cur++) {
		$end = $all_ends[$cur];						#Get the current end
		my $r = $remote_blast->submit_blast($end);	#Blast it
		if($v){
			$resNum++;
			print "Retrieving result $resNum of ", scalar(@all_ends);
		}
		#Try and get the results
		while ( my @rids = $remote_blast->each_rid ) {
			foreach my $rid ( @rids ) {
				my $rc = $remote_blast->retrieve_blast($rid);
				if( !ref($rc) ) {
					if($v){print ".";}
					if( $rc < 0 ) {	#If $rc < 0, there has been an error
						$beenError = "true";
						$remote_blast->remove_rid($rid);
						if($v){print "\nUnable to retrieve $resNum!";}
					}
					sleep $waitTime unless ($beenError eq "true");
				}
				else {
					if($v){print "|\n";}
					if($toWrite =~ /b/i){
						$fName = $cur+1 . "-" . $outFname . ".blast";
						$remote_blast->save_output($fName);
						if($v){print "BLAST results written to $fName\n";}
					}
					$lastGood = $resNum;
					while(my $result = $rc->next_result()){
						push(@blast_results, $result);	#store the output
					}
					$remote_blast->remove_rid($rid);
				}
			}
		}
	}
	if($beenError eq "true"){
		if($v){
			print "\nError in retrieving results\n";
			print "This is most likely do to a timeout of the connection to the\n";
			print "NCBI BLAST server.\n";
			print "Please try the follwoing:\n";
			print "Use \"-w:\" and increase the wait time. 20-30 is recommended\n";
			print "Use \"-f: b\" to output intermediate results. Then, call the\n";
			print "\tprogram with a second \"-i:\" equal to the last BLAST file output.\n";
		}
		die "Sorry for the inconvenience\n";
	}
}
#############################################################################

#---Parse blast output-------------------------------------------------------
#::::Load previous BLAST results (as needed):::::::::::::::::::::::::::::::::
if(scalar(@infiles) > 1){
	if($v){print "Preparing to read $blastFile from 1 to $numBfiles...\n";}
	for($f = 1; $f <= $numBfiles; $f++){
		$file = $f . $blastFile;
		$searchio = new Bio::SearchIO( -file => $file );
		while( my $result = $searchio->next_result ){
			push(@blast_results, $result);
		}
	}
}
#::::Parse old and new BLAST results:::::::::::::::::::::::::::::::::::::::::
if($v){print "Parsing blast results...\n";}
foreach $blastResult (@blast_results) {
	while( my $blastHit = $blastResult->next_hit ){
		$qName = $blastResult->query_name;
		#hsp: highest scoring pair
		while( my $hsp = $blastHit->next_hsp ){
			$curHit = {};
			$curHit->{sbjct} = $blastHit->description;
			$curHit->{query} = $qName;
			$curHit->{strand} = $hsp->strand('sbjct');
			$curHit->{start} = $hsp->start('sbjct');
			$curHit->{end} = $hsp->end('sbjct');
			push(@hits, $curHit);
		}#end while hsp
	}#end while blastHit
}#end foreach
if($v){print "Finished parsing blast results.\n";}
#############################################################################

#---Find matches-------------------------------------------------------------
#::::Prepare hash of hits by subject identifier::::::::::::::::::::::::::::::
if($v){print "Searching for paired matches...\n";}
for($sbjctCount = 0; $sbjctCount < $#hits; $sbjctCount++){
	push(@all_sbjcts, $hits[$sbjctCount]{sbjct});
}
%seenSbjcts = ();
foreach $curSbjct (@all_sbjcts) {
	push(@uniq_sbjcts, $curSbjct) unless $seenSbjcts{$curSbjct}++;
}
$u_sbjcts = {};
foreach $uSbjct (@uniq_sbjcts) {
	foreach $hit (@hits) {
		%hashhit = %{ $hit };
		if($hashhit{sbjct} eq $uSbjct){
				push(@sbjct_matches, $hit);
		}
	}
	if(scalar(@sbjct_matches) > 1){
		#sort the sbjct_matches before storing them
		@sbjct_matches = sort { #min of a and min of b
				($a->{start} < $a->{end} ? $a->{start} : $a->{end}) 
					<=> 
				($b->{start} < $b->{end} ? $b->{start} : $b->{end})
			} @sbjct_matches;
		$u_sbjcts->{$uSbjct} = [@sbjct_matches];
	}
	@sbjct_matches=();
}
#::::Search each subject for hits within range of one another::::::::::::::::
%hashsbjcts = %{ $u_sbjcts };
foreach $key ( keys %hashsbjcts ) {
	@cur_matches = @{ $hashsbjcts{$key} };
	for($matchCount = 0; $matchCount < $#cur_matches - 1; $matchCount++){
		%testor = %{ $cur_matches[$matchCount] };
		$testeeCount = $matchCount+1;
		%testee = %{ $cur_matches[$testeeCount] };
		$ss = $testee{start} - $testor{start};
		$se = $testee{start} - $testor{end};
		$es = $testee{end} - $testor{start};
		$ee = $testee{end} - $testor{end};
		@all_dists = ($ss, $se, $es, $ee);
		@all_dists = sort {$a <=> $b} @all_dists;
		$inRange = $all_dists[0] <= $gapSz;
		while($inRange && $testeeCount < scalar(@cur_matches)){
			$torName = $testor{query};
			$teeName = $testee{query};
			$torEnd = substr($torName, -1);
			$teeEnd = substr($teeName, -1);
			$torName = substr($torName, 0, length($torName)-2);
			$teeName = substr($teeName, 0, length($teeName)-2);
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
			$torStr = $testor{strand};
			$teeStr = $testee{strand};
#</LBC>
			$overlapping = $all_dists[0] + $overlapSz <= 0;
			$nameSame = $torName eq $teeName;
			$endsNstrands = ($torEnd ne $teeEnd || $torStrand ne $teeStrand) &&
				!($torEnd eq $teeEnd && $torStrand eq $teeStrand);
			if(!$nameSame && !$overlapping && $endsNstrands){
				$pair = {};
				$pair->{left} = $torName;
				$pair->{right} = $teeName;
				$pair->{lstr} = $torStr;
				$pair->{rstr} = $teeStr;
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
				$torStrand = $torStr > 0 ? "(3')" : "(5')";
				$teeStrand = $teeStr > 0 ? "(5')" : "(3')";
#<LBC comment="not needed?">	$torStr = $torStr > 0 ? "+" : "-";
#<LBC comment="not needed?">	$teeStr = $teeStr > 0 ? "+" : "-";
				$pair->{revid} = $teeName . $teeStrand . "," . $torName . $torStrand;
				$pair->{forfound} = 1;
				$pair->{revfound} = 0;
#</LBC>
				$pair->{id} = $torName . $torStrand . "," . $teeName . $teeStrand;
				$pair->{hosts} = [ $key ];
				$pair->{gaps} = [ $all_dists[0] ];
				$pair->{avgap} = $all_dists[0];
				push(@all_pairs, $pair);
			}
			if(++$testeeCount < $#cur_matches){
				%testee = %{ $cur_matches[$testeeCount] };
				$ss = $testee{start} - $testor{start};
				$se = $testee{start} - $testor{end};
				$es = $testee{end} - $testor{start};
				$ee = $testee{end} - $testor{end};
				@all_dists = ($ss, $se, $es, $ee);
				@all_dists = sort {$a <=> $b} @all_dists;
				$inRange = $all_dists[0] + $overlapSz <= $gapSz;
			}
		}
	}
}
if($v){print "Finished searching for pairs.\n";}
#############################################################################

#---Find primers-------------------------------------------------------------
#::::Find unique pairs:::::::::::::::::::::::::::::::::::::::::::::::::::::::
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
if($v){print "Parsing for unique pairs....\n";}
#</LBC>
%seenPairs = ();
foreach $curPair (@all_pairs) {
	if($seenPairs{$curPair->{id}}++){
		$i = 0; 
		while($i < $#uniq_pairs && $uniq_pairs[$i++]->{id} ne $curPair->{id}){}
		push(@{ $uniq_pairs[$i]->{hosts} }, ${ $curPair->{hosts} }[0]);
		push(@{ $uniq_pairs[$i]->{gaps} }, ${ $curPair->{gaps} }[0]);
		$uniq_pairs[$i]->{avgap} = (($uniq_pairs[$i]->{avgap} * 
			$#{ $uniq_pairs[$i]->{gaps} }) + ${ $curPair->{gaps} }[0]) 
			/ scalar(@{ $uniq_pairs[$i]->{gaps} });
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
		$uniq_pairs[$i]->{forfound}++;
	} elsif($seenPairs{$curPair->{revid}}++){		$i = 0; 
		while($i < $#uniq_pairs && $uniq_pairs[$i++]->{id} ne $curPair->{id}){}
		push(@{ $uniq_pairs[$i]->{hosts} }, ${ $curPair->{hosts} }[0]);
		push(@{ $uniq_pairs[$i]->{gaps} }, ${ $curPair->{gaps} }[0]);
		$uniq_pairs[$i]->{avgap} = (($uniq_pairs[$i]->{avgap} * 
			$#{ $uniq_pairs[$i]->{gaps} }) + ${ $curPair->{gaps} }[0]) 
			/ scalar(@{ $uniq_pairs[$i]->{gaps} });
		$uniq_pairs[$i]->{revfound}++;
#</LBC>
	} else {
		push(@uniq_pairs, $curPair);
	}
}
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
if($v){print "Finished parsing for unique pairs.";}
#</LBC>
#::::Design primers::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
if($v){print "Designing primers...\n";}
$summary = "";
$primer = "";
foreach $uPair (@uniq_pairs) {
	if($v){
		print $uPair->{id}, ": ", scalar(@{ $uPair->{hosts} });
		print "  |  ", $uPair->{avgap}, "\n";
	}
	if($toWrite =~ /s/i){
		$summary .= $uPair->{id} . ": " . scalar(@{$uPair->{hosts}});
		$summary .= "  |  " . $uPair->{avgap} . "\n";
		for($i = 0; $i <= $#{ $uPair->{hosts} }; $i++){
			$summary .= $i+1 . ": h:\t" . substr(${$uPair->{hosts}}[$i],0,50);
			$summary .= "\tg:   " . ${$uPair->{gaps}}[$i] . "\n";
		}
	}
	$l = -1;
	while($all_contigs[++$l]->display_id ne $uPair->{left}){}
	if($uPair->{lstr} > 0){$uPair->{lcont} = $all_contigs[$l];}
	else{
		$seq = $all_contigs[$l]->seq;
		$seq =~ tr/ACGT/TGCA/;
		$seq = reverse $seq;
		$uPair->{lcont} = Bio::Seq->new(
			-seq              => $seq,
			-display_id       => $all_contigs[$l]->display_id(),
			-accession_number => $all_contigs[$l]->accession_number(),
			-alphabet         => $all_contigs[$l]->alphabet() );
	}
	$r = -1;
	while($all_contigs[++$r]->display_id ne $uPair->{right}){}
	if($uPair->{rstr} > 0){$uPair->{rcont} = $all_contigs[$r];}
	else{
		$seq = $all_contigs[$r]->seq;
		$seq =~ tr/ACGT/TGCA/;
		$seq = reverse $seq;
		$uPair->{rcont} = Bio::Seq->new(
			-seq              => $seq,
			-display_id       => $all_contigs[$r]->display_id(),
			-accession_number => $all_contigs[$r]->accession_number(),
			-alphabet         => $all_contigs[$r]->alphabet() );
	}
	$lSeq = $uPair->{lcont}->seq();
	$rSeq = $uPair->{rcont}->seq();
	$nSeq = 'N' x POSIX::floor( $uPair->{avgap} );
	$uPair->{p3seq} = Bio::Seq->new(
									-seq		=>	$lSeq . $nSeq . $rSeq ,
									-display_id	=>	$uPair->{id} ,
									-alphabet	=>	'dna' );
	push(@out_pairs, $uPair->{p3seq});
	my $primer3 = Bio::Tools::Run::Primer3->new(-seq=>$uPair->{p3seq});
	$primer3->program_dir($p3path);
#<LBC edit_ver="1.01.04.07.09" edit_date="2004-07-09">
	for($i = 1; $i <= $maxPrimerTries; $i++){
		$minProdIncl[$i] = POSIX::floor($minProdIncl[($i-1)] / 2);
		$maxProdIncl[$i] = $maxProdIncl[($i-1)] + ($maxProdIncl[0] - $minProdIncl[0]);
	}
	$nPT = 1;	#number Primer Tries
#<LBC edit_date="2004-07-28">
	$mPR = 0;	#max Primers Returned
	while(($mPR == 0) && ($nPT <= $maxPrimerTries)){
		$inclSz = ($maxProdIncl[$nPT] - $minProdIncl[$nPT] + 1);
		$inclReg = length($lSeq) - (1 + $maxProdIncl[$nPT]);
		$inclLen = length($lSeq) + length($nSeq) + $minProdIncl[$nPT] + $inclSz - $inclReg;
		$inclReg .= "," . $inclLen;
		$primer3->add_targets(INCLUDED_REGION=> $inclReg);
		$target = length($lSeq) - $minProdIncl[$nPT] + 1;
		$targetSz = length($lSeq) + length($nSeq) + $minProdIncl[$nPT] - $target;
		$target .= "," . $targetSz;
		$primer3->add_targets(TARGET=>$target);
		$szRange = length($nSeq) + (2 * $minProdIncl[$nPT]) - 2;
		if($szRange < 50){$szRange = 50;} #should be Primer size, but 50 works for now
#<ADE 7-27-04 ... order of operations>
		$szRange .= "-" . (length($nSeq) + (2 * $maxProdIncl[$nPT]) - 2);
		#$szRange = $szRange . "-" . length($nSeq) + (2 * $maxProdIncl[$nPT]);
#</ADE>
		$primer3->add_targets(PRIMER_PRODUCT_SIZE_RANGE=>$szRange);
		#print "$inclReg | $inclReg2 | $target |--> $szRange <--|";
		#print "$inclReg | $target |--> $szRange <--|";
		$uPair->{primers} = $primer3->run;
		$mPR = $uPair->{primers}->{maximum_primers_returned};
		#print "$mPR\n";
		$nPT++;
	}
	$primer .= ">" . $uPair->{id} . "\n";
	if(defined($uPair->{primers}->{results}->{PRIMER_ERROR})){
		#$primer .= "Error while attempting to generate primers.\n"
		$primer .= $uPair->{primers}->{results}->{PRIMER_ERROR} . "\n";
		#print $uPair->{primers}->{results}->{PRIMER_ERROR} . "\n";
	} else {
#</LBC edit_date="2004-07-28" edit_date="2004-07-09">
		$numPrimers = $uPair->{primers}->{maximum_primers_returned};
		$primer .= $numPrimers . " primers were generated.\n";
		for($i = 0; $i < $numPrimers; $i++){
			$primer .= "Right $i: ";
			$primer .= $uPair->{primers}->{results_by_number}->{$i}->{PRIMER_RIGHT_SEQUENCE} . "\n";
			$primer .= "Left  $i: ";
			$primer .= $uPair->{primers}->{results_by_number}->{$i}->{PRIMER_LEFT_SEQUENCE} . "\n";
		}
	}
}
if($v){print "Primers designed.\n";}
if($toWrite =~ /s/i){
	if($v){print "Writing summary to summary-$outFname.txt...\n";}
	open(SUMMARY, ">summary-$outFname.txt");
	print SUMMARY $summary;
}
if($toWrite =~ /p/i){
	if($v){print "Writing primers to primers-$outFname.txt...\n";}
	open(PRIMERS, ">primers-$outFname.txt");
	print PRIMERS $primer;
}
#::::Prepare file of matches:::::::::::::::::::::::::::::::::::::::::::::::::
if($toWrite =~ /m/i){
	if(!$outFile){$outFile = $outFname . '.' . $outFtype;}
	if($v){print "Writing to matches-$outFile...\n";}
	
	write_sequence(">matches-$outFile",$outFtype,@out_pairs);
}
#############################################################################