#!/usr/bin/perl
#tblastn_pipe.pl version 0.0
#author: Tobias Mede & Sonja J. Prohaska
#  date: Friday 13.01.2012

use strict;
use lib "/homes/brauerei/tobias/PKR/INFGamma/Tummelplatz";
use interval;

my $usage= << "JUS";
  usage: tblastn_pipe.pl [OPTIONS] TBLASTN_OUTPUT_FILE(S)
  needs: 
options: -d       debug
         -h       help
JUS

my @options = grep {substr($_, 0, 1) eq '-'} @ARGV;
my @others  = grep {substr($_, 0, 1) =~ /^[^-]/} @ARGV;

##default settings:
my $debug  = 0 ;         #default: no debuging
my $savedir = "./";
my $max_evalue = 0;
##----------

##adjust to your usage:
#segment 1
my $a1= 32;
my $a2= 98;
my $evala= 1.0e-5;
#segment 2
my $b1= 130;
my $b2= 172;
my $evalb= 1.0e-5;
##----------

##setting parameters:
foreach my $option (@options){
  if    ($option eq "-d")          { $debug = 1; }
  elsif ($option eq "-h")          { print $usage }
  else  { die "Invalid command line option: \'$option\'!\n$usage"; }
}

# read in files and prepare data structure ####################################################

foreach my $file (@others) {
  #prepare directory for saving the results
  my ($path,$basename,$ext)=parse_filename($file);
  $savedir= $path . "locus/";
  mkdir($savedir);
  
  if ($debug) {
    print STDERR "opening file \'$file\'\n";
  }
  open (IN,"<$file") or die "Can\'t open file \'$file\'!\n";
  
  my @data;
  my @Cfilt;
  my $vieh=$basename;
  
  while (<IN>)  {	
    if ($_ =~ /^#/) {
      next;
    }
    
    my @l = split;
    # for tblastn output format: 7 
    my $ori = "";
    my $s = "";
    my $e = "";
    
    # extract orientation information
    if ($l[8]<$l[9]){
      $ori = "+";
      $s = $l[8];
      $e = $l[9];
    }
    else{
      $ori = "-";
      $s = $l[9];
      $e = $l[8];
    }
    
    # data: [0]=chrom [1]=e-value [2]=prot_from [3]=prot_to [4]=genomic_from [5]=genomic_to [6]=orientation +|-
    push @data, [$l[1], $l[10], $l[6], $l[7], $s, $e, $ori];
    if ($debug){
      print STDERR "adding data ($l[1], $l[10], $l[6], $l[7], $s, $e, $ori)\n";
    }
  }
  close IN;

# process data structure: first, sort by chromosome; second, sorted by genomic start position 
  if ($debug){
    print STDERR "sorting data...";
  }
  my @sorted = sort { $a->[0] cmp $b->[0] || $a->[4] <=> $b->[4] } @data;
  if ($debug){
    print STDERR "done.\n";
    print_data_struct(\@sorted,\*STDERR);
  }

  my $bound;
  my $zaehl=-1;
  my @C0;
  foreach my $hit (@sorted) {
    $zaehl++;
    my $chrom = $hit->[0];

    # fill all hits of one chromosome into cluster C0
    if (defined($bound)==0 || $bound eq $chrom) {
      push @C0, $hit;
      $bound=$chrom;
      if($zaehl==$#sorted) {
	  goto TAG}
    }
### filter hits according to given constraints ################################################
    else {
    TAG: $bound=$chrom;
      if ($debug) {
      print STDERR "Filtering " . scalar(@C0) . " hits (bound=$bound, chrom=$chrom).\n";
      }
      for (my $i=0;$i<@C0;$i++) { 

 	#calculate coverage of exons with boudnaries a1-a2 and e-value threshold $evala
	my $intersect=interval::intersect($a1,$a2,@C0->[$i][2],@C0->[$i][3]);
	if ($intersect != 0 ) {	
	  my $cover=($intersect->[1]-$intersect->[0]+1)/($a2-$a1+1);
	  my $cover2=($intersect->[1]-$intersect->[0]+1)/(@C0->[$i][3]-@C0->[$i][2]+1);
	  if ($cover >= 0.9 && $cover2 >=0.9) {
	  
	    if (@C0->[$i][1] <= $evala ) {
	      push @Cfilt, [$C0[$i][0], $C0[$i][1],$C0[$i][2],$C0[$i][3],$C0[$i][4],$C0[$i][5],$C0[$i][6]];
	    }
	  }
	}
	#calculate coverage of exons with boudnaries b1-b2 and e-value threshold $evalb
	my $intersect=interval::intersect($b1,$b2,@C0->[$i][2],@C0->[$i][3]);
	if ($intersect != 0 ) {		
	  my $cover=($intersect->[1]-$intersect->[0]+1)/($b2-$b1+1);
	  my $cover2=($intersect->[1]-$intersect->[0]+1)/(@C0->[$i][3]-@C0->[$i][2]+1);
	  if ($cover >= 0.9 && $cover2 >=0.9) {
	  
	    if (@C0->[$i][1] <= $evalb ) {
	      push @Cfilt, [$C0[$i][0], $C0[$i][1],$C0[$i][2],$C0[$i][3],$C0[$i][4],$C0[$i][5],$C0[$i][6]];
	    }
	  }
	}
      }
      @C0=();
      
##### resolution of redundant hits ############################################################

      unless (scalar(@Cfilt)==0) {
	my @grenzen; #array: 1 times two elements (begin,end,proteinpos1-2,highest e-value,strang) for the filtered C0 cluster Cfilt
	my $genloca = $Cfilt[0][4]; #open cluster (start)
	my $genlocb = $Cfilt[0][5]; #close cluster (end)
	my $protloca =  $Cfilt[0][2]; #prot (start)
	my $protlocb = $Cfilt[0][3]; #prot (end)
	my $eva = $Cfilt[0][1]; #best e-value
	my $strand = $Cfilt[0][6]; # + ... plus-strand; - ... minus-strand;

	for (my $j=0;$j<@Cfilt;$j++)
	{
	  my $overlap=interval::overlap($genloca,$genlocb,@Cfilt->[$j][4],@Cfilt->[$j][5]);
	  if ($overlap != 0) { #there is an overlap
	    $genloca = $overlap->[0]; # auch funktionieren sollte: @Cfilt->[$j][4];
	    $genlocb = $overlap->[1]; # auch funktionieren sollte: @Cfilt->[$j][5];
	    my $protover=interval::overlap($protloca,$protlocb,@Cfilt->[$j][2],@Cfilt->[$j][3]);
	    $protloca = $protover->[0];
	    $protlocb = $protover->[1];
	    
	    if ($eva > @Cfilt->[$j][1]) {
	      $eva = @Cfilt->[$j][1];
	    }
	    #look for errors
	    if  ($strand =~/^\+/ & $Cfilt[$j][6] =~/^\+/) {
	      #both plus orientation -> ok
	    }
	    elsif ($strand =~/^-/ & $Cfilt[$j][6]=~/^-/ ) {
	      #both minus orientation -> ok
	    }
	    else {
	      #should not happen (that overlapping hits are of different orientation)!
	      print STDERR "Warning! Bad cluster!\n";
	    }
	  }
	  else 
	  {
            #close cluster and open new cluster
	    my $ori;
	    if ($strand =~ /^\+/) {
	      $ori="+";
	    }
	    elsif ($strand =~/^-/) {
	      $ori="-";
	    }
	    push @grenzen, [$genloca,$genlocb,$protloca,$protlocb,$eva,$ori];
	    $genloca = @Cfilt->[$j][4];
	    $genlocb = @Cfilt->[$j][5];
	    $protloca = @Cfilt->[$j][2];
	    $protlocb = @Cfilt->[$j][3];
	    $eva =  @Cfilt->[$j][1];
	    $strand = $Cfilt[$j][6]; 
	  }
	}
	# handling of last cluster
	my $ori;
	if ($strand =~/^\+/) {$ori="+"}
	elsif ($strand =~/^-/) {$ori="-"}


####### check relative localization of segments to one another & print output ##################

	push @grenzen, [$genloca,$genlocb,$protloca,$protlocb,$eva,$ori];

	print "---\n";
	print "# Number of hits: " . scalar(@grenzen) . "\n";
	
	#separates hits on the plus strand from hits on the minus strand
	my @grenzenplus;
	my @grenzenminus;
	for (my $i=0; $i<scalar(@grenzen); $i++){
	  if (${$grenzen[$i]}[5] eq "+") {push @grenzenplus, @grenzen[$i];}
	  else {push @grenzenminus, @grenzen[$i];}
	}
	
	if (scalar(@grenzenplus)==1) { #if less than 2 hits, this is not the locus we want
	  print "# $vieh $Cfilt[0][0] @{$grenzenplus[0]}\n";
	}
	else {
	  my $first;
	  my $second;
	  for (my $i=0; $i<scalar(@grenzenplus); $i++) {
	    if (defined($first)) {
	      my $isex2=isexon2(\@{$grenzenplus[$i]});
	      if ($isex2==1) {
		
		$second=$i;
		# print first pair
		print "$vieh $Cfilt[0][0] @{$grenzenplus[$first]}\n";
		print "$vieh $Cfilt[0][0] @{$grenzenplus[$second]}\n";
		# set first and second back
		undef($first);
		undef($second);
	      }
	      else {
		print "# $vieh $Cfilt[0][0] @{$grenzenplus[$first]}\n";
		$first=$i;
	      }
	    }
	    else { #first is undef
	      if (isexon2(\@{$grenzenplus[$i]})==1) {
		print "# $vieh $Cfilt[0][0] @{$grenzenplus[$i]}\n";
		next;
	      }
	      else {
		$first=$i;
	      }
	    }
	  }
	}
	if (scalar(@grenzenminus)== 1) { #if less than 2 hits, this is not the locus we want
	  print "# $vieh $Cfilt[0][0] @{$grenzenminus[0]}\n";
	}
	else {
	  my $first;
	  my $second;
	  for (my $i=0; $i<scalar(@grenzenminus); $i++) {
	    if (defined($first)) {
	      my $isex2=isexon2(\@{$grenzenminus[$i]});
	      if ($isex2==1) {
		print "# $vieh $Cfilt[0][0] @{$grenzenminus[$first]}\n";
		$first=$i;
		
	      }
	      else {
		$second=$i;
		# print first pair
		print "$vieh $Cfilt[0][0] @{$grenzenminus[$second]}\n";
		print "$vieh $Cfilt[0][0] @{$grenzenminus[$first]}\n";
		# set first and second back
		undef($first);
		undef($second);
	      }
	    }
	    else { #first is undef
	      if (isexon2(\@{$grenzenminus[$i]})==1) {
		$first=$i;
	      }
	      else {
		print "# $vieh $Cfilt[0][0] @{$grenzenminus[$i]}\n";
		next;
	      }
	    }
	  }
	}
      }
      @Cfilt =();
    }
  }
}

###############################################################################################
  
sub isexon2 {
  my $exon=shift;
  
  my $intersect=interval::intersect($b1,$b2,$exon->[2],$exon->[3]);
  my $bool=0;
  
  if ($intersect != 0 ) {						
    my $cover=($intersect->[1]-$intersect->[0]+1)/($b2-$b1+1);
    if ($cover >= 0.9) {
      $bool=1;   
    }
  }
  # print "exons $exon->[2] $exon->[3] overlap $a1 $a2 --> $bool\n";
  return $bool;
}


#split into directory, basename and extention
#takes: filename with or without path
#returns: array holding three elements
#         [0] path
#         [1] basename
#         [2] extentions (starting with the dot)
sub parse_filename {
my $file=shift;
  my $path;
  my $filename;
  my $basename;
  my $extentions;
  if ($file =~ /\//) { # file has path and filename
  #  print STDERR "has path <$file>";
    if ($file =~ /^(.*\/)([^\/]*)$/) {
      $path=$1;
      $filename=$2;
   #   print STDERR "<$path>\n";
    }
  }
  else { # file has no path
    $filename=$file;
  }
  if ($filename =~ /\./) {# has an extension
    if($filename=~ /^([^\.]+)(\..*)/) {# get basename
      $basename=$1;
      $extentions=$2;
    }
  }
  else {
    $basename=$filename;
  }
  return $path,$basename,$extentions;
}

sub print_data_struct {
  my ($ref, $fh)=@_;
  foreach my $i (@$ref) {
    foreach my $j (@$i) {
      print $fh "$j ";
    }
    print $fh "\n";
  }
}
