perl/etc/WoodyMiaoLin/Tiger/StatisticPstIcoverage.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4 use threads;
   5
   6
   7 sub SiteCover {
   8   my $file = shift;
   9   my %count;
  10   open my $in, "-|", "samtools view -f 64 -F 1796 $file" or die "Error opening $file: $!\n";
  11   while (<$in>) {
  12     my @align = split /\t/;
  13     my @match = $align[5] =~ /(\d+)M/g;
  14     my $sum_match = 0;
  15     $sum_match += $_ foreach @match;
  16     next if $sum_match < 64;
  17     @match = sort {$b <=> $a} @match;
  18     next if $match[0] < 30;
  19     if ($align[1] & 0x10) {
  20       next if $align[5] =~ /[NHP=X]/;
  21       next unless $align[5] =~ /(\d+)M$/;
  22       next unless $1 >= 5;
  23       my @del = $align[5] =~ /(\d+)D/g;
  24       my $sum_del = 0;
  25       $sum_del += $_ foreach @del;
  26       my $c = $align[3] + $sum_match + $sum_del - 1;
  27       $count{$align[2]}{$c}[0] = 1;
  28     } else {
  29       next unless $align[5] =~ /^(\d+)M/;
  30       next unless $1 >= 5;
  31       $count{$align[2]}{$align[3]+3}[1] = 1;
  32     }
  33   }
  34   close $in;
  35   my $t = localtime;
  36   print STDERR "$t\tCount file $file complete!\n";
  37   return \%count;
  38 }
  39
  40
  41 my %thread; #Creat a hash of threads.
  42 while (<>) {
  43   chomp;
  44   $thread{$_} = threads->new (\&SiteCover, $_);
  45 }
  46 my $t = localtime;
  47 print STDERR "$t\tDestribute threads completed!\n";
  48
  49
  50 my %all; #Threads return a hash of all individual's coverage.
  51 foreach (keys %thread) {
  52   $all{$_} = $thread{$_}->join;
  53 }
  54 $t = localtime;
  55 print STDERR "$t\tRead all file completed!\n";
  56
  57
  58 my (%stat_site, %stat_indiv); #Count sites' coverage and individuals' coverage.
  59 foreach my $i (keys %all) {
  60   $i =~ /^(\d\d)/;
  61   my $bit = 18 - $1;
  62   foreach my $s (keys %{$all{$i}}) {
  63     foreach my $c (keys %{$all{$i}{$s}}) {
  64       $stat_site{$s}{$c}[0] = 0 unless defined $stat_site{$s}{$c}[0];
  65       $stat_site{$s}{$c}[1] = 0 unless defined $stat_site{$s}{$c}[1];
  66       if (defined $all{$i}{$s}{$c}[0]) {
  67         ++$stat_indiv{$i}[0];
  68         $stat_site{$s}{$c}[0] = $stat_site{$s}{$c}[0] | ($all{$i}{$s}{$c}[0] <<= $bit);
  69       } else {
  70         $all{$i}{$s}{$c}[0] = 0;
  71       }
  72       if (defined $all{$i}{$s}{$c}[1]) {
  73         ++$stat_indiv{$i}[1];
  74         $stat_site{$s}{$c}[1] = $stat_site{$s}{$c}[1] | ($all{$i}{$s}{$c}[1] <<= $bit);
  75       } else {
  76         $all{$i}{$s}{$c}[1] = 0;
  77       }
  78       $all{$i}{$s}{$c}[2] = $all{$i}{$s}{$c}[0] | $all{$i}{$s}{$c}[1];
  79       $all{$i}{$s}{$c}[3] = $all{$i}{$s}{$c}[0] & $all{$i}{$s}{$c}[1];
  80       ++$stat_indiv{$i}[2] if $all{$i}{$s}{$c}[2];
  81       ++$stat_indiv{$i}[3] if $all{$i}{$s}{$c}[3];
  82     }
  83   }
  84 }
  85 $t = localtime;
  86 print STDERR "$t\tDo stastistics completed!\n";
  87
  88
  89 my @stat_time; #Print site coverage and count coverage of times.
  90 open Osite, ">", "PstIcoverage.lst";
  91 print Osite "#SCAFFOLD\tCOORDINATE\tREVERSE\tFORWARD\tUNION\tINTERSECTION\n";
  92 foreach my $s (sort keys %stat_site) {
  93   foreach my $c (sort {$a <=> $b} keys %{$stat_site{$s}}) {
  94     $stat_site{$s}{$c}[2] = $stat_site{$s}{$c}[1] | $stat_site{$s}{$c}[0];
  95     $stat_site{$s}{$c}[3] = $stat_site{$s}{$c}[1] & $stat_site{$s}{$c}[0];
  96     printf Osite "%s\t%d\t%018b\t%018b\t%018b\t%018b\n", $s, $c, @{$stat_site{$s}{$c}};
  97     foreach my $r (0..3) {
  98       my @n = split //, sprintf("%b", $stat_site{$s}{$c}[$r]);
  99       my $n = 0;
 100       $n += $_ foreach @n;
 101       ++$stat_time[$n][$r];
 102     }
 103   }
 104 }
 105 close Osite;
 106 $t = localtime;
 107 print STDERR "$t\tPrint to PstIcoverage.lst completed!\n";
 108
 109
 110 open Oindiv, ">", "PstIindivCoverage.lst"; #Print individual and time coverage.
 111 print Oindiv "#INDIV.\tREVERSE\tFORWARD\tUNION\tINTERSECTION\n";
 112 foreach my $i (sort keys %stat_indiv) {
 113   printf Oindiv "%s\t%d\t%d\t%d\t%d\n", $i, @{$stat_indiv{$i}};
 114 }
 115 print Oindiv "#TIME\tREVERSE\tFORWARD\tUNION\tINTERSECTION\tTIME\tREVERSE\tFORWARD\tUNION\tINTERSECTION\n";
 116 my @accum;
 117 foreach my $n1 (1..18) {
 118   foreach my $n2 (0..3) {
 119     foreach my $n3 ($n1..18) {
 120       $accum[$n1][$n2] += $stat_time[$n3][$n2];
 121     }
 122   }
 123   printf Oindiv "%s\t%d\t%d\t%d\t%d\t%s\t%d\t%d\t%d\t%d\n", ">=$n1", @{$accum[$n1]}, "=$n1", @{$stat_time[$n1]};
 124 }
 125 close Oindiv;
 126 $t = localtime;
 127 print STDERR "$t\tPrint to PstIindivCoverage.lst completed!\n";
 128
 129
 130 $t = localtime;
 131 print STDERR "$t\tDone!\n";
 132
 133