#!/usr/bin/perl -w
#Author: Feng Cheng
use warnings;
use strict;

my $syn = $ARGV[0];
my $chrM = $ARGV[1];
my $out1 = $syn.".seg";
my $out2 = $out1.".merge";
my $out3 = $out2.".pos";

#system("awk '{print \"a:\"\$1\"\\ta:\"\$2\"\\t\"\$3\"\\t\"\$4\"\\tb:\"\$5\"\\tb:\"\$6\"\\t\"\$7\"\\t\"\$8}' $syn > $syn.re");
&processSyn($syn,$syn.".re");

&process($syn.".re",$chrM,$out1,$out2,$out3);

sub process {
  my ($syn,$chrM,$out1,$out2,$out3) = @_;

  my (%gene2pos,%a2b,%agenes2i,%bgenes2i,@agenes,@bgenes);
  &readSyn($syn,\%a2b,\%agenes2i,\%bgenes2i,\@agenes,\@bgenes,\%gene2pos,$chrM);

  my $itnum = 150;    # b abnormal location
  my @data;
  &pickSeg($out1,\%a2b,\%agenes2i,\%bgenes2i,\@agenes,\@bgenes,\%gene2pos,$itnum,\@data);

  my $spnum = 150;
  my $mingF = 10;
  &merge($out2,$out3,$spnum,\@data,\@agenes,\@bgenes,\%gene2pos,\%a2b,$mingF);

}

sub merge {
  my ($out2,$out3,$spnum,$data,$agenes,$bgenes,$gene2pos,$a2b,$mingF) = @_;

  my (%i2rec,%i2data);
  for(my $i=0;$i<@$data;$i+=1) {
    if(!exists($i2rec{$i})) {
      my @data1 = split(/\t/,$data->[$i]);
      for(my $j=$i+1;$j<@$data;$j+=1) {
        if(!exists($i2rec{$j})) {
          my @data2 = split(/\t/,$data->[$j]);
          if($data1[1] eq $data2[1] && $data1[5] eq $data2[5]) {
            my $dist1 = $data2[2] - $data1[3] -1;
            my $dist2 = 0;
            my $flag = 0;
            if($data1[7]<$data2[6]) {
              $dist2 = $data2[6] - $data1[7] -1;
              if($dist1<$spnum && $dist2<$spnum) {
              #if($dist1<$spnum) {
                if($data1[2]>$data2[2]) {
                  $data1[2] = $data2[2];
                }
                if($data1[3]<$data2[3]) {
                  $data1[3] = $data2[3];
                }
                $data1[7] = $data2[7];
                $flag = 1;
              }
            }
            elsif($data1[6]>$data2[7]) {
              $dist2 = $data1[6] - $data2[7] -1;
              if($dist1<$spnum && $dist2<$spnum) {
              #if($dist1<$spnum) {
                if($data1[2]>$data2[2]) {
                  $data1[2] = $data2[2];
                }
                if($data1[3]<$data2[3]) {
                  $data1[3] = $data2[3];
                }
                $data1[6] = $data2[6];
                $flag = 1;
              }
            }
            elsif($dist1<$spnum) {
              if($data1[2]>$data2[2]) {
                $data1[2] = $data2[2];
              }
              if($data1[3]<$data2[3]) {
                $data1[3] = $data2[3];
              }
              if($data1[6]>$data2[6]) {
                $data1[6] = $data2[6];
              }
              if($data1[7]<$data2[7]) {
                $data1[7] = $data2[7];
              }
              $flag = 1;
            }
            if($flag==1) {
              $i2rec{$j} = "Y";
              $data1[4] += $data2[4];
              $data1[8] += $data2[8];
              $j = $i+1;                       #-recheak from head
            }
          }
        }
      }
      $i2data{$i} = "$data1[1]\t$data1[2]\t$data1[3]\t$data1[4]\t$data1[5]\t$data1[6]\t$data1[7]\t$data1[8]";
    }
  }

  open(my $FW1,">$out2");
  open(my $FW2,">$out3");
  my $index = 1;
  foreach my $ii (sort {$a<=>$b} keys %i2data) {
    my @data = split(/\t/,$i2data{$ii});
    if($data[3]>$mingF) {
      print $FW1 "$index\t$i2data{$ii}\n";

      my ($agene1,$agene2) = ($agenes->[$data[1]],$agenes->[$data[2]]);
      my @agpos1 = split(/\t/,$gene2pos->{$agene1});
      my @agpos2 = split(/\t/,$gene2pos->{$agene2});
      if($agpos1[0] ne $data[0] || $agpos2[0] ne $data[0]) {
        print "What??? $data[0]\t$agpos1[0]\t$agpos2[0]\n";
      }
      my ($bgene1,$bgene2) = ($bgenes->[$data[5]],$bgenes->[$data[6]]);
      my @bgpos1 = split(/\t/,$gene2pos->{$bgene1});
      my @bgpos2 = split(/\t/,$gene2pos->{$bgene2});
      if($bgpos1[0] ne $data[4] || $bgpos2[0] ne $data[4]) {
        print "What??? $data[4]\t$bgpos1[0]\t$bgpos2[0]\n";
      }

      my ($bag1,$bag2) = ($a2b->{$agene1},$a2b->{$agene2});
      my @bagp1 = split(/\t/,$gene2pos->{$bag1});
      my @bagp2 = split(/\t/,$gene2pos->{$bag2});
      if($bagp1[0] ne $data[4] || $bagp2[0] ne $data[4]) {
        print "What??? $data[4]\t$bagp1[0]\t$bagp2[0]\n";
      }

      if($bagp1[1] < $bagp2[1]) {
        print $FW2 "FragA$index\t$data[0]\t$agpos1[1]\t$agpos2[2]\tFragB$index\t$data[4]\t$bgpos1[1]\t$bgpos2[2]\n";
      }
      else {
        print $FW2 "FragA$index\t$data[0]\t$agpos2[2]\t$agpos1[1]\tFragB$index\t$data[4]\t$bgpos1[1]\t$bgpos2[2]\n";
      }
      $index += 1;
    }
  }
  close($FW1);
  close($FW2);
}

sub pickSeg {
  my ($out,$a2b,$agenes2i,$bgenes2i,$agenes,$bgenes,$gene2pos,$itnum,$data) = @_;

  my $flag = 0;
  my $index = 1;
  my ($cag,$cbg) = (0,0);
  my ($startai,$stopai,$startbi,$stopbi,$pchra,$pchrb,$agene,$bgene,$posa,$posb,$achr,$bchr,$agi,$bgi,$ritnum);
  open(my $FW,">$out");
  for(my $i=0;$i<@$agenes;$i+=1) {
    ($agene,$bgene) = ($agenes->[$i],$a2b->{$agenes->[$i]});
    ($posa,$posb) = ($gene2pos->{$agene},$gene2pos->{$bgene});
    my @pas = split(/\t/,$posa);
    my @pbs = split(/\t/,$posb);
    ($achr,$bchr) = ($pas[0],$pbs[0]);
    ($agi,$bgi) = ($agenes2i->{$agene},$bgenes2i->{$bgene});

    if($flag==1) {
      $ritnum = 0;
      #if($startbi<=$stopbi) {
        if($bgi<$startbi) {
          $ritnum = $startbi - $bgi - 1;
        }
        elsif($bgi>$stopbi) {
          $ritnum = $bgi - $stopbi - 1;
        }
      #}
      #elsif($startbi>$stopbi) {
      #  if($bgi>$startbi) {
      #    $ritnum = $bgi - $startbi - 1;
      #  }
      #  elsif($bgi<$stopbi) {
      #    $ritnum = $stopbi - $bgi - 1;
      #  }
      #}
      

      if($pchra ne $achr || $pchrb ne $bchr || $ritnum > $itnum) {
        $data->[$index-1] = "$index\t$pchra\t$startai\t$stopai\t$cag\t$pchrb\t$startbi\t$stopbi\t$cbg";
        print $FW "$data->[$index-1]\t**$pchra ne $achr || $pchrb ne $bchr || $ritnum > $itnum\n";
        $index += 1;
        $cag = 0;
        $cbg = 0;
        $flag = 0;
      }
      else {
        $stopai = $agi;
        if($bgi<$startbi) {
          $startbi = $bgi;
        }
        elsif($bgi>$stopbi) {
          $stopbi = $bgi;
        }
        $cag += 1;
        $cbg += 1;
      }
    }

    if($flag==0) {
      $startai = $agi;
      $stopai = $agi;
      $startbi = $bgi;
      $stopbi = $bgi;
      $pchra = $achr;
      $pchrb = $bchr;
      $cag = 1;
      $cbg = 1;
      $flag = 1;
    }
  }

  $data->[$index-1] = "$index\t$pchra\t$startai\t$stopai\t$cag\t$pchrb\t$startbi\t$stopbi\t$cbg";
  print $FW "$data->[$index-1]\t**$pchra ne $achr || $pchrb ne $bchr || $ritnum > $itnum\n";


  close($FW);
}

sub processSyn {
  my ($in,$out) = @_;

  open(my $FW,">$out");
  open(my $FR,$in);
  while($_=<$FR>) {
    chomp;
    my @temp = split(/\t/);
    $temp[1] =~ s/^SL4\.0//;
    $temp[1]=~/^([^\d]+)(\d+)$/;     #-it depends--
    my $chr = int($2);
    print $FW "a:$temp[0]\ta:$1";
    if($chr<10) {
      print $FW "00$chr";
    }
    elsif($chr<100) {
      print $FW "0$chr";
    }
    else {
      print $FW "$chr";
    }
    print $FW "\t$temp[2]\t$temp[3]\t";

    $temp[5]=~/^([^\d]+)(\d+)$/;
    $chr = int($2);
    print $FW "b:$temp[4]\tb:$1";
    if($chr<10) {
      print $FW "00$chr";
    }
    elsif($chr<100) {
      print $FW "0$chr";
    }
    else {
      print $FW "$chr";
    }
    print $FW "\t$temp[6]\t$temp[7]\n";
  }
  close($FR);
  close($FW);
}

sub readSyn {
  my ($syn,$a2b,$agenes2i,$bgenes2i,$agenes,$bgenes,$gene2pos,$chrM) = @_;

  open(my $FR,$syn);
  my $index = 0;
  while($_=<$FR>) {
    #if(/$chrM/ && !/NW/) {                    # it depends---
    if(!/Scaffold/) {                    # it depends---
      chomp;
      my @temp = split(/\t/);
      $a2b->{$temp[0]} = $temp[4];
      $agenes2i->{$temp[0]} = $index;
      $agenes->[$index] = $temp[0];
      $index += 1;
      
      $gene2pos->{$temp[0]} = $temp[1]."\t".$temp[2]."\t".$temp[3];
      $gene2pos->{$temp[4]} = $temp[5]."\t".$temp[6]."\t".$temp[7];
    }
  }
  close($FR);

  system("sort -k6,6 -k7,7n $syn > $syn.sort");
  open($FR,$syn.".sort");
  $index = 0;
  while($_=<$FR>) {
    #if(/$chrM/ && !/NW/) {                    # it depends---
    if(!/Scaffold/) {                    # it depends---
      chomp;
      my @temp = split(/\t/);
      $bgenes2i->{$temp[4]} = $index;
      $bgenes->[$index] = $temp[4];
      $index += 1;
    }
  }
  close($FR);
  system("rm $syn.sort");
}


