#!/usr/bin/perl -w

#  script to convert blastall output to blink input.  

$version = "v0.37 30 Oct 2006";	
# updated user interface
# fixed bug:  earlier versions did not work properly if phi > 0 and sigma = 0.  


$Usage ="
Usage:
		
./blast2blink.pl 
	-i infile  
		11 columns, achieved using blast's -m 8 option
		first column is unique id OR if it contains (gi|### OR gi###) then gi### will be used.
			program will attempt to read format.
		beginning and end positions of the hit are in the 7th, 8th, 9th, 10th fields.
		NOTE: assumes that all hits between two gis are sequential in blast ouput -- a time saver,
			but assumption should be checked
	-o outfile
		two-column file of ids that were accepted by all criteria, ready for blink.
			

    optional commands:  
	-p phi_threshhold
		0 to 1
		for a pair of sequences to be considered linked, they must have this 
		  minimum number of bases spanned by all hits between the two sequences, 
		  scaled to the total number of bases spanned by both sequences, including unaligned \"sticky ends\".
		default = 0.
	-s2 two-way_sigma_treshhold 
		0 to 1
		for a pair of seuences to be considered linked, they must have this 
		  minimum proportion of each sequence that hits to the other sequence.
		default = 0.
	-s1 one-way_sigma_threshhold
		0 to 1
		for a pair of sequences to be considered linked, they must have this 
		  minimum proportion of one of the sequences that hits to the other sequence.
		  (the other sequence may hit to the first with a lower proportion of its bases).
		default = 0.   
	-t table_file
		if values are given for sigma or phi, this file is required.
		format 1 (program will try to detect format):
		  fields must be delimited by bars |.  
		  first field must contain the gi in gixxx format.
		  third field must contain the length of the sequence.
		format 2:
		  fields delimited by space
		  first field contains unique id
		  second field contains length of sequence
			
	-g filename
		file that contains a list of gis.  
		if specified, the output will contain only those edges involving gis on this list.
		  gis must be gixxx or gi|xxxx and can be in any column but only one per line. 

	Calculation of phi and sigma:
					a 	  b1   	c1   b2   c2  ... bn     cn   d		    
					-----------------------------------------------
						  |||||||    |||\\\\       \\\||||
					     ----------------------------------------------------
					     w	  x1	y1   x2	    y2	    xn   yn		z


				phi = 	max ((cn-b1),(yn-x1)) /
					( max((b1-a),(x1-w)) + max ((cn-b1),(yn-x1)) + max((z-yn),(d-cn)) )

				sigma1 = (sum(ci-bi+1))/d), i = 1 to n 	
				sigma2 = (sum(yi-xi+1))/z), i = 1 to n	

			NOTE: this description is not quite correct:  the sums are actually set unions, because hits can overlap.  

	examples:

		./blast2blink.pl -i paps.blast.out -o paps.hitsp90s90 -t paps.table -s 0.9 -p 0.9
			makes a hitlist where each edge consists of sequences that are similar enough to hit with at least 90% 
			  their bases and similar enough in length that the hitting region is at least 90% of the total 
			  region spanned.
		
		./blast2blink.pl -i paps.blast.out -o paps.sampleHitsp50s20 -t paps.table -s 0.2 -p 0.5 -g paps.sampleGis
			makes a hitlist that involve only the gis in the sample file, for which 20% or more of each seq hits
			  to the other, and where the hitting region spans at least 50% of the total region spanned. 

		./blast2blink.pl -i paps.blast.out -o paps.hitsp0s90 -m 1 -s 0.9
			makes a hitlist where at least 90% of at least one of the sequences hits to the other, and no restrictions
			  are placed on where those hits are in the total region spanned.

#  	NOTE: skips the edge if one strand is reversed relative to the other

";

# command line input ----------------------------------
while ($fl = shift @ARGV)
  {
  $par = shift @ARGV;
  die ("$Usage") if (! $par);
  if ($fl =~ /-i/) {$infile = $par;}
  if ($fl =~ /-o/) {$outfile = $par;}
  if ($fl =~ /-t/) {$tablefile = $par;}
  if ($fl =~ /-g/) {$gifile = $par;}
  if ($fl =~ /-p/) {$phiT = $par;}
  if ($fl =~ /-s1/) {$sigmaOneT = $par;}
  if ($fl =~ /-s2/) {$sigmaTwoT = $par;}
  if ($fl =~ /^-s$/) {$oldCmd = $par;}
  }

# comment out if not runing interactively:

  if (-e $outfile) {print "Output file already exists. Continue?\n"; $ans = (<STDIN>); die if (! ($ans =~ /y/i));}

die ("Use of the command \"-s\" is no longer supported.\n  Please use \"-s1\" (one-way sigma) or \"-s2\" (two-way sigma).\n\n")
	if ($oldCmd);
die ("$Usage") if (! ($infile && $outfile));
die ("Missing infile $infile.\n") if (! -e $infile);
die ("Missing table file.\n") if ($tablefile && (! -e $tablefile));
die ("Missing gi file.\n") if ($gifile && (! -e $gifile));
if ($sigmaOneT || $sigmaTwoT || $phiT) {die ("Must supply table if using sigma and phi\n") if (! $tablefile);}
if (! $sigmaOneT) {$sigmaOneT = 0;}
if (! $sigmaTwoT) {$sigmaTwoT = 0;}
if (! $phiT) {$phiT = 0;}
die ("Please choose -s1 or -s2 (one-way or two-way sigma threshhold), not both.\n") if ($sigmaTwoT && $sigmaOneT);
if ($sigmaTwoT) {$mode=2; $sigmaT=$sigmaTwoT;}
elsif ($sigmaOneT) {$mode=1; $sigmaT=$sigmaOneT;}
else {$mode=0; $sigmaT=0;}

$v = 0;			# turn these flags on for debugging
$v1 = 0;


# print version and settings ----------------------------

print "\nblast2blink.pl\n$version\nM. McMahon\n\n";
print "Settings:\n";
print "  Blast results file:\t$infile\n";
($bform, $exId) = blastFormat($infile);
print "    format detected, the following is an example sequence id:  $exId\n";
if ($gifile) {print "    only sequences referenced in the file $gifile will be considered\n";}
print "  blast2blink.pl output:\t$outfile\n";
print "  sigma threshhold:\t$sigmaT\n";
if ($mode) 
  {
  if ($mode == 1) {print "    at least one sequence will be required to have this proportion hit to the other sequence\n";}
  elsif ($mode  == 2) {print "    both sequences will be required to have this proportion involved in the hit\n";}
  }
print "  phi threshhold:\t$phiT\n";
if ($sigmaT || $phiT) 
  {
  print "  length information for calculating sigma and/or phi will be taken from: $tablefile\n";
  $tform = tableFormat($tablefile);
  if ($tform == 1) {print "    table file appears to have format 1 (bar-delimited, gi required, third field contains length).\n\n";}
  else {print "    table file appears to have format 2 (two columns: unique_id length).\n\n";}
  }


# make length hash if phi or sigma  ---------------------
if ($phiT || $sigmaT)
  {
  open FH3, "<$tablefile";
  while (<FH3>)
    {
    next if (/^\s*$/);
    if ($tform == 1)
      {
      @fields = split /\|/;
      if ($fields[0] =~ /(gi\d+)/) {$gi = $1;}
      if ($fields[2] =~ /(\d+)/) {$giLengthH{$gi} = $1;}
      else {die ("Problem with format of table file (column 3 is not a number).\n");}
      }
    elsif ($tform == 2)
      {
      @f = split;
      $gi = shift @f;
      $len = shift @f;
      die ("Problem with format of table file (column 2 is not a number).\n") if (! ($len =~ /^\d+$/));
      $giLengthH{$gi} = $len;
      }
    }
  close FH3;
  }

# read gi file if necessary ----------------------------
if ($gifile)
  {
  $targs = 1;
  open FHg, "<$gifile";
  while (<FHg>)
    {
    ($gi) = /gi\|?(\d+)/;
    $gi = "gi$gi";
    if (! exists $trgGi{$gi}) {$trgGi{$gi} = 1;}
    }
  close FHg;
  }
else {$targs = 0;}

# parse blast output -----------------------------------
open FH1, "<$infile";
open FH2, ">$outfile";
if (! ($sigmaT || $phiT))
  {
  while (<FH1>)
    {
    @d = split;
    ($def1, $def2, $b1, $e1, $b2, $e2) = ($d[0], $d[1], $d[6], $d[7], $d[8], $d[9]);  
    if (($b1>$e1) || ($b2>$e2)) {print "Reversal in hit $def1 $def2 (skipping hit)\n"; next;}
    if (($bform eq "fst") || ($bform eq "notFstGi"))    
      {
      ($gi1) = ($def1 =~ /gi\|?(\d+)\|?/);
      ($gi2) = ($def2 =~ /gi\|?(\d+)\|?/);
      $gi1 = "gi$gi1";
      $gi2 = "gi$gi2";
      }
    else {$gi1 = $def1; $gi2 = $def2;}
    if ($targs) {next if ( (! exists $trgGi{$gi1}) || (! exists $trgGi{$gi2}) );}
    print FH2 "$gi1\t$gi2\n";
    }
  }
else  				
  {
  $oldGi1=""; $oldGi2="";
  while (<FH1>)
    {
    @d = split;
    if ($v1) {print "Considering line:\n$_";}
    ($def1, $def2, $b1, $e1, $b2, $e2) = ($d[0], $d[1], $d[6], $d[7], $d[8], $d[9]);  
    if (($b1>$e1) || ($b2>$e2)) {print "Reversal in hit $def1 $def2 (skipping hit)\n"; next;}
    if (($bform eq "fst") || ($bform eq "notFstGi"))
      {
      ($newGi1) = ($def1 =~ /gi\|?(\d+)\|?/);
      ($newGi2) = ($def2 =~ /gi\|?(\d+)\|?/);
      $newGi1 = "gi$newGi1"; $newGi2 = "gi$newGi2";
      }
    else {$newGi1 = $def1; $newGi2 = $def2;}
    if (($newGi1 eq $oldGi1) && ($newGi2 eq $oldGi2)) 		# assumes that all hits between two gis are sequential in blast ouput
      {
      # add to existing record if same as before ------
      push @{$r1ref}, [ $b1, $e1 ];
      push @{$r2ref}, [ $b2, $e2 ];
      }
    else
      {
      # process old record, if exists -----------------
      if ($oldGi1 =~ /\w/)
        {
        if ($v1) {print "\ncalculating stats for pair $oldGi1 and $oldGi2:\n";} 
        ($sigma1, $sigma2, $phi) =  processRecord ($oldGi1, $oldGi2, $r1ref, $r2ref);
	$accept = acceptRecord ($sigma1, $sigma2, $phi);     
        if ($accept) {print FH2 "$oldGi1\t$oldGi2\n";}
        if ($v) {print "the edge has acceptance:$accept\tbecause:  $sigma1, $sigma2, $phi\n";}
        if ($v) {print "Continue?\n"; $answer = (<STDIN>); die if ($answer =~ /n/i);}
        $oldGi1 = ""; $oldGi2 = "";
	}
      # start new record -------------------------------
      if (  ($targs && (exists $trgGi{$newGi1}) && (exists $trgGi{$newGi2}) ) || (! $targs)  )
        {
        $r1ref = [ [ $b1, $e1 ] ];
        $r2ref = [ [ $b2, $e2 ] ];
        $oldGi1=$newGi1; $oldGi2=$newGi2;
        }
      } # end if new edge
    } # end while reading lines, process last set ------

    if (  ($targs && (exists $trgGi{$newGi1}) && (exists $trgGi{$newGi2}) ) || (! $targs)  )
      {
      if ($v1) {print "calculating stats for last pair, $oldGi1 and $oldGi2:\n";}
      ($sigma1, $sigma2, $phi) =  processRecord ($oldGi1, $oldGi2, $r1ref, $r2ref);
      $accept = acceptRecord ($sigma1, $sigma2, $phi);
      if ($accept) {print FH2 "$oldGi1\t$oldGi2\n";}
      }
  } # end if length requirements phi or sigma
close FH1;
close FH2;

# subroutines --------------------------------------------------------------------------
#  findPhi, findSumAndRange, acceptRecord, processRecord 

# sub processRecord ----------------------------------------------
#	takes array
#		($gi1, $gi2, $r1ref, $r2ref) where rXref points 
#		  to an array of start/stop pairs for the gi.
#	must have access to giLengthH
#	calls sub findSumAndRange and findPhi
#	returns array
#		($sigma1, $sigma2, $phi) 

sub processRecord 
{
my ($gi1, $gi2, $r1ref, $r2ref) = @_;
my $length1 = $giLengthH{$gi1}; 
my $length2 = $giLengthH{$gi2};
die ("Problem:  $gi1 or $gi2 has length 0\nEither the table file is in an incorrect format\n  or the gis are not in the table file.\nDied in subroutine processRecord.\n") if (! ($length1 && $length2));
my ($sum1, $first1, $last1) = findSumAndRange ($r1ref);
my ($sum2, $first2, $last2) = findSumAndRange ($r2ref);
my $sigma1 = $sum1/$length1; 
my $sigma2 = $sum2/$length2;
if ($v1) {print "sigma1: $sigma1 = $sum1/$length1\nsigma2: $sigma2 = $sum2/$length2\n";}
my $phi = findPhi ($length1, $first1, $last1, $length2, $first2, $last2);
return ($sigma1, $sigma2, $phi);
}

# subroutine findPhi ---------------------------------------------
#	takes array of edge data:
#		($length2, $first1, $last1, $length2, $first2, $last2)
#	returns phi
 
sub findPhi
{
my ($length1, $first1, $last1, $length2, $first2, $last2) = @_;
my $l1 = $first1 - 1;
my $l2 = $first2 - 1;
my $r1 = $length1 - $last1;
my $r2 = $length2 - $last2;
my $h1 = $last1 - $first1 + 1;
my $h2 = $last2 - $first2 + 1;
my @lefts = sort { $a <=> $b } ($l1, $l2);
my @hits = sort { $a <=> $b } ($h1, $h2);
my @rights = sort { $a <=> $b } ($r1, $r2);
if ($v1) {print "region array:  ($l1, $l2, $h1, $h2, $r1, $r2)\n";}
if ($v1) {print "sorted info: @lefts, @hits, @rights\n";}
my $phi = ($hits[1]) / ($lefts[1] + $hits[1] + $rights[1]);
if ($v1) {print "phi: $phi = ($hits[1]) / ($lefts[1] + $hits[1] + $rights[1]) \n";} 
return $phi;
}

# subroutine findSumAndRange -------------------------------------
#	takes array ref that points to array of arrays
#		($r1ref)
#	returns array: sum of hits, start of first hit, end of last hit
#		($sum, $first, $last) 

sub findSumAndRange
{
my ($aref) = @_;
my ($sum, $first, $last, $min, $max) = (0, 0, 0, 0, 0);
my $bb = 0; my $eb = 0;
my @refsA = @{$aref};
my $sRefs = ();
my $n = scalar @refsA;

if ($n == 1)
  {
  $sum = ($refsA[0][1] - $refsA[0][0] + 1);
  ($first, $last) = ($refsA[0][0], $refsA[0][1]);
  }
else 
  {
  @sRefs = sort { ${$a}[0]  <=> ${$b}[0] } @refsA;
  ($min, $max) = @{$sRefs[0]};
  $first = $min;
  for $i (1..($n-1))
    {
    if ($sRefs[$i][0] > $max)  			# ranges do not intersect
      {
      $sub = ($max - $min + 1); 
      $sum += $sub;  
      ($min, $max) = @{$sRefs[$i]};
      }
    else
      {
      if ($sRefs[$i][1] > $max)			# range can be expanded
        {$max = $sRefs[$i][1];}
      } 
    }
  $last = $max;
  $sub = ($max - $min + 1);
  $sum += $sub;
  }
return ($sum, $first, $last);
}

# sub acceptRecord -----------------------------------------------
#	takes array
#		($sigma1, $sigma2, $phi)
#	must have access to all cutoffs and mode value.
#	returns 1 if the record passes sigma and phi criteria
#		otherwise returns 0.

sub acceptRecord
{
my ($sigma1, $sigma2, $phi) = @_;
my ($acceptableSum, $acceptablePhi, $pass) = (0, 0, 0);
if ($mode == 2)
  {
  if (($sigma1 >= $sigmaT) && ($sigma2 >= $sigmaT)) 
    {$acceptableSum = 1;}
  else {$acceptableSum = 0;}
  }
elsif ($mode == 1)
  {
  if (($sigma1 >= $sigmaT) || ($sigma2 >= $sigmaT)) 
    {$acceptableSum = 1;}
  else {$acceptableSum = 0;} 
  }
elsif ($sigmaT == 0)
  {$acceptableSum = 1;}
if ($phi >= $phiT)
  {$acceptablePhi = 1;}
else {$acceptablePhi = 0;}
if ($acceptableSum && $acceptablePhi)
  {$pass = 1;}
else {$pass = 0;}
return $pass;
}

# sub tableFormat -------------------------------------------------
#	takes table_file_name
#	returns format 1 or 2

sub tableFormat
{
my ($tfile) = @_;
my $fmt = 0;
open FHtemp, "<$tfile" or die ("Trouble opening table file to detect its format.\n");
while (<FHtemp>)
  {
  if (/\|/) {$fmt=1;}
  else {$fmt=2;}
  last;
  }
close FHtemp;
return $fmt;
}

# sub blastFormat -------------------------------------------------
#	takes blast_output_file_name
#	returns fst (if gi|\d+.+) or notFstGi (if gi\d+.+) or else uniqueId
#	  and an example (first in file)

sub blastFormat
{
my ($bfile) = @_;
my $fmt = 0;
my $example = 0;
my $id = 0;
open FHtemp, "<$bfile" or die ("Trouble opening blast file to detect its format.\n");
while (<FHtemp>)
  {
  my @f = split;;
  $id = shift @f;
  if ($id =~ /^\s*gi\|(\d+)\|/) {$example = "gi$1"; $fmt="fst";}
  elsif ($id =~ /^\s*(gi\d+)/) {$example = $1; $fmt="notFstGi";}
  else {$example = $id; $fmt = "uniqueId";}
  last;
  }
close FHtemp;
return ($fmt, $example);
}
