perl/etc/justonce/trireader.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4 use List::Util 1.26 qw(sum0);
   5 use Data::Dump qw(ddx);
   6
   7 my $fhead = 'snp/bam.head';
   8 my @files = qw[snp/SNP.out snp/tomor1.snp.out snp/tumour2.snp.out];
   9
  10 sub getDat($) {
  11         my $fh = $_[0];
  12         while(<$fh>) {
  13                 chomp;
  14                 my @dat = split /\t/;
  15                 next unless $dat[6] eq 'PASS';
  16                 my @dep1 = split /,/,$dat[9];
  17                 my @dep2 = split /,/,$dat[10];
  18                 my $sum1 = sum0(@dep1);
  19                 my $sum2 = sum0(@dep2);
  20                 next unless ($sum1+$sum2) > 20;
  21                 #return [@dat[0,1,7,6,9,10],$sum1+$sum2];
  22                 return [@dat[0,1,7,3]];
  23         }
  24         return ["\t",-1,"NN"];
  25 }
  26
  27 my (@ChrIDs,%ChrLen,%Chrank);
  28 open H,'<',$fhead or die $!;
  29 while (<H>) {
  30         chomp;
  31         my @dat = split /\t/;
  32         next if $dat[0] ne '@SQ';
  33         my $cid = (split /:/,$dat[1])[1];
  34         my $clen = (split /:/,$dat[2])[1];
  35         push @ChrIDs,$cid;
  36         $ChrLen{$cid} = $clen;
  37 }
  38 close H;
  39 %Chrank = map { $ChrIDs[$_], $_ + 1 } 0 .. $#ChrIDs;
  40 #ddx \%Chrank,\@ChrIDs,\%ChrLen;
  41
  42 my ($fha,$fhb,$fhc);
  43 open $fha,'<',$files[0] or die $!;
  44 open $fhb,'<',$files[1] or die $!;
  45 open $fhc,'<',$files[2] or die $!;
  46 <$fha>;<$fha>;<$fhb>;<$fhb>;<$fhc>;<$fhc>;
  47
  48 my @IDs = qw[A B C];
  49 my %FHs = (
  50         'A' => [$fha,getDat($fha)],
  51         'B' => [$fhb,getDat($fhb)],
  52         'C' => [$fhc,getDat($fhc)],
  53 );
  54 #print $FHs{'A'}->[1][0],"--\n";
  55
  56 open O,'>','res.tsv' or die $!;
  57 while (($FHs{'A'}->[1][0] ne "\t") and ($FHs{'B'}->[1][0] ne "\t") and ($FHs{'C'}->[1][0] ne "\t")) {
  58         #my %ChrIDs = map { $_ => $FHs{$_}->[1][0] } @IDs;
  59         my @aID = sort { $Chrank{$FHs{$b}->[1][0]} <=> $Chrank{$FHs{$a}->[1][0]} || $FHs{$b}->[1][1] <=> $FHs{$a}->[1][1] } @IDs;       # desc
  60         #ddx \%FHs,\@aID;
  61         for my $i (1 .. $#aID) {
  62                 while (
  63                         (($Chrank{$FHs{$aID[0]}->[1][0]} == $Chrank{$FHs{$aID[$i]}->[1][0]}) and ($FHs{$aID[0]}->[1][1] > $FHs{$aID[$i]}->[1][1]))
  64                         or ($Chrank{$FHs{$aID[0]}->[1][0]} > $Chrank{$FHs{$aID[$i]}->[1][0]})
  65                 ) {
  66                         $FHs{$aID[$i]}->[1] = getDat($FHs{$aID[$i]}->[0]);
  67                 }
  68         }
  69         if ($FHs{'A'}->[1][1] == $FHs{'B'}->[1][1] and $FHs{'A'}->[1][1] == $FHs{'C'}->[1][1]) {
  70                 #ddx \%FHs;
  71                 my $type='N';
  72                 if ( $FHs{'A'}->[1][2] eq $FHs{'B'}->[1][2] and $FHs{'A'}->[1][2] eq $FHs{'C'}->[1][2] ) {
  73                         $type='AAA';
  74                 } elsif ( $FHs{'A'}->[1][2] eq $FHs{'B'}->[1][2] and $FHs{'A'}->[1][2] ne $FHs{'C'}->[1][2] ) {
  75                         $type='AAB';
  76                 } elsif ( $FHs{'A'}->[1][2] ne $FHs{'B'}->[1][2] and $FHs{'A'}->[1][2] eq $FHs{'C'}->[1][2] ) {
  77                         $type='ABA';
  78                 } elsif ( $FHs{'A'}->[1][2] ne $FHs{'B'}->[1][2] and $FHs{'B'}->[1][2] eq $FHs{'C'}->[1][2] ) {
  79                         $type='ABB';
  80                 } elsif ( $FHs{'A'}->[1][2] ne $FHs{'B'}->[1][2] and $FHs{'A'}->[1][2] ne $FHs{'C'}->[1][2] and $FHs{'B'}->[1][2] ne $FHs{'C'}->[1][2] ) {
  81                         $type='ABC';
  82                 }
  83                 print O join("\t",$FHs{'A'}->[1][0],$FHs{'A'}->[1][1],$FHs{'A'}->[1][3],$type,join(',',$FHs{'A'}->[1][2],$FHs{'B'}->[1][2],$FHs{'C'}->[1][2])),"\n";
  84                 for my $f (@aID) {
  85                         $FHs{$f}->[1] = getDat($FHs{$f}->[0]);
  86                 }
  87         } else {
  88                 $FHs{'A'}->[1] = getDat($FHs{'A'}->[0]);
  89         }
  90 }
  91 close O;
  92 close $fha;close $fhb;close $fhc;