perl/etc/WoodyMiaoLin/Tiger/StatisticSNPcoverage.pl

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