perl/etc/justonce/oyk/bin/heterozygosity.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 my ($list,$store,$output) = @ARGV;
   6
   7 open LI,"<$list" or die($!);
   8 open OUT,">$output" or die($!);
   9
  10 while (<LI>){
  11         chomp;
  12         my @ID = split /\s+/,$_;
  13         my @tag = qw{F M C};
  14         foreach my $i (0..2){
  15         my $id = $ID[$i];
  16         my $total = 0;
  17         my $hete = 0;
  18
  19         if ($tag[$i] eq 'C'){
  20                 my $file = "$store/r$ID[2].cpie";
  21                 open IN,"<$file" or die($!);
  22                 while (my $line = <IN>){
  23                         chomp($line);
  24                         next if ($line =~ /^#/);
  25                         my @data = split /[\t\/;]/,$line;
  26                         $total+=1;
  27                         unless ($data[10] eq $data[11]){
  28                                 $hete++;
  29                         }
  30                 }
  31         }else{
  32                 my $file = "$store/p$ID[2].$tag[$i].tsv";
  33                 open IN,"<$file" or die($!);
  34                 while (my $line = <IN>){
  35                         chomp($line);
  36                         my @data = split /\t/,$line;
  37                         next if ($data[3] eq '.' or $data[3] < 100);
  38                         my @tM = splice @data,4;
  39                         my %Dep;
  40                         my $depcheck = 1;
  41                         for (@tM){
  42                                 my $depsum = 0;
  43                                 my @Depinfo = split /[;,]/,$_;
  44                                 for my $i (1..scalar @Depinfo - 1){
  45                                         if ($Depinfo[$i] eq '.'){$Depinfo[$i] = 0;}
  46                                         $Dep{$i - 1} += $Depinfo[$i];
  47                                         $depsum += $Depinfo[$i];
  48                                 }
  49                                 if ($depsum <= 50){
  50                                         $depcheck *= 0;
  51                                 }else{
  52                                         $depcheck *= 1;
  53                                 }
  54                         }
  55                         next if ($depcheck == 0);
  56                         $total += 1;
  57
  58                         my @values = sort {$b<=>$a} values %Dep;
  59                         if (scalar @values > 1){
  60                                 if ($values[1] > $values[0] * 0.1){
  61                                         $hete++;
  62                                 }
  63                         }
  64                 }
  65                 close IN;
  66         }
  67         if ($total > 0){
  68                 print OUT "$id\t$hete/$total\t",$hete/$total,"\n";
  69         }else{
  70                 print OUT "$id\t$hete/$total\t","NA","\n";
  71         }
  72         }
  73 }
  74 close LI;
  75 close OUT;