tools/linkage/clinkage.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 unless (@ARGV > 0) {
   6     print "perl $0 <In_file> <Out_File>\n";
   7     exit 0;
   8 }
   9
  10 # choose Kosambi, which is slower than Haldane, for density markers
  11 my $Error=0;
  12 sub corRIL($) {
  13         my ($bigRth)=@_;
  14         my $smallr;
  15         $smallr=0.5*$bigRth/(1-$bigRth);
  16         if ($smallr >= 0.5) {
  17                 $Error=1;
  18                 return $smallr;
  19         } else {
  20                 $Error=0;
  21         }
  22         # applies the distance function
  23         my $RILdist = 25 * log((1 + 2 * $smallr) / (1 - 2 * $smallr));
  24         return $RILdist;
  25 }
  26
  27 open I,'<',$ARGV[0] or die $!;
  28 open O,'>',$ARGV[1] or die $!;
  29 my $t=<I>;
  30 $t=~s/\tPos\t/\tPos\tcM\tΔcM\tkB\/cM\t/;
  31 #print O $t;
  32 print O "#ChrID\tPos\tcM\tΔcM\tkB/cM\tR_RIL\n";
  33 my (@Pos,%Dat);
  34 my ($chr,$pos,$dat);
  35 while (<I>) {
  36         chomp;
  37         ($chr,$pos,$dat)=split /\t/;
  38         my @a=split / /,$dat;
  39         $Dat{$pos}=\@a;
  40         push @Pos,$pos;
  41 }
  42 my $LastcM=0;
  43 my $ratio;
  44 print O "$chr\t$Pos[0]\t$LastcM\t0\tNA\t0\n";   #,join(' ',@{$Dat{$Pos[0]}}),"\n";
  45 for my $i (1..$#Pos) {
  46         my $ary_a=$Dat{$Pos[$i-1]};
  47         my $ary_b=$Dat{$Pos[$i]};
  48         my ($cross,$cnt)=(0,0);
  49         for my $p (0..$#$ary_a) {
  50                 next if $ary_a->[$p] eq '-' or $ary_b->[$p] eq '-';
  51                 ++$cross if $ary_a->[$p] ne $ary_b->[$p];
  52                 ++$cnt;
  53         }
  54         my $cross0 = $cross/$cnt;
  55         $cross = &corRIL($cross0);
  56         $LastcM += $cross;
  57         $ratio=$cross? int(0.5+100*($Pos[$i]-$Pos[$i-1])/(1000*$cross)) : 'NA';
  58         $ratio='Error' if $Error == 1;
  59         print O "$chr\t$Pos[$i]\t",int(0.5+100*$LastcM)/100,"\t",int(0.5+100*$cross)/100,"\t"
  60                 ,$ratio,"\t",100*$cross0,"\n";  #,join(' ',@{$Dat{$Pos[$i]}}),"\n";
  61 }
  62
  63 close O;
  64 close I;
  65 __END__
  66 cat chrorder | while read a; do ./clinkage.pl $a.ccq cm$a.ccq;done &
  67
  68 cat chrorder | while read a; do ./linkage/clinkage.pl dat20101214/$a.marker dat20101214/$a.cm;done &