projects/recombineMap/simulate.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 <Sperm count> <Site count>\n" if @ARGV<2;
  11 my $indi=shift;
  12 my $site=shift;
  13
  14 my $RateCMperMb=2.1;
  15 my $ChrLen=249250621;
  16 $ChrLen=24925062;
  17
  18 my $Rate = ($ChrLen/$site/1000000)*$RateCMperMb/100;
  19 my (@Indi,%S,@Parent);
  20 open O,'>',"gt_${indi}_${site}.gt" or die;
  21 open OR,'>',"gt_${indi}_${site}.dat" or die;
  22 print O "#Rate=$Rate, Distence=",int($ChrLen/$site),"\n#recombRate\tRecombined\tSNPid\t\@Genotype_of_${indi}\n";
  23 print OR "#Rate=$Rate, Distence=",int($ChrLen/$site),"\n#Parent\tRecombined\tSNPid\t\@Genotype_of_${indi}\n";
  24 push @{$Indi[0]},0 for (0 .. $indi-1);
  25 my $NewRate = $Rate;
  26 my $SumRate=0;
  27 for my $snp (1..$site) {
  28         my @t=();
  29         my $sumC=0;
  30         $Parent[$snp]=int(rand(2));
  31         $NewRate = ($Rate*$snp)-$SumRate;       # a feedback is better, EP(?)
  32         for my $i (0 .. $indi-1) {
  33                 my $x=(rand(1)>$NewRate)?0:1;
  34                 if ($x) {       # 1:recombine
  35                         push @t,($Indi[$snp-1][$i])?0:1;
  36                 } else {
  37                         push @t,($Indi[$snp-1][$i])?1:0;
  38                 }
  39                 $sumC += $x;
  40         }
  41         $SumRate += $sumC/$indi;
  42         #$NewRate = ($Rate*($snp+1))-$SumRat;
  43         print O join("\t",sprintf("% .6f",$NewRate),$sumC,$snp,@t),"\n";
  44         ++$S{$sumC};
  45         $Indi[$snp]=\@t;
  46         print OR join("\t",$Parent[$snp],$sumC,$snp,map {$_^$Parent[$snp]} @t),"\n";
  47 }
  48 close O;
  49 close OR;
  50
  51 my $t=0;
  52 my $cnt=$indi*$site;
  53 open O,'>',"gt_${indi}_${site}.stat" or die;
  54 print O "Rate=$Rate, Distence=",int($ChrLen/$site),"\ncM pre Mbp:",$RateCMperMb,"\nRecombine Times and Rate:\n";
  55 for my $c (sort {$a<=>$b} keys %S) {
  56         print O "$c\t$S{$c}\t",$c/$indi,"\n";
  57         $t += $c*$S{$c};
  58 }
  59 print O "Average Rate: ",$t/$cnt,"\n";
  60 print O "Average Recombine Rate: ",100*$t/$indi/($ChrLen/1000000)," (cM pre Mbp)\n\n";
  61
  62 my %Sgt;
  63 print O "GenoTypes:\n";
  64 for my $snp (1..$site) {
  65         for my $i (0 .. $indi-1) {
  66                 ++$Sgt{$Indi[$snp][$i]};
  67         }
  68 }
  69 print O "0: $Sgt{0}, 1: $Sgt{1}.\nRatio: ",$Sgt{0}/$Sgt{1},"\n";
  70 close O;