perl/etc/gcmtdiv.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "usage: $0 <file prefix>\n" unless @ARGV == 1;
   6
   7 my (@cnt,@ref,@div);
   8 open CNT,'<',$ARGV[0].'.sam.gc.stat' or die $!;
   9 open REF,'<',$ARGV[0].'.ref.gc.stat' or die $!;
  10
  11 while(<CNT>) { last if /^#COL=higher/; }
  12 while(<CNT>) {
  13         chomp;
  14         next if /^#/;
  15         #last if /^$/;
  16         my @t=split /\t/,$_;
  17         my $id=shift @t;
  18         $cnt[$id]=[@t];
  19 }
  20 close CNT;
  21
  22 while(<REF>) { last if /^#COL=higher/; }
  23 while(<REF>) {
  24         chomp;
  25         next if /^#/;
  26         my @t=split /\t/,$_;
  27         my $id=shift @t;
  28         $ref[$id]=[@t];
  29 }
  30 close REF;
  31
  32 $cnt[$_][$_] *= 2  for (0..100);
  33 $ref[$_][$_] *= 2  for (0..100);
  34
  35 sub dePepperC($$) {
  36         my ($h,$range)=@_;
  37         return 0 if $range > 2;
  38         my ($points,$sum)=(0,0);
  39         for my $i ($h-$range .. $h+$range) {
  40                 next unless $div[$i][$i];
  41                 next if $i==$h;
  42                 ++$points;
  43                 $sum+=$div[$i][$i];
  44         }
  45         if ($points >= 2) {
  46                 for my $i ($h-$range .. $h+$range) {
  47                         print $div[$i][$i],' ';
  48                 }
  49                 print "\n",$sum / $points,"\n",'*' x 75,"\n";
  50                 return $sum / $points;
  51         } else {
  52                 return &dePepperC($h,++$range);
  53         }
  54 }
  55 sub dePepperA($$) {
  56         my ($h,$l)=@_;
  57         my $range=1;
  58         my ($points,$sum)=(0,0);
  59         for my $i ($h-$range .. $h+$range) {
  60                 for my $j ($l-$range .. $l+$range) {
  61                         next unless $div[$i][$j];
  62                         next if $i==$h and $j==$l;
  63                         ++$points;
  64                         $sum+=$div[$i][$j];
  65                 }
  66         }
  67         if ($points >= 2) {
  68                 for my $i ($h-$range .. $h+$range) {
  69                         for my $j ($l-$range .. $l+$range) {
  70                                 print $div[$i][$j],' ';
  71                         }
  72                         print "\n";
  73                 }
  74                 print $sum / $points,"\n",'=' x 75,"\n";
  75                 return $sum / $points;
  76         } else {
  77                 return 0;
  78         }
  79 }
  80 sub dePepper($$$){
  81         my ($h,$l,$range)=@_;
  82         my ($points,$sum)=(0,0);
  83         for my $i ($h-$range .. $h+$range) {
  84                 for my $j ($l-$range .. $l+$range) {
  85                         next unless $ref[$i][$j];
  86                         next if $i==$h and $j==$l;
  87                         ++$points;
  88                         $sum+=$ref[$i][$j];
  89                 }
  90         }
  91         if ($points >= 2) {
  92                 for my $i ($h-$range .. $h+$range) {
  93                         for my $j ($l-$range .. $l+$range) {
  94                                 print $ref[$i][$j],' ';
  95                         }
  96                         print "\n";
  97                 }
  98                 print $sum / $points,"\n";
  99                 return $sum / $points;
 100         } else {
 101                 return &dePepper($h,$l,++$range);
 102         }
 103 }
 104
 105 my $refvalue=0;
 106 for my $h (0..100) {
 107         for my $l (0..100) {
 108                 #print "$h,$l\tCnt:$cnt[$h][$l]\tRef:$ref[$h][$l]\n";
 109                 #next unless $cnt[$h][$l];
 110                 $refvalue = $ref[$h][$l];
 111                 unless ($refvalue) {
 112                         if ($cnt[$h][$l]) {
 113                                 #$ref[$h][$l]=-1;
 114                                 $refvalue = &dePepper($h,$l,1);
 115                                 print '-' x 75,"\n";
 116                         } else {$refvalue=1;}
 117                 }
 118                 $div[$h][$l] = $cnt[$h][$l] / $refvalue;
 119         }
 120 }
 121 my @divOut;
 122 for my $h (0..100) {
 123         for my $l (0..100) {
 124                 if ($h<$l) {
 125                         $divOut[$h][$l] = 0;
 126                 } elsif ($div[$h][$l]) {
 127                         $divOut[$h][$l] = $div[$h][$l];
 128                 } else {
 129                         $divOut[$h][$l] = &dePepperA($h,$l);
 130                         if ($h == $l and $divOut[$h][$l] == 0) {
 131                                 $divOut[$h][$l] = &dePepperC($h,1);
 132                         }
 133                 }
 134         }
 135 }
 136
 137 open OUT,'>',$ARGV[0].'.div.gc.stat' or die $!;
 138 print OUT "#COL=higher_GC%_read, ROW=lower_GC%_read\n";
 139 print OUT join("\t",'#',0..100),"\n";
 140 print OUT "$_\t",join("\t",@{$divOut[$_]}),"\n" for (0..100);
 141 close OUT;
 142 #print "$_: [",scalar @{$cnt[$_]},']',join(' ',@{$cnt[$_]}),"\n" for (0 .. $#cnt);
 143 #print "$_: [",scalar @{$ref[$_]},']',join(' ',@{$ref[$_]}),"\n" for (0 .. $#ref);
 144