perl/etc/WoodyMiaoLin/Tiger/StatisticSNPcoverage.16tigers.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5
   6 my $t = localtime;
   7 print STDERR "$t\tStart!\n";
   8
   9
  10 my (@indiv, @stat_time, %nSNP, %iSNP);
  11 while (<>) {
  12   next if /^#/;
  13   chomp;
  14   my @line = split /\t/;
  15   next if $line[7] =~ /^INDEL/;
  16   next if $line[5] < 20;
  17   $line[7] =~ /;FQ=([0-9\-.]+)/;
  18   next unless $1 > 0;
  19   my $a = 0;
  20   my (@a, @b);
  21   foreach (9..26) {
  22     next if $_ == 24;
  23     next if $_ == 25;
  24     my @sample = split /:/, $line[$_];
  25     if ($sample[4] >= 20 and $sample[2] >= 1) {
  26       ++$a;
  27       push @a, $sample[0];
  28       push @b, $_;
  29     } else {
  30       next;
  31     }
  32   }
  33   next if $a < 6;
  34   next if (!grep {$_ ne $a[0]} @a);
  35   ++$stat_time[$a];
  36 #  $nSNP{$line[0]}{$line[1]} = $a;
  37 #  $iSNP{$line[0]}{$line[1]} = 0;
  38   foreach (@b) {
  39     $_ = 24 if $_ == 26;
  40     ++$indiv[$_];
  41 #    $iSNP{$line[0]}{$line[1]} |= (1 << (26-$_));
  42   }
  43 }
  44 $t = localtime;
  45 print STDERR "$t\tRead file completed!\n";
  46
  47
  48 #open OSNP, ">", "SNPcoverage.lst";
  49 #print OSNP "#SCAFFOLD\tCOORDINATE\tINVID.\tNO.\n";
  50 #foreach my $c (sort keys %nSNP) {
  51 #  foreach my $d (sort {$a <=> $b} keys %{$nSNP{$c}}) {
  52 #    printf OSNP "%s\t%d\t%018b\t%d\n", $c, $d, $iSNP{$c}{$d}, $nSNP{$c}{$d};
  53 #  }
  54 #}
  55 #close OSNP;
  56 #$t = localtime;
  57 #print STDERR "$t\tPrint to SNPCoverage.lst completed!\n";
  58
  59
  60 my @id = qw/01_JHH001 02_GZXJ03 03_GZXJ05 04_GZXJ26 05_GZXJ27 06_GZXJ28 07_GZXJ29 08_GZXJ30 09_GZXJ33 10_BHX011 11_BHX019 12_GZXJ04 13_GZXJ06 14_GZXJ31 15_GZXJ32 18_GZXJ38/;
  61 open Oindiv, ">", "SNPindivCoverage.16tigers.lst";
  62 print Oindiv "#INDIV.\tNO.\n";
  63 foreach (9..24) {
  64   print Oindiv "$id[$_-9]\t$indiv[$_]\n";
  65 }
  66 print Oindiv "#TIME\tNO.\tTIME\tNO.\n";
  67 my @accum;
  68 foreach my $e (6..16) {
  69   foreach ($e..16) {
  70     $accum[$e] += $stat_time[$_];
  71   }
  72   printf Oindiv "%s\t%d\t%s\t%d\n", ">=$e", $accum[$e], "=$e", $stat_time[$e];
  73 }
  74 close Oindiv;
  75 $t = localtime;
  76 print STDERR "$t\tPrint to SNPindivCoverage.lst completed!\n";
  77
  78
  79 $t = localtime;
  80 print STDERR "$t\tDone!\n";
  81
  82
  83
  84
  85