perl/etc/justonce/oyk/bin/contamination.MC.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 my ($list,$store,$output) = @ARGV;
   6
   7 my %stat;
   8 open LI,"<$list" or die($!);
   9 open OUT,">$output" or die($!);
  10 while (my $info = <LI>){
  11         chomp($info);
  12         my ($M,$F,$C) = split /\s+/,$info;
  13         my %hete = &get_hete($C);
  14         my $total = 0;
  15         my $mismatch = 0;
  16         open FC,"<$store/p$C.C.tsv" or die($!);
  17         while (my $line = <FC>){
  18                 chomp($line);
  19                 my @data = split /\t/,$line;
  20                 my @Bases = split /,/,$data[2];
  21                 next unless (defined $hete{$data[0]});
  22                 next if ($data[3] eq '.' or $data[3] < 100);
  23                 my @tmpInfo = splice @data,4;
  24                 my %Dep;
  25                 my $depcheck = 1;
  26                 my $sum = 0;
  27                 for (@tmpInfo){
  28                         my $depsum = 0;
  29                         my @Depinfo = split /[;,]/,$_;
  30                         for my $i (1..scalar @Depinfo - 1){
  31                                 if ($Depinfo[$i] eq '.'){$Depinfo[$i] = 0;}
  32                                 $Dep{$i - 1} += $Depinfo[$i];
  33                                 $depsum += $Depinfo[$i];
  34                         }
  35                         $sum += $depsum;
  36                         if ($depsum < 100){
  37                                 $depcheck *= 0;
  38                         }else{
  39                                 $depcheck *= 1;
  40                         }
  41                 }
  42                 next if ($depcheck == 0);
  43
  44                 $total++;
  45                 my @dKeys = sort { $Dep{$b} <=> $Dep{$a} } keys %Dep;
  46                 if (@dKeys>1 and $Dep{$dKeys[1]}  >= $Dep{$dKeys[0]} * 0.1){
  47                         my @rKeys = sort {$a<=>$b} @dKeys[0,1];
  48                         my $gt = join('/',$Bases[$rKeys[0]],$Bases[$rKeys[1]]);
  49                         unless ($gt eq $hete{$data[0]}){
  50                                 $mismatch++;
  51                                 print "$C\t$data[0]\t$gt\t$hete{$data[0]}\n";
  52                         }
  53                 }else{
  54                         $mismatch++;
  55 #                       print "$C\t$data[0]\n";
  56                 }
  57         }
  58         close FC;
  59         my $freq;
  60         if ($total == 0){
  61                 $freq = "NA";
  62         }else{
  63                 $freq = $mismatch/$total;
  64         }
  65         print OUT "p$C\t$mismatch/$total\t$freq\n";
  66 }
  67 close LI;
  68
  69 sub get_hete {
  70         my ($input) = @_;
  71         my %temp;
  72         open TE,"<$store/p$input.M.tsv" or die($!);
  73         while (my $line = <TE>){
  74                 chomp($line);
  75                 my @data = split /\t/,$line;
  76                 next if ($data[3] eq '.' or $data[3] < 100);
  77                 my @Depinfo = split /[;,]/,$data[4];
  78                 my @Bases = split /,/,$data[2];
  79                 my %Dep;
  80                 my $depsum = 0;
  81                 for my $i (1..scalar @Depinfo - 1){
  82                         if ($Depinfo[$i] eq '.'){$Depinfo[$i] = 0;}
  83                         $Dep{$i - 1} = $Depinfo[$i];
  84                         $depsum += $Depinfo[$i];
  85                 }
  86                 next if ($depsum <= 50);
  87                 my @dKeys = sort { $Dep{$b} <=> $Dep{$a} } keys %Dep;
  88                 if (@dKeys>1 and $Dep{$dKeys[1]}  >= $Dep{$dKeys[0]} * 0.1){
  89                         my @rKeys = sort {$a<=>$b} @dKeys[0,1];
  90                         my $gt = join('/',$Bases[$rKeys[0]],$Bases[$rKeys[1]]);
  91                         $temp{$data[0]} = $gt;
  92                 }else{
  93                         next;
  94                 }
  95         }
  96         close TE;
  97         return %temp;
  98 }