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

   1 use strict;
   2 use warnings;
   3
   4 my $qc_dir = shift;
   5 my $output = shift;
   6 open OUT,">$output" or die($!);
   7
   8 $qc_dir =~ s/\/$//g;
   9 my @path = split /\//,$qc_dir;
  10 my $now = $path[scalar @path - 2];
  11 my @useless = splice @path,scalar @path - 2;
  12 my $total_dir = join "/",@path;
  13 my @content = readpipe("ls -F $total_dir");
  14
  15 my %chips;
  16 for (@content){
  17         chomp($_);
  18         if ($_ =~ /^\d+CG\S+\/$/ or $_ =~ /^\d+CG\S+\@$/){
  19                 $_ =~ s/\/$//;
  20                 $_ =~ s/\@$//;
  21                 my ($date) = $_ =~ /^(\d+)CG/;
  22                 $chips{$_} = $date;
  23         }
  24 }
  25
  26 my @dirs = ($now);
  27 if ($now =~ /^\d+CG/){
  28         my ($nowtime) = $now =~ /^(\d+)CG/;
  29         foreach my $dir (sort {$chips{$b}<=>$chips{$a}} keys %chips){
  30                 next if ($dir eq $now);
  31                 next unless (-e "$total_dir/$dir/6record");
  32                 if ($chips{$dir} eq $nowtime){
  33                         push @dirs,$dir;
  34                 }elsif ($chips{$dir} < $nowtime && scalar @dirs < 3){
  35                         push @dirs,$dir;
  36                 }
  37         }
  38 }
  39
  40 my (%nFs,%aFs,%store);
  41 foreach my $take_dir (@dirs){
  42         my $list = "$total_dir/$take_dir/family.lst";
  43         open LI,"<$list" or die($!);
  44         while (my $info = <LI>){
  45                 chomp($info);
  46                 next unless ($info =~ /\S+/);
  47                 my ($M,$F,$C) = split /\s+/,$info;
  48                 $aFs{$M} = "p$C.M";
  49                 $aFs{$F} = "p$C.F";
  50                 if ($take_dir eq $now){
  51                         $nFs{$M} = "p$C.M";
  52                         $nFs{$F} = "p$C.F";
  53                 }
  54                 $store{$M} = "$total_dir/$take_dir/4tsv";
  55                 $store{$F} = "$total_dir/$take_dir/4tsv";
  56         }
  57         close LI;
  58 }
  59
  60 my @now_fathers = sort keys %nFs;
  61 my @all_fathers = sort keys %aFs;
  62 print OUT "\t";
  63 print OUT join("\t",@all_fathers),"\n";
  64 foreach my $F1 (@now_fathers){
  65         print OUT "$F1\t";
  66         my @value;
  67         foreach my $F2 (@all_fathers){
  68                 my (%one,%two);
  69                 &get_geno(\%one,$F1);
  70                 &get_geno(\%two,$F2);
  71                 my ($err,$total) = (0,0);
  72                 foreach my $locus (keys %one){
  73                         if (defined $two{$locus}){
  74                                 $total++;
  75                                 unless ($one{$locus} eq $two{$locus}){
  76                                         $err++;
  77                                 }
  78                         }
  79                 }
  80                 if ($total > 0){
  81                         push @value,$err / $total;
  82                 }else{
  83                         push @value,"NA";
  84                 }
  85         }
  86         print OUT join("\t",@value),"\n";
  87 }
  88 close OUT;
  89
  90 ##########################################
  91 sub get_geno {
  92         my ($hash,$tag) = @_;
  93         my $file = "$store{$tag}/$aFs{$tag}.tsv";
  94         open TE,"<$file" or die($!);
  95         while (<TE>){
  96                 chomp;
  97                 my @data = split /\t/,$_;
  98                 next unless ($data[3] > 100);
  99                 my @geno = splice @data,4;
 100
 101                 my %check;
 102                 for (@geno){
 103                         my @info = split /;/,$_;
 104                         my @alleles = split /\//,$info[0];
 105                         my @sort = sort @alleles;
 106                         my $allele = join "/",@sort;
 107                         $check{$allele}++;
 108                 }
 109                 foreach my $allele (keys %check){
 110                         if ($check{$allele} == scalar @geno){
 111                                 $hash->{$data[0]} = $allele;
 112                         }
 113                 }
 114         }
 115         close TE;
 116 }