#!/usr/local/bin/perl 
#######################################################################
# This software has been created by Genome Research Limited (GRL).    # 
# GRL hereby grants permission to use, copy, modify and distribute    # 
# this software and its documentation for non-commercial purposes     # 
# without fee at the user's own risk on the basis set out below.      #
# GRL neither undertakes nor accepts any duty whether contractual or  # 
# otherwise in connection with the software, its use or the use of    # 
# any derivative, and makes no representations or warranties, express #
# or implied, concerning the software, its suitability, fitness for   #
# a particular purpose or non-infringement.                           #
# In no event shall the authors of the software or GRL be responsible # 
# or liable for any loss or damage whatsoever arising in any way      # 
# directly or indirectly out of the use of this software or its       # 
# derivatives, even if advised of the possibility of such damage.     #
# Our software can be freely distributed under the conditions set out # 
# above, and must contain this copyright notice.                      #
#######################################################################

package ByTempl;
use FileHandle;

##########################################################################################
#
# Pre process the templates program output (rmott) 
#
# separates template data into array of arrays of good and bad stuff....
# these are returned to calling program



my @Contig_data;		# lines of CAF derived data 
my %CReads;
my @byCont_output;
my @byAssem_output;
my %Contig_names;               # contigs ..
my %ContComp;
my %Templates;			# names of templates used in assembly
my @Assmebly;			# ordered Contig names .. -> assembly
my @HasFwdRev;
my %Offsets;			# positions of starts of contigs from start of assembly
my $minDist;
my %ContLr;
my %Off;		     
my %ContLines;
my @Table;

##########################################################################################
## Get subs ...
##########################################################################################
sub Get_Offset{
   %Offsets;
}
sub Get_CLR{
    %ContLr;
}
##########################################################################################
sub Get_CR{
    %CReads;
}
##########################################################################################
sub Get_Off{
    %Off;
}
##########################################################################################
sub Get_ContigOrd{
  @Assembly;
}
##########################################################################################
sub Get_ContComp{
    %ContComp;
}
##########################################################################################
sub Get_Assem{
    @byAssem_output;
}
##########################################################################################
sub GetMaxBase{
 
}
##########################################################################################
sub Get_Assem_Max{
    my $max = $Offset{"end"};
    $max;
}
##########################################################################################
# returns template list

sub Get_Templ_Array{
    keys(%Templates);
}
##########################################################################################
# read the templates file into table structure
sub Filter{
    my ($f) = @_;
    my @newData;

    return if($f ==0 );
    if($f == 1)
    {
	my @conts;
	@newData = grep(((abs($_->[18] - $_->[17]) >= 1000) &&((abs($_->[20] - $_->[19]) >= 1000))),
			@byAssem_output);
	undef(@byAssem_output);
	@byAssem_output = @newData;
	return;
    }
}
sub Read_Templ{

    my ($filename , $dist ,$tempfile) = @_;      # filename .. min separation
    my $line;
    my $Cont_Keys;
    my $Cont;
    my %A;
    my $j = 0;
    my $i=0;
    
    $minDist = $dist;		      # set global dist		
   
    if($tempfile)
    {
	open(CONTLIST, $filename);
    }else{
	open(CONTLIST,"templates < $filename|");
    }

    $|=1;
    &Read(CONTLIST);
  
  


}
#######################################################################################################
sub Read{
    my ($CONTLIST) = @_;
   
    print "Crow : Reading Data \n";
    while($line = <$CONTLIST>)	
    {  	   
            
	my @matched = split('\s+',$line);
	if(!defined($Contig_names{$matched[0]})){$Contig_names{$matched[0]}= $j++;}   # Contig names
	push(@{$ContLines{$matched[0]}},$i); 
        push(@{$Templates{$matched[2]}},$i);  # builds list of all templates used in assembly... rows in cont data 
	$Contig_data[($i++)] = \@matched;    # builds array of lines of raw CAF data
	push (@{$CReads{$matched->[0]}},$matched->[2]);
    }
    close($CONTLIST);
    if(($#Contig_data < 0)){print "Crow : Hmmm no data ??!! Bye !!  \n";exit(-1);}

    foreach (keys(%Contig_names)) # all names
    {
        $ContComp{$_} = 1;
	$Assembly[$Contig_names{$_}] = $_; # build a list of the contigs and their order from file !!!
    }

    print "Crow : Calculating \n";				
    &Calc_offsets(@Assembly);	# calculate indentation of each contig from left most in assem
    &Build_CLR(@Assembly);
    &Build_CReads;
    print "Crow : PreProcessing Reads \n";
    &summary_byAssem;		# build summary table ..
    if($#byAssem_output < 0){print "Crow : Hmmm no data ??!! Bye !!  \n";exit(-1);}
    print "Crow : Checking Read Consistency\n"; 
    &Assem_Consist;             # return summary table  
}
########################################################################################################
sub RevComp{
    my ($c, $cont) = @_;

    if($ContComp{$cont} ==1 ){$ContComp{$cont} =0;}else {$ContComp{$cont} =1;}
    print "Crow : Reverse Compliment $cont \n";

    my @Affreads = grep(($_->[0] eq $cont),@Contig_data);

    my $cr = $ContLr{$cont}->[1] - $ContLr{$cont}->[0];
    my $temp;
    foreach $read (@Affreads)
    { 
     $read->[4] = $cr - $read->[4];	
     $read->[5] = $cr - $read->[5];
     if($read->[8] eq "negative"){$read->[8] = "positive";}else{$read->[8] = "negative";} 
    }		
    &summary_byAssem;
    print "Crow : Checking Read Consistency\n"; 
    &Assem_Consist; # return summary table
    &main::Draw_Reads;
  
}				
#########################################################################################################
# build list of reads for given contig
sub Build_CLR{

    my (@Ass) = @_;
   
    my $rt =0;
    my $max;
    my $cname;
    my $summax=0;
    
    
    foreach $cname (@Ass)
    {
	$max = &Contig_max($cname);
	$summax+=$max;
        push(@{$ContLr{$cname}},$rt, $summax);  # contig max
     	$rt = $summax;       
    }
   
    %ContLr;			# return 

}
##########################################################################################################
sub Build_CReads{

    my $cname;
    my $line;
    my $on =0;

    foreach $cname (@Contig_data)   # each contig
    {
      push (@{$CReads{$cname->[0]}},$cname->[2]); # ditto       
    }

}
##########################################################################################
sub ReadTable{
    my ($contig) = @_;
   
    my $sorted;
    my $found =0;
    undef(@Table);
    @Table = grep(($_->[0] eq $contig),@Contig_data);
  
    @Table;
}
##########################################################################################
# Change Contig Order ....
#
sub NewContOrder{

    my (@Ass) = @_;

    undef(@Assembly);
    @Assembly = @Ass;   
    &Calc_offsets(@Assembly);
    &summary_byAssem;		# build summary table ..
    print "Crow : Checking Read Consistency\n"; 
    &Assem_Consist;             # return summary table
}

#############################################################################################################
sub Calc_offsets{
   
    my (@Ass) = @_;
    my $l =1;
    my $rt;
    my $cname;

    $Offset{"start"} = $l;   # left most is 1
    
    foreach $cname (@Ass)
    {
        $rt = $l + &Contig_max($cname) +1;  # contig max
     	$Offset{$cname} = $l;
      	$l = $rt;       
    }
   
    $Offset{"end"} = $l;

    %Offset;			# return 
}
#############################################################################################################

##########################################################################################
# build summary table

sub summary_byAssem{

    my $tname;
    
    my $forwards_right=0;	#forward strand reads ... -> in assembly
    my $forwards_left=0;	#forward strand reads ... <- in assembly 
    my $reverse_right=0;	#reverse strand reads ... -> in assembly
    my $reverse_left=0;		#reverse strand reads ... <- in assembly 
    my $nName;			# 
    my $ls1=0;			# 
    my $ls2=0;			# 
    my $rs1=0;			# 
    my $rs2=0;			# 
    my $l;			# 
    my $rr;			# 
    my $cname;			# 
    my $nf;			# 
    my $nr;			# 
    my $offs;			# 
    my $i =0;
    my $limit = &main::GetLim;

   
    # if num read / num templates is large this is quick else slow
    # or lots of piddly contigs
    foreach(@byAssem_output){undef($_);}
    undef(@byAssem_output);
    
    foreach $tname (keys(%Templates))   # each template
    {
	my @reads;
        my @fwd;		# fwd for this temp
        my @rev;		# rev ditto
        my $lr;
        my @notdone = @{$Templates{$tname}};
       
        
	foreach $row (@notdone)
	{
	   
	    $nName = $Contig_data[$row];
	 
	    if($$nName[2] eq $tname)  # lines for this contig and templ. 
	    { 
		# count occurence of read types
		 #  print "$tname $$nName[2] $row\n";		
      		if($$nName[3] eq "forward")
		{
		    if($$nName[4] < $$nName[5]){$forwards_right++;}else{$forwards_left++;}
		    push(@fwd,[@{$nName}]);   		   
		} 
		if($$nName[3] eq "reverse")
		{
		    if($$nName[4] < $$nName[5]){$reverse_right++;}else{$reverse_left++;}		
		    push (@rev,[@{$nName}]);
		}
	       
	    }
		
	}

	undef(@notdone);
        $nf = $#fwd;		# number of fwd
        $nr = $#rev;		# ditto rev
        my $pt =0;
        my $hfr =0;		# BOOL has fwd and rev
        my $inlim = 0;		# BOOL is <= sep limit
        my $fr =0;
        my $row;
        my @l;			# <-
        my @r;			# ->
        my @lread;
        my @rread;
    
        foreach $row (@fwd)	# adjust the co-ordinates from contig to assem
	{
         $$row[4] += ($Offset{$$row[0]});
         $$row[5] += ($Offset{$$row[0]});             
	}
        foreach $row (@rev)	# ditto rev
	{
         $$row[4] += ($Offset{$$row[0]});
         $$row[5] += ($Offset{$$row[0]});
	}
     
        if(($nf > -1) && ($nr > -1)) # have fwd and rev
          {
	      $hfr =1;		# true
	      @lread = &leftread(\@fwd,\@rev);  # get leftmost of set
	      @rread = &rightread(\@fwd,\@rev);	# get right most
           
	      $ls1 = $lread[4];	# get co-ord of leftmost
	      $ls2 = $lread[5]; # 
	      $rs1 = $rread[4]; # ditto right
	      $rs2 = $rread[5]; #
	      
	      $pt = &Point_to($ls1,$ls2,$rs1,$rs2);               # BOOL point towards -> <-
	      $inlim = &OK_apart($ls1,$ls2,$rs1,$rs2,$limit);       # check separation
	      
	   }
        if(($nf > -1) && ($nr == -1))        ## fwd only
          {
	      # as above but fwd only
	      $hfr =0;
	     
	      @lread = &leftread(\@fwd,\@fwd);
	      @rread = &rightread(\@fwd,\@fwd);
	      $ls1 = $lread[4];
	      $ls2 = $lread[5];
	      $rs1 = $rread[4];
	      $rs2 = $rread[5];
	      $pt =0;
	      $inlim = &OK_apart($ls1,$ls2,$rs1,$rs2,$limit);
	      # print "==> $#fwd $lread[0] $lread[2] $ls1 $ls2  $rs1 $rs2 pt == $pt inlim == $inlim\n";
          }
        if(($nf == -1) && ($nr > -1))        ## rev only
         {
           # as above rev only
	     $hfr =0;
	     @lread = &leftread(\@rev,\@rev);
	     @rread = &rightread(\@rev,\@rev);
	     $ls1 = $lread[4];
	     $ls2 = $lread[5];
	     $rs1 = $rread[4];
	     $rs2 = $rread[5];
	     $pt =0;
	     $inlim = &OK_apart($ls1,$ls2,$rs1,$rs2,$limit);
	     # print "==> $#rev $lread[0] $lread[2] $ls1 $ls2  $rs1 $rs2 pt == $pt inlim == $inlim\n";
         }
        # now prob of biggest ....
	
        # temp name sum_fr sum_fl sum_rr sum_rl hasfwdrev pointtowards isinlimit 0 co-ord_ls1 ls2 rs1 rs2

        push(@byAssem_output,[$tname,$forwards_right,$forwards_left,$reverse_right,
			      $reverse_left,$hfr,$pt,
			      $inlim,0,$ls1,$ls2,$rs1,$rs2,
			      $lread[0],$rread[0],$lread[3],$rread[3],
			      $ContLr{$lread[0]}->[0],$ContLr{$lread[0]}->[1],
			      $ContLr{$rread[0]}->[0],$ContLr{$rread[0]}->[1]]);
      
	# reset some values 
     
	$forwards_right=0;
	$forwards_left=0;
	$reverse_right=0;
	$reverse_left=0;
    }
   

}
########################################################################################
sub QuickReport{

    my ($fname , $limt) = @_;
    my $rname = $ENV{HOME}."/crow.report";
    my @thisCont;
    my ($cnt,$cstart,$cstop,$size,$nr,$nt,$ng,$nb,$nl,$nf,$nrv);
    my $tot =0;     
    my $nrtot=0;
    my $nftot=0;
    my $nrvt =0; 
    my $nttot =0;

    print "Crow : Writing Report $rname\n";

 
    # crow reporting 
format STDOUT_TOP =
Crow        :Quick Summary
FILE        : @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<   
              $fname
SEP LIMIT   : @<<<<<<<<<<<< 
              $limit                

Contigs are in current assembly order
Good ,Bad and Lone Groups refer to visible groups which may or may not contain reads in more than one Contig
 
Contig Name  Contig_Start Contig_Stop  Size    Num Reads (f)   (r)    Num Templ  Num Good Grps  Num Bad Grps  Num Lone Grps
----------------------------------------------------------------------------------------------------------------------------
#
.


write ;
    
format STDOUT =
@|||||||||||@||||||||||||@||||||||||||@||||||||@|||||||||@||||||@|||||@||||||||@||||||||||||||@|||||||||||||||@|||||||||||
   $cnt,       $cstart,     $cstop,     $size,   $nr,      $nf,   $nrv,  $nt,        $ng,           $nb,            $nl
.
 

    my @t;
   

    foreach $cont (@Assembly)
    {
	my %count;
	$cnt = $cont;
        $cstart = $ContLr{$cnt}->[0];
	$cstop = $ContLr{$cnt}->[1];
	$size = $cstop-$cstart;	
	$tot+= $size;
	$nr = grep(($_->[0] eq $cnt),@Contig_data);
	$nrtot+= $nr;
	foreach ( @{$CReads{$cnt}})
	{
	    $count{$_}++;
	}
        @t = keys(%count);
	$nt = $#t +1;
	$nttot += $nt;
        $nf = grep(($_->[0] eq $cnt)&& ($_->[3] eq "forward"),@Contig_data);
	$nftot += $nf;
	$nrv = grep(($_->[0] eq $cnt)&& ($_->[3] eq "reverse"),@Contig_data);
	$nrvt+= $nrv;
	$ng = grep(((($_->[13] eq $cnt) || ($_->[14] eq $cnt)) 
		    && ($_->[5]) &&($_->[6]) && ($_->[7]) && (!$-->[8])),@byAssem_output);
	$nb = grep(((($_->[13] eq $cnt) || ($_->[14] eq $cnt)) 
		    && ($_->[5]) && ((!$_->[6]) || (!$_->[7]))),@byAssem_output);
	$nl = grep(((($_->[13] eq $cnt) || ($_->[14] eq $cnt)) 
		    && (!$_->[5])),@byAssem_output);
	write;
    }
my $y ="";
format STDT =

---------------------------------------------------------------------------------------------------------------------------
Total
@|||||||||||@||||||||||||@||||||||||||@||||||||@|||||||||@||||||@|||||@||||||||
    $y,         $y,            $y,        $tot,  $nrtot,  $nftot,$nrvt, $nttot
---------------------------------------------------------------------------------------------------------------------------
.     
 
#    close(REPORT);
    STDOUT->format_name("STDT");
   write;
   print "Crow : Report written\n";
}
########################################################################################
# change min allowed separation of good fr reads

sub ChangeLimit{
    my ($lim) = @_;
    my $line;
    my $inlim;

    foreach $line (@byAssem_output)
    {
     $inlim = &OK_apart($$line[9],$$line[10],$$line[11],$$line[12],$lim);  # test separation
     $$line[7] = $inlim;	# adjust summary table
    }
   
}
########################################################################################
sub leftread{

# if left most fwd read is left of left most rev read ... left read = leftmost fwd read else leftmost rev

  my ($fwd,$rev) = @_;
  my $fr =0;
  my @f = @{$fwd};
  my @r = @{$rev};
  my @lf = left_most(@f);
  my @lr = left_most(@r);
  my @leftread;
 
  if(&fr($lf[4],$lf[5],$lr[4],$lr[5])){@leftread = @lf;}else{@leftread = @lr;}

  @leftread;
}
#######################################################################################
sub rightread{

# if right most f is left of right most rev rightmost = rightmost rev

  my ($fwd,$rev) = @_;
  my $fr =0;
  my @f = @{$fwd};
  my @r = @{$rev};
  my @rf = right_most(@f);
  my @rr = right_most(@r);
  my @rightread;

  if(&fr($rf[4],$rf[5],$rr[4],$rr[5])){@rightread = @rr;}else{@rightread = @rf;}
  
  @rightread;
}
######################################################################################
sub left_most{
  my (@fr) =@_;
  my @ls1;
  my $lefts1;
  my @ls2;
  my $lefts2;
#  my @left_most;
  my @least;

  @fr = sort by_s1 @fr;
  @ls1 = @{@fr[0]}; 

  @fr = reverse sort by_s2 @fr;
  @ls2 = @{@fr[0]};
 
  if($ls1[4] < $ls2[5]){@least = @ls1;}else {@least = @ls2;}
  @least;       # leftmost read for template in assembly !!
}
####################################################################################################
sub right_most{
  my (@fr) =@_;
  my @ls1;
  my @ls2;
  my @left_most;
  my @most;
  my @fwdrev;

  @fr = reverse sort by_s1 @fr;
  @ls1 = @{@fr[0]};

  @fr =  sort by_s2 @fr;
  @ls2 = @{@fr[0]};
 
  if($ls1[4] > $ls2[5]){@most = @ls1;}else {@most = @ls2;}
 
  @most;       # rightmost read for template in assembly !!
}
#####################################################################################################
# do two given reads point towards

sub Point_to{

my ($ls1,$ls2,$rs1,$rs2) = @_;
my $pt =0;                             #point same way

   if(($ls2 > $ls1) && ($rs2 < $rs1))  #leftmost points right and rightmost points left
    {
     $pt =1;           
    }
  if(($ls2 < $ls1) && ($rs2 > $rs1))   #leftmost points left and rightmost points right
    {
     $pt = -1;           
    }
 $pt;
}
######################################################################################################
# is their separation <= limit

sub OK_apart{

my ($ls1,$ls2,$rs1,$rs2,$lim) =@_;
my $left;
my $right;
my $OK =0;

if($ls1 < $ls2){$left = $ls1;}else{$left = $ls2;}
if($rs1 < $rs2){$right = $rs2;}else{$right = $rs1;}

if((abs($right - $left)) <= $lim){$OK = 1;}

$OK;

}
#####################################################################################################
# are they internally consistent

sub Assem_Consist{
my $outline;

foreach $outline (@byAssem_output)   # use summary table
{ 
    
   if((($$outline[1] > 0) && ($$outline[2] > 0)) || (($$outline[3] > 0) && ($$outline[4] > 0)))
   {
       $$outline[8] =1;
   #    print "internally inconsistent reads -><-... $$outline[0]\n";
   }

   #if a template has forward and reverse reads in the same contig then they should be -><-

   if((($$outline[1] > 0) || ($$outline[2] > 0)) && (($$outline[3] > 0) || ($$outline[4] > 0)))
    {
	if((($$outline[1] >0) && ($$outline[3] >0)) || (($$outline[2] > 0) && ($$outline[4] > 0)))
	{
            $$outline[8] = 1;
	#    print "template with forward and rev reads with differing abs dirs .. $$outline[0] \n";
	}	
    } 
}	  

}
#########################################################################################
# sorting

sub by_s1{
    $$a[4]<=>$$b[4];
}
sub by_s2{
    $$b[5]<=>$$a[5];
}
##########################################################################################
sub Template_xtents_by_Cont{

    my ($cname,$tname) = @_;
    my $line;
    my @fwdrev;
    # template

    foreach $line  (@Contig_data)
    {
     
	if(($$line[0] eq $cname) && ($$line[2] eq $tname))
	{
	    push(@fwdrev,$$line[4]);
	    push(@fwdrev,$$line[5]); 
	}
    }
    @fwdrev = sort by_val @fwdrev;
   
   # if($cname eq "Contig7"){foreach (@fwdrev){print"$_\t";}}
   # print "\n";
 @fwdrev;
}
##########################################################################################
sub by_val{
    $a<=>$b;
}
#############################################################################################
sub Contig_min{
    my ($cont)= @_;
    my $nLine;
    my @vals;
    my $n =0;

    foreach $nLine (@Contig_data)
    {
	if($$nLine[0] eq $cont)
	{
	    $vals[$n++] = $$nLine[4];
            $vals[$n++] = $$nLine[5];
	}
    }				
    
    @vals = sort by_val @vals;

    @vals[0];
}
###########################################################################################
sub Contig_max{
    my ($cont) =@_;
    my $nLine;
    my @lines = $ContLines{$cont};
    my @vals;
    my $n =0;
    my @line;
    my $i;

    foreach $nLine (@{$ContLines{$cont}})
    {
     push(@vals,$Contig_data[$nLine]->[4],$Contig_data[$nLine]->[5]);        
    }
	
    @vals = reverse sort by_val @vals;
 #   if($cont eq "Contig14.0"){print "$cont   $vals[0]  \n";}
    $vals[0];
}				# 
##########################################################################################

#################################################################################################################
sub fr{
 my ($fs1,$fs2,$rs1,$rs2) = @_;
 my $fr =0;

 if(($fs1 <= $rs2) && ($fs2 <= $rs1)){$fr =1;} # fwd is to left of rev
 if(($rs1 <= $fs2) && ($rs2 <= $fs1)){$fr =0;} # rev left of fwd ... this line not necessary ..

 $fr;
}	
################################################################################################################
sub Point_towards{		# priority to write this !!!!
    my ($fs1,$fs2,$rs1,$rs2) = @_;

    my $pt =0;
    #print "$fs1\t$fs2\t$rs1\t$rs2\n";
    
    if(($fs2 > $fs1) && ($rs2 < $rs1))  #fwd points right and rev points left
    {
      if(&fr($fs1,$fs2,$rs1,$rs2)) #fwd is to left of rev
      {
       
	  $pt = 1;
	
      }
    }

    if(($rs2 > $rs1) && ($fs2 < $fs1))  #rev points right and fwd points left
    {
	if(!&fr($fs1,$fs2,$rs1,$rs2)) #rev is to left of fwd
	{	    	    
	    $pt = 1;
	}
    }				
       $pt;
}
##########################################################################################################
sub by_temp_name{
    $$a[2] cmp $$b[2];
}
sub by_cont_name{
    $$a[0] cmp $$b[0];
}
##########################################################################################################
1;


