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

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 my ($list,$store,$outdir) = @ARGV;
   6
   7 my (%stat,%dep);
   8 my @fams;
   9 open LI,"<$list" or die($!);
  10 while (my $info = <LI>){
  11         chomp($info);
  12         next if ($info =~ /^#/);
  13         my ($M,$F,$C) = split /\s+/,$info;
  14         &output($C,"M",$M);
  15         &output($C,"F",$F);
  16 }
  17 close LI;
  18
  19 #########################################################
  20 sub output {
  21         my ($C,$tag,$id) = @_;
  22         open TE,"<$store/p$C.$tag.tsv" or die($!);
  23         open TOUT,">$outdir/$id.txt" or die($!);
  24         print TOUT "#Sample ID: $id\n";
  25         print TOUT "MarkerID\tGenoType\tQUAL\tACGTdepths\n";
  26
  27         my @AllBase = qw{A C G T};
  28         while (my $line = <TE>){
  29                 chomp($line);
  30                 my @data = split /\t/,$line;
  31                 next if ($data[3] eq '.' or $data[3] < 100);
  32                 my @tM = splice @data,4;
  33                 my %Dep;
  34                 my $depcheck = 1;
  35                 my @bases = split /,/,$data[2];
  36                 for (@tM){
  37                         my $depsum = 0;
  38                         my @Depinfo = split /[;,]/,$_;
  39                         for my $i (1..scalar @Depinfo - 1){
  40                                 if ($Depinfo[$i] eq '.'){$Depinfo[$i] = 0;}
  41                                 $Dep{$bases[$i - 1]} += $Depinfo[$i];
  42                                 $depsum += $Depinfo[$i];
  43                         }
  44                         if ($depsum <= 50){
  45                                 $depcheck *= 0;
  46                         }else{
  47                                 $depcheck *= 1;
  48                         }
  49                 }
  50                 next if ($depcheck == 0);
  51
  52                 my @sort = sort {$Dep{$b}<=>$Dep{$a}} keys %Dep;
  53                 my $genotype;
  54                 if (scalar @sort == 1){
  55                         $genotype = join "/",$sort[0],$sort[0];
  56                 }elsif (scalar @sort > 1 && $Dep{$sort[1]} <= $Dep{$sort[0]} * 0.01){
  57                         $genotype = join "/",$sort[0],$sort[0];
  58                 }elsif (scalar @sort > 1 && $Dep{$sort[1]} >= $Dep{$sort[0]} * 0.1){
  59                         my @rKeys = sort @sort[0,1];
  60                         $genotype = join "/",@rKeys;
  61                 }
  62
  63                 for (@AllBase){
  64                         unless (defined $Dep{$_}){
  65                                 $Dep{$_} = ".";
  66                         }
  67                 }
  68                 if (defined $genotype && $genotype =~ /[ACGT]+/){
  69                         print TOUT "$data[0]\t$genotype\t$data[3]\t$Dep{A},$Dep{C},$Dep{G},$Dep{T}\n";
  70                 }
  71         }
  72         close TE;
  73         close TOUT;
  74 }