perl/etc/justonce/tmp/cutmsmcdat.pl

   1 #!/usr/bin/perl -w
   2
   3 use strict;
   4
   5 my ($fname) = @ARGV;
   6
   7 $fname =~ /\b(chr\w+)\b/;
   8 my $chrid = $1;
   9
  10 my @indcnt = 1 .. 9;
  11
  12 warn "($chrid,$fname)\n";
  13 my (@OutFH,@remained);
  14
  15 open IN,'<',$fname or die $!;
  16
  17 for (@indcnt) {
  18         my $t;
  19         open $t,'>',"dat/${chrid}.s$_.msmcdat" or die $!;
  20         $OutFH[$_] = $t;
  21         $remained[$_] = 0;
  22 }
  23
  24 while (<IN>) {
  25         chomp;
  26         my @dat = split /\t/,$_;
  27         #print "[$dat[-1]]\n";
  28         #next if length $dat[-1] != 54;
  29         next unless $dat[2] > 0;
  30         unless ($dat[-1] =~ /^[ATCG]+$/) {
  31                 warn "[!]$dat[-1].\n";
  32                 $remained[$_] += $dat[2] for @indcnt;
  33                 next;
  34         }
  35         for my $c (@indcnt) {
  36                 my $cnt = 2*$c;
  37                 my $str = substr $dat[-1],0,$cnt;
  38                 my %test;
  39                 ++$test{$_} for (split //,$str);
  40                 if ( (sort {$a<=>$b} values %test)[0] == $cnt ) {
  41                         $remained[$c] += $dat[2];
  42                         next;
  43                 }
  44                 my $FH = $OutFH[$c];
  45                 my $value3 = $dat[2] + $remained[$c];
  46                 $remained[$c] = 0;
  47                 print $FH join("\t",@dat[0,1],$value3,$str),"\n";
  48         }
  49 }
  50
  51 close IN;
  52
  53 close $OutFH[$_] for @indcnt;
  54
  55 __END__
  56 ls -1 20141118_high_coverage.chr*.msmcdat|xargs -n1 ./cutmsmcdat.pl
  57
  58 msmc2 --fixedRecombination -o out.s5.msmc2 chr1s5.msmcdat
  59
  60 msmc --fixedRecombination -o out1/s5 dat/chr*.s5.msmcdat
  61 msmc2 -o out2/s5 dat/chr*.s5.msmcdat
  62