#!/usr/bin/perl -w



use strict;
use POSIX qw(ceil floor);


die "Usge: $0  eland_file  tag_length  base_output_name\n" unless (@ARGV == 3);
my ($elandFile, $tagLen, $baseName, $startType, $isSub) = @ARGV;


# parse Eland mapping file.
# %tagMap -- each key is chr, each corresponding value is an anonymous array of tags
my %tagMap;
parseElandFile($elandFile, \%tagMap, $tagLen, $baseName);


#
my (%tagStackLocs, %tagSignalMap);
processTagMapping(\%tagMap, \%tagStackLocs, \%tagSignalMap);




# 
my ($tagClustersTxtFileName, $tagSignalWigFileName);
$tagClustersTxtFileName = "${baseName}_tag-clusters_orig.txt";
$tagSignalWigFileName = "${baseName}_tag-signal_orig.wig";



# 
open(TXT, ">$tagClustersTxtFileName");
print TXT "#chr\toutStart\toutEnd\toutLen\toutCount\tinStart\tinEnd\tinLen\tinCount\n";

foreach my $chr (sort {$a <=> $b} keys %tagStackLocs)
{
  my @tagRefs = @{ $tagStackLocs{$chr} };
  foreach my $tagRef (sort {$a->[1] <=> $b->[1]} @tagRefs)
  {
    print TXT join("\t", @{$tagRef}), "\n";
  }
}


# pirnt the tag signal map
if (0)
{
  open(WIG, ">$tagSignalWigFileName");
  
  foreach my $chr (sort {$a <=> $b} keys %tagSignalMap)
  {
    my @tagRefs = @{ $tagSignalMap{$chr} };
    
    if    ($chr == 23) { $chr = "X"; }
    elsif ($chr == 24) { $chr = "Y"; }
    $chr = "chr".$chr;
    
    print WIG "track type=wiggle_0\n";
    foreach my $tagRef (sort {$a->[1] <=> $b->[1]} @tagRefs)
    {
      print WIG join("\t", ($chr, @{$tagRef}[1..3])), "\n";
    }
  }
}






sub processSignalMap
{
  my ($tagSignalMapRef) = @_;
  my %tagSignalMap = %{ $tagSignalMapRef };

  my %tagClusters;
  
  foreach my $chr (keys %tagSignalMap)
  {
    my @stepRefs = sort {$a->[1] <=> $b->[1]} @{ $tagSignalMap{$chr} };
    push @stepRefs, [$chr, 1000000000, 1000000005, 0];

    my @clusterStack;
    for (my $i = 0; $i < @stepRefs; $i++  )
    {
      if (@clusterStack == 0)
      {
        push @clusterStack, [ @{$stepRefs[$i]} ];
      }
      else
      {
        if ($clusterStack[-1][2] >= $stepRefs[$i][1])
        {
          push @clusterStack, [ @{$stepRefs[$i]} ];
        }
        else # overlap breaks, pop the stack
        {	  
	  my @bases;
	  foreach my $stepRef (@clusterStack)
	  {
	    my ($sChr, $sStart, $sEnd, $sHeight) = @$stepRef;
	    push @bases, ( ($sHeight) x ($sEnd-$sStart+1) );
	  }
	  my ($bestOverlapIndStart, $bestOverlapIndEnd, $highestNumOfOverlappingTags) = getBestOverlap(@bases);

	  my ($outStart, $outEnd) = ($clusterStack[0][1], $clusterStack[-1][2]);
	  my $outLen = $outEnd - $outStart + 1;
	  my $outCount = ceil(sum(@bases)/$tagLen);
	  my ($inStart, $inEnd) = ($bestOverlapIndStart+$outStart, $bestOverlapIndEnd+$outStart);
	  my $inLen = $inEnd - $inStart + 1;
	  my $inCount = $highestNumOfOverlappingTags;
	  push @{ $tagClusters{$chr} }, [$chr, $outStart, $outEnd, $outLen, $outCount, $inStart, $inEnd, $inLen, $inCount];
	  
          @clusterStack = ();
          push @clusterStack, [ @{$stepRefs[$i]} ];
        }
      }
    }
  }

  return(\%tagClusters);
}





# %tagMap -- input data: chr, start, end of each mapped tag
# %tagStackLocs -- output data: chr, start, end, outer, ... of each tag cluster
# %tagSignalMap -- output data: chr, start, end, stack height of each step in eacg tag cluster

sub processTagMapping
{
  my ($tagMapRef, $tagStackLocsRef, $tagSignalMapRef) = @_;
  my  %tagMap = %{ $tagMapRef };
  
  my $count = 0;
  foreach my $chr (sort keys %tagMap)
  {
    my @overlapTagStack; # holding overlapping tags
    
    # sort (must!) on the starts of the tags
    my @tags = sort {$a->[0] <=> $b->[0]} @{ $tagMap{$chr} };
    push @tags, [ 1000000000, 1000000005 ]; # artificial tag for pop the last real tag stack
    
    for (my $i = 0; $i < @tags; $i++  )
    {
      if (@overlapTagStack == 0)
      {
        push @overlapTagStack, [ @{$tags[$i]} ];
      }
      else
      {
        if ($overlapTagStack[-1][1] >= $tags[$i][0])
        {
          push @overlapTagStack, [ @{$tags[$i]} ];
        }
        else # overlap breaks, pop the stack
        {
          $count += @overlapTagStack;
          processTagStack($chr, \@overlapTagStack, $tagStackLocsRef, $tagSignalMapRef);
          @overlapTagStack = ();
          push @overlapTagStack, [ @{$tags[$i]} ];
        }
      }
    }
  }
  
  # sort tag grous on each chr
  foreach my $chr (sort keys %{$tagStackLocsRef})
  {
    my @tagStackRefs = @{ $tagStackLocsRef->{$chr} };
    @tagStackRefs = sort {$a->[1] <=> $b->[1]} @tagStackRefs;
    $tagStackLocsRef->{$chr} = [ @tagStackRefs ];
  }
}


sub processTagStack
{
  my ($chr, $overlapTagStackRef, $tagStackLocsRef, $tagSignalMapRef) = @_;
  my @overlapTagStack = @{ $overlapTagStackRef };

  my $stackStart = $overlapTagStack[0][0]; # the start of the first
  my @sorted = sort {$a->[1] <=> $b->[1]} @overlapTagStack; # sort on the end
  my $stackEnd = $sorted[-1][1]; # the end of the sorted last
  
  my $endInd = $stackEnd - $stackStart;
  my @temp;
  @temp[0..$endInd] = map{0} (0..$endInd); # initialize all elements to 0
  foreach my $ref (@overlapTagStack) # stack
  {
    my ($tagStart, $tagEnd) = @$ref;
    my $indStart = $tagStart - $stackStart;
    my $indEnd   = $tagEnd   - $stackStart;
    @temp[$indStart..$indEnd] = map {$_+1} @temp[$indStart..$indEnd];
  }
  
  my ($bestOverlapIndStart, $bestOverlapIndEnd, $highestNumOfOverlappingTags) = getBestOverlap(@temp);
  my $bestOverlapGenomicStart = $bestOverlapIndStart + $stackStart;
  my $bestOverlapGenomicEnd   = $bestOverlapIndEnd   + $stackStart;
  
  push @{ $tagStackLocsRef->{$chr} }, [ $chr, $stackStart, $stackEnd, $stackEnd-$stackStart+1, scalar(@overlapTagStack),
      $bestOverlapGenomicStart, $bestOverlapGenomicEnd, $bestOverlapGenomicEnd-$bestOverlapGenomicStart+1, $highestNumOfOverlappingTags ];

  my @stepStacks = getSignalSteps(@temp);
  foreach my $stepStackRef (@stepStacks)
  {
    my ($stepStackStartInd, $stepStackEndInd, $stepStackSig) = @$stepStackRef;
    push @{ $tagSignalMapRef->{$chr} }, [ $chr, $stepStackStartInd+$stackStart, $stepStackEndInd+$stackStart+1, $stepStackSig ];
  }

}


sub getSignalSteps
{
  my @stepStacks;
  my ($stepStackStartInd, $stepStackEndInd, $stepStackSig) = (0, 0, 0);
  
  my @bases = (@_, 0); # add 0, which doesn't occur in @_, to force a possible last highest overlap
  for (my $i = 0; $i < @bases-1; $i++)
  {
    if ($bases[$i] != $bases[$i+1])
    {
      $stepStackEndInd = $i;
      $stepStackSig = $bases[$i];
      push @stepStacks, [$stepStackStartInd, $stepStackEndInd, $stepStackSig];
      ($stepStackStartInd, $stepStackEndInd, $stepStackSig) = ($i+1, 0, 0);
    }
  }
  
  return(@stepStacks);
}


sub getBestOverlap
{
  my @bases = (@_, 0); # add 0, which doesn't occur in @_, to force a possible last highest overlap
  my $highestNumOfOverlappingTags = (sort {$b <=> $a} @_)[0];
  
  my %highestOverlaps;
  my $start = -1;
  for (my $i = 0; $i < @bases; $i++)
  {
    if ($bases[$i] == $highestNumOfOverlappingTags)
    {
      if ($start == -1)
      {
        $start = $i;
      }
    }
    else
    {
      if ($start != -1)
      {
        $highestOverlaps{$start} = $i-1; # end
        $start = -1;
      }
    }
  }
  
  my @sortedStarts = sort {$highestOverlaps{$b}-$b <=> $highestOverlaps{$a}-$a} keys %highestOverlaps; # sort on length
  return($sortedStarts[0], $highestOverlaps{$sortedStarts[0]}, $highestNumOfOverlappingTags);
}


sub addingMinDistancesFromMrnas
{
  print STDERR "Adding minimum distances from mrnas ...\n";
  my ($geneFile, $tagStackLocsRef) = @_;
  
  my %txstarts;
  open(GEN, $geneFile);
  foreach my $line (<GEN>)
  {
    $line =~ s/^\s+//;
    $line =~ s/\s+$//;
    my @fs = split(/\t/, $line);
    
    my $ori = $fs[9];
    my $chr = $fs[14];
    $chr =~ s/chr//i;
    $chr = 23 if ($chr eq "x" or $chr eq "X");
    $chr = 24 if ($chr eq "y" or $chr eq "Y");

    if ($ori eq "+") { push @{ $txstarts{$chr} }, [ 1, $fs[16] ]; }
    else             { push @{ $txstarts{$chr} }, [ 2, $fs[17] ]; }
  }

  my @minDistancesFromGenes;
  foreach my $chr (keys %{$tagStackLocsRef})
  {
    my @txstarts = @{ $txstarts{$chr} };
    my @tagStackRefs = @{ $tagStackLocsRef->{$chr} };
    my @tagStackMids;
    for (my $i = 0; $i < @tagStackRefs; $i++)
      { $tagStackMids[$i] = [ (-1)*$i, int(($tagStackRefs[$i][1]+$tagStackRefs[$i][2])/2) ]; }
    my @txsTagStackMix = sort {$a->[1] <=> $b->[1]} (@txstarts,@tagStackMids); # sort on location

    for (my $i = 0; $i < @txsTagStackMix; $i++)
    {
      if ($txsTagStackMix[$i][0] <= 0) # this is a tag stack
      {
        my $upMrnaStart = 0;
        for (my $j = $i; $j >= 0; $j--) # search upstream for '-' mrna
        {
          if ($txsTagStackMix[$j][0] == 2)
          {
            $upMrnaStart = $txsTagStackMix[$j][1];
            last;
          }
        }
        
        my $dnMrnaStart = 0;
        for (my $j = $i; $j < @txsTagStackMix; $j++) # search upstream for '+' mrna
        {
          if ($txsTagStackMix[$j][0] == 1)
          {
            $dnMrnaStart = $txsTagStackMix[$j][1];
            last;
          }
        }
        
        my $midPoint = $txsTagStackMix[$i][1];
        if ($upMrnaStart != 0 and $dnMrnaStart != 0)
        {
          push @{ $tagStackRefs[ (-1)*$txsTagStackMix[$i][0] ] }, min(($midPoint-$upMrnaStart), ($dnMrnaStart-$midPoint));
        }
        elsif ($upMrnaStart == 0 and $dnMrnaStart != 0)
        {
          push @{ $tagStackRefs[ (-1)*$txsTagStackMix[$i][0] ] }, $dnMrnaStart-$midPoint;
        }
        elsif ($upMrnaStart != 0 and $dnMrnaStart == 0)
        {
          push @{ $tagStackRefs[ (-1)*$txsTagStackMix[$i][0] ] }, $midPoint-$upMrnaStart;
        }
        else
        {
          push @{ $tagStackRefs[ (-1)*$txsTagStackMix[$i][0] ] }, -1111;
        }
      }
    }

    $tagStackLocsRef->{$chr} = [ @tagStackRefs ];
  }
}


sub sum
{
    my @nums = @_;
    my $sum = 0;

    foreach my $num (@nums)
    { $sum += $num; }

    return $sum;
}


sub min
{
    my @nums = @_;
    @nums = sort {$a <=> $b} @nums;

    return $nums[0];
}


sub round
{
    my($number) = shift;
    return int($number + .5 * ($number <=> 0));
}



sub parseElandFile
{
  my ($elandFile, $hashRef, $tagLen, $baseName) = @_;

  #open(RDS, ">${baseName}_reads.bed");
  #open(TGS, ">${baseName}_tags.bed");

  open(ELD, $elandFile) or die "cannot read $elandFile: $!\n";

  foreach my $line (<ELD>)
  {
    chomp $line;
    my @fs = split(/\t/, $line);
    if ($fs[2] =~ /U\d/) # unique match
    {
      my $chr = $fs[6];
      $chr =~ s/\.fa//;
      
      my $pos = $fs[7];
      my $ori = $fs[8];
      
      my ($rstart, $rend, $tstart, $tend);
      
      # reads directly from eland file
      if ($ori eq "F") {$rstart = $pos; $rend   = $rstart + 27 - 1;}
      else             {$rend   = $pos; $rstart = $rend   - 27 + 1;}
      
      # tags -- extended reads
      if ($ori eq "F") {$tstart = $pos; $tend   = $tstart + $tagLen - 1;}
      else             {$tend   = $pos; $tstart = $tend   - $tagLen + 1;}
      
      $chr =~ s/chr//i;
      $chr = 23 if ($chr eq "x" or $chr eq "X");
      $chr = 24 if ($chr eq "y" or $chr eq "Y");
      
      #print RDS "chr$chr\t$rstart\t$rend\n";
      #print TGS "chr$chr\t$tstart\t$tend\n";
      
      push @{ $hashRef->{$chr} }, [ $tstart, $tend ]; # chr, ( [start, end], ... )
    }
  }
}