projects/recombineMap/dynResolve.pl

   1 #!/bin/env perl
   2 =pod
   3 Author: Hu Xuesong @ BIOPIC <galaxy001@gmail.com>
   4 Version: 1.0.0 @ 20120412
   5 =cut
   6 use strict;
   7 use warnings;
   8 #use Data::Dump qw(ddx);
   9
  10 die "Usage: $0 <in file> <out file>\n" if @ARGV<2;
  11 my $in=shift;
  12 my $out=shift;
  13
  14 my (@Indi,@Parents);
  15 my ($indi,$site)=(0,0);
  16
  17 open I,'<',$in or die;
  18 while (<I>) {
  19         next if /^#/;
  20         chomp;
  21         my ($parent,undef,$id,@GT)=split /\t/;
  22         $Indi[$id] = \@GT;      # @Indi[SNP][0..$indi-1]
  23         push @Parents,$parent;
  24         die if $indi && $indi != @GT;
  25         $indi = scalar @GT;
  26 }
  27 $site = $#Indi; # count of data lines, 1st is 0. @Indi[1..SNP][0..$indi-1]
  28 push @{$Indi[0]},0 for (0 .. $indi-1);  # also set @Indi[SNP==0][0..$indi-1]=0
  29 print "Loaded: $site SNP sites * $indi Sperms.\n";
  30
  31 my @a=split //,('0' x $site);
  32 my @b=split //,('1' x $site);
  33 my (@matrix,@path);     # 2D array in perl can be malloced auto.ly. @matrix[0..SNP][0..SNP]
  34 my ($i,$j,$sc,$mstep);
  35 my %COLOR=('1'=>'32','2'=>'33','3'=>'36','0'=>'0');
  36
  37 sub scoring($$$) {      # Recombin_Count >> 0, so 0 is reserved for NULL.
  38         my ($i,$LR,$lastGrid)=@_;
  39         my ($lasti,$lastj)=@$lastGrid;
  40         my $lastLR;
  41         if ($path[$lasti][$lastj] == 1) {
  42                 $lastLR = 0;
  43         } else { $lastLR = 1; }
  44         my $RC=0;
  45         for my $n (0..$indi-1) {
  46                 if ( ($LR^$Indi[$mstep][$n]) != ($lastLR^$Indi[$mstep-1][$n]) ) {       # remember the higher precedence of '!=' to '^'
  47                         ++$RC;
  48                 }
  49         }
  50         return $matrix[$lasti][$lastj]+$RC;
  51 }
  52
  53 # Initialization
  54 $matrix[0][0]=0;
  55 $path[0][0]=0;
  56 for ($i=0;$i<=$site;$i++) {
  57         for ($j=0;$j<=$site;$j++) {
  58                 $matrix[$i][$j] = 0;
  59         }
  60 }
  61 #for ($j=1;$j<=$#b;$j++) {$matrix[0][$j] = $j*$INDEL;$path[0][$j]=3;}
  62 =pod
  63  ij                                    j => R
  64          00, 0, -      01, 1, R      02, 2, RR     03, 3, RRR    04, 4, RRRR
  65          10, 1, L      11, 2, LR/RL  12, 3, NNN    13, 4, NNNN   (14, 5)
  66  i => L  20, 2, LL     21, 3, NNN    22, 4, NNNN   (23, 5)       (24, 6)
  67          30, 3, LLL    31, 4, NNNN   (32, 5)       (33, 6)       (34, 7)
  68          40, 4, LLLL   (41, 5)       (42, 6)       (43, 7)       (44, 8)
  69 =cut
  70 for ($mstep=1;$mstep<=$site;$mstep++) { # starts from 01 & 10.
  71         for ($i=0;$i<=$mstep;$i++) {
  72                 $j = $mstep-$i; # only half matrix needed.
  73                 my ($sc1,$sc2)=(0,0);   # better to set NULL as -1.
  74                 if ($j>0) {     # has Up, can be L
  75                         $sc1=scoring($i,0,[$i,$j-1]);
  76                 }
  77                 if ($i>0) {     # has Left, can be R
  78                         $sc2=scoring($i,1,[$i-1,$j]);
  79                 }
  80                 if ($sc1) {
  81                         if ($sc2) {     # both $sc1 and $sc2
  82                                 if ($sc1 < $sc2) {
  83                                         $path[$i][$j]=1;
  84                                         $sc = $sc1;
  85                                 } elsif ($sc1 == $sc2) {
  86                                         $path[$i][$j]=3;        # half L => 0.5, no difference.
  87                                         $sc = $sc1;
  88                                 } else {        # $sc1 > $sc2
  89                                         $path[$i][$j]=2;
  90                                         $sc = $sc2;
  91                                 }
  92                         } else {        # only $sc1
  93                                 $path[$i][$j]=1;
  94                                 $sc = $sc1;
  95                         }
  96                 } elsif ($sc2) {        # only $sc2
  97                         $path[$i][$j]=2;
  98                         $sc = $sc2;
  99                 } else {        # $sc1==$sc2==0
 100                         #die "Impossible for $i==$j==0. $sc1,$sc2.";
 101                         $path[$i][$j]=0;        # for path, 00 <=> 11
 102                         $sc = 0;
 103                 }
 104                 $matrix[$i][$j] = $sc;
 105         }
 106 }
 107 my ($min,$mini,$minj)=($matrix[0][$site],0,$site);
 108 for ($i=1;$i<=$site;$i++) {     # init value is from $i==0
 109         $j = $site-$i;
 110         if ($min > $matrix[$i][$j]) {
 111                 $min=$matrix[$i][$j];
 112                 ($mini,$minj)=($i,$j);
 113         }
 114 }
 115 # Traceback
 116 my ($sumxorP,@resultP,@xorP)=(0);
 117 while ($mini>=0 and $minj>=0 and $mini+$minj) {
 118         push @resultP,$path[$mini][$minj];
 119         if ($path[$mini][$minj] & 1) {
 120                 --$minj;
 121         } else { --$mini; }
 122 }
 123 @resultP = reverse @resultP;
 124 for (0..$#resultP) {
 125         my $t = ($resultP[$_]^$Parents[$_])&1;
 126         push @xorP,$t;
 127         $sumxorP += $t;
 128 }
 129
 130 print "Parents: ",join(',',@Parents),"\n";
 131 print "Results: ",join(',',@resultP),"\n";
 132 print "bitXORs: ",join(',',@xorP),"\nSumbXOR: $sumxorP\n";
 133 print '-' x 78,"\nColor:";
 134 for (sort keys %COLOR) {
 135         print "\033[",$COLOR{$_},";1m $_";
 136 }
 137 print "\033[0;0m\n";
 138
 139 open O,'>',$out or die;
 140 print O "Loaded: $site SNP sites * $indi Sperms.\n";
 141
 142 for ($i=0;$i<=$site;$i++) {
 143         for ($j=0;$j<=$site-$i;$j++) {
 144                 print "\033[",$COLOR{$path[$i][$j]},";1m";
 145                 printf("%4d ",$matrix[$i][$j]);
 146                 print O $matrix[$i][$j],',',$path[$i][$j],"\t";
 147         }
 148         print "\033[0;0m\n";
 149         print O "\n";
 150 }
 151
 152 print O '-' x 78,"\n";
 153 print O "Parents: ",join(',',@Parents),"\n";
 154 print O "Results: ",join(',',@resultP),"\n";
 155 print O "bitXORs: ",join(',',@xorP),"\nSumbXOR: $sumxorP\n";
 156 close O;