perl/etc/WoodyMiaoLin/Tiger/GZXJ036parenthood.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 open I, "-|", "bcftools view tigers.bcgv.bcf";
   6 open NOTDAD, ">", "notDAD.vcf";
   7 open NOTMUM, ">", "notMUM.vcf";
   8 open LOG, ">", "parenthood.log";
   9
  10 my ($yDad, $nDad, $yMum, $nMum) = (0, 0, 0, 0);
  11 while (<I>) {
  12   next if /^#/;
  13   my @line = split /\t/;
  14   next if $line[7] =~ /^INDEL/;
  15   next if $line[5] < 999;
  16   my @JHH001 = split /:/, $line[9];
  17   next if $JHH001[4] < 99;
  18   my @BXH011 = split /:/, $line[18];
  19   next if $BXH011[4] < 99;
  20   my @GZXJ36 = split /:/, $line[25];
  21   next if $GZXJ36[4] < 99;
  22   my @aJ = split /\//, $JHH001[0];
  23   my @aB = split /\//, $BXH011[0];
  24   my @aG = split /\//, $GZXJ36[0];
  25   my ($GB, $GJ);
  26   foreach my $nG (0..1) {
  27     foreach my $nB (0..1) {
  28       ++$GB if $aG[$nG] eq $aB[$nB];
  29     }
  30     foreach my $nJ (0..1) {
  31       ++$GJ if $aG[$nG] eq $aJ[$nJ];
  32     }
  33   }
  34   if ($GB) {
  35     ++$yDad;
  36   } else {
  37     ++$nDad;
  38     foreach (0..8) {
  39       print NOTDAD $line[$_], "\t";
  40     }
  41     print NOTDAD $line[9], "\t", $line[18], "\t", $line[25], "\n";
  42   }
  43   if ($GJ) {
  44     ++$yMum;
  45   } else {
  46     ++$nMum;
  47     foreach (0..8) {
  48       print NOTMUM $line[$_], "\t";
  49     }
  50     print NOTMUM $line[9], "\t", $line[18], "\t", $line[25], "\n";
  51   }
  52 }
  53
  54 close I;
  55 close NOTDAD;
  56 close NOTMUM;
  57
  58 print LOG "yBXH011\tnBXH011\tyJHH001\tnJHH001\n";
  59 print LOG "$yDad\t$nDad\t$yMum\t$nMum\n";
  60
  61 close LOG;