perl/etc/justonce/dephist.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 my %dat=(
   6 'Drone1' => 'd1.depth.xz',
   7 'Scouts1' => 's1.depth.xz',
   8 'Scouts2' => 's2.depth.xz',
   9 'Recruit2' => 'r2.depth.xz',
  10 'Recruit3' => 'r3.depth.xz',
  11 );
  12 my @id=sort keys %dat;
  13
  14 my $EffLen=199723619;
  15 my $Ncnt= 219629612 - $EffLen;
  16
  17 sub openfile($) {
  18     my ($filename)=@_;
  19     my $infile;
  20     if ($filename=~/.xz$/) {
  21             open( $infile,"-|","xz -dc $filename") or die "Error opening $filename: $!\n";
  22     } elsif ($filename=~/.gz$/) {
  23         open( $infile,"-|","gzip -dc $filename") or die "Error opening $filename: $!\n";
  24     } else {open( $infile,"<",$filename) or die "Error opening $filename: $!\n";}
  25     return $infile;
  26 }
  27
  28 my (%Hist,%Stack);
  29 sub dohist($$) {
  30         my ($id,$value)=@_;
  31         $value=0 if $value == 65535;
  32         if ($value <= 150) {
  33                 ++$Hist{$value}{$id};   # 0 .. 100
  34         } elsif ($value <= 1500) {
  35                 ++$Hist{200 + int(($value-100)/10)}{$id};       # 150~159->205, 1490~1499->339, 1500->340
  36         } else { ++$Hist{999}{$id}; }   # >1000
  37 }
  38 sub getV($){
  39         my $v=$_[0];
  40         if ($v<=150) {
  41                 return $v;
  42         } elsif ($v<=340) {
  43                 return '~'.(109 + ($v-200)*10);
  44         } else {
  45                 return '>1500';
  46         }
  47 }
  48
  49 for my $f (@id) {
  50         my $fh=openfile($dat{$f});
  51         $Hist{0}{$f} = -$Ncnt;
  52         while (<$fh>) {
  53                 next if /^>/;
  54                 chomp;
  55                 my @d = split /\s+/;
  56                 dohist($f,$_) for @d;
  57         }
  58 }
  59
  60 my %Sum;
  61 for my $v (sort {$b<=>$a} keys %Hist) {
  62         for my $id (@id) {
  63                 unless (exists $Hist{$v}{$id}) {
  64                         $Hist{$v}{$id}=0;
  65                 }
  66                 $Sum{$id} += $Hist{$v}{$id};
  67                 $Stack{$v}{$id} = $Sum{$id};
  68         }
  69 }
  70 open O,'>',"depth.hist";
  71 print O join("\t",'Depth',@id),"\n";
  72 for my $v (sort {$a<=>$b} keys %Hist) {
  73         print O getV($v),"\t";
  74         for my $id (@id) {
  75                 if (exists $Hist{$v}{$id}) {
  76                         print O $Hist{$v}{$id},', ',$Stack{$v}{$id},',',int(10000*$Stack{$v}{$id}/$EffLen)/10000,"\t";
  77                 }
  78         }
  79         print O "\n";
  80 }
  81 close O;