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

   1 #!/usr/bin/env perl
   2 =pod
   3 Author: HU Xuesong @ BGI <huxuesong@genomics.org.cn>, LI Bowen <libowen@genomics.cn>
   4 Version: 1.0.0 @ 20180624
   5 =cut
   6 use strict;
   7 use warnings;
   8 use POSIX;
   9
  10 use FindBin qw($RealBin);
  11 if ($FindBin::VERSION < 1.51) {
  12         warn "[!]Your Perl is too old, thus there can only be ONE `bsuit` file in your PATH. [FindBin Version: $FindBin::VERSION < 1.51]\n\n"
  13 }
  14 FindBin::again();
  15 use lib "$RealBin/../";
  16
  17 my @tModes = qw(CHIP PCR);
  18 my %tMode = map { $_ => 1 } @tModes;
  19 my @tParentages = qw(DUO TRIO);
  20 my %tParentage = map { $_ => 1 } @tParentages;
  21
  22 my $thetMode = uc shift;
  23 unless (exists $tMode{$thetMode}) {
  24         die "[x]mode can only be:[",join(',',@tModes),"].\n";
  25 }
  26 my $thetParentage = uc shift;
  27 unless (exists $tParentage{$thetParentage}) {
  28         die "[x]Parentage can only be:[",join(',',@tParentages),"].\n";
  29 }
  30
  31 my $qc_dir = shift;
  32 my $output = shift;
  33 open OUT,">$output" or die($!);
  34
  35 $qc_dir =~ s/\/$//g;
  36 my @path = split /\//,$qc_dir;
  37 my $now = $path[scalar @path - 2];
  38 my @useless = splice @path,scalar @path - 2;
  39 my $total_dir = join "/",@path;
  40 my @content = readpipe("ls -F $total_dir");
  41
  42 my %chips;
  43 for (@content){
  44         chomp($_);
  45         if ($_ =~ /^\d+CG\S+\/$/ or $_ =~ /^\d+CG\S+\@$/){
  46                 $_ =~ s/\@$//;
  47                 $_ =~ s/\/$//;
  48                 my ($date) = $_ =~ /^(\d+)CG/;
  49                 $chips{$_} = $date;
  50         }
  51 }
  52
  53 my @dirs = ($now);
  54 if ($now =~ /^\d+CG/){
  55         my ($nowtime) = $now =~ /^(\d+)CG/;
  56         foreach my $dir (sort {$chips{$b}<=>$chips{$a}} keys %chips){
  57                 next unless (-e "$total_dir/$dir/6record");
  58                 next if ($dir eq $now);
  59                 if ($chips{$dir} eq $nowtime){
  60                         push @dirs,$dir;
  61                 }elsif ($chips{$dir} < $nowtime && scalar @dirs < 3){
  62                         push @dirs,$dir;
  63                 }
  64         }
  65 }
  66
  67 my @fams;
  68 my (%Fs,%store);
  69 foreach my $take_dir (@dirs){
  70         my $list = "$total_dir/$take_dir/family.lst";
  71         open LI,"<$list" or die($!);
  72         while (my $info = <LI>){
  73                 chomp($info);
  74                 my ($M,$F,$C) = split /\s+/,$info;
  75                 if ($take_dir eq $now){
  76                         push @fams,"p$C";
  77                 }
  78                 $Fs{$M} = "p$C.M";
  79                 $Fs{$F} = "p$C.F";
  80                 $store{$M} = "$total_dir/$take_dir/4tsv";
  81                 $store{$F} = "$total_dir/$take_dir/4tsv";
  82         }
  83         close LI;
  84 }
  85
  86 print join("\n",@fams),"\n";
  87 my @fathers = sort keys %Fs;
  88 print OUT "\tTotal\t";
  89 print OUT join("\t",@fathers),"\n";
  90 foreach my $family (@fams){
  91         my $Mfile = "$total_dir/$now/4tsv/$family.M.tsv";
  92         my $Cfile = "$total_dir/$now/4tsv/$family.C.tsv";
  93         print OUT "$family\t";
  94         my @outputs;
  95         my @totals;
  96         foreach my $father (@fathers){
  97                 my $total = 0;
  98                 my $mismatch = 0;
  99                 my $Ffile = "$store{$father}/$Fs{$father}.tsv";
 100                 system("mkdir -p $total_dir/$now/4tsv/temp") unless (-e "$total_dir/$now/4tsv/temp");
 101                 my @info = readpipe("perl $RealBin/oykn.pl $thetMode $thetParentage $Mfile $Ffile $Cfile $total_dir/$now/4tsv/temp/tm$family.$father");
 102                 foreach my $line (@info){
 103                         chomp($line);
 104                         next if ($line =~ /^#/);
 105                         my @data = split /\t/,$line;
 106                         next unless (defined $data[7]);
 107                         $total++;
 108                         if ($data[7] == 0.0001){
 109                                 $mismatch++;
 110                         }
 111                 }
 112                 push @outputs,$mismatch;
 113                 push @totals,$total;
 114         }
 115         my @sort = sort {$a<=>$b} @totals;
 116         print OUT "$sort[0]-$sort[scalar @sort - 1]\t";
 117         print OUT join("\t",@outputs),"\n";
 118 }
 119 close OUT;
 120
 121 system("rm -rf $total_dir/$now/4tsv/temp") if (-e "$total_dir/$now/4tsv/temp");