perl/etc/WoodyMiaoLin/Tiger/maternity.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 open I, "-|", "bcftools view tigers.bcgv.bcf";
   6
   7 my (@isMom, @noMom);
   8 while (<I>) {
   9   next if /^#/;
  10   my @line = split /\t/;
  11   next if $line[7] =~ /^INDEL/;
  12   next if $line[5] < 20;
  13   my @JHHgt = split /:/, $line[9];
  14   next if $JHHgt[4] < 20;
  15   my ($Ja1, $Ja2) = split /\//, $JHHgt[0];
  16   foreach (10..26) {
  17     my @ANYgt = split /:/, $line[$_];
  18     next if $ANYgt[4] < 20;
  19     my ($Aa1, $Aa2) = split /\//, $ANYgt[0];
  20     my $a;
  21     ++$a if $Aa1 eq $Ja1;
  22     ++$a if $Aa1 eq $Ja2;
  23     ++$a if $Aa2 eq $Ja1;
  24     ++$a if $Aa2 eq $Ja2;
  25     $_ -= 10;
  26     if ($a) {
  27       ++$isMom[$_];
  28     } else {
  29       ++$noMom[$_];
  30     }
  31   }
  32 }
  33
  34 close I;
  35 open O, ">", "maternity.log";
  36
  37 print O "#\tGZXJ03\tGZXJ05\tGZXJ26\tGZXJ27\tGZXJ28\tGZXJ29\tGZXJ30\tGZXJ33\tBHX011\tBHX019\tGZXJ04\tGZXJ06\tGZXJ31\tGZXJ32\tGZXJ36\tGZXJ37\tGZXJ38\n";
  38 print O "isMom\t";
  39 print O join "\t", @isMom;
  40 print O "\nnoMom\t";
  41 print O join "\t", @noMom;
  42 print O "\n";
  43
  44 close O;