trick/delly2treeomicsInput.pl at master · ruping/trick · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
use strict;

my $delly = shift;
my $sampleRename = shift;
my $pvaf = shift;
my $skipLC = shift;
if ($pvaf eq '') {
  $pvaf = "FALSE";
}
if ($skipLC eq '') {
  $skipLC = "TRUE";
}

my %sampleRename;
open IN, "$sampleRename";
while ( <IN> ) {
  chomp;
  my ($sn, $mapping) = split("\t", $_);
  $sampleRename{$sn} = $mapping;
}
close IN;

open IN, "$delly";
my %colnames;
my %colindex;
while ( <IN> ) {
  chomp;
  next if /^##/;
  if ($_ =~ /^#CHROM/) {
    s/^#//;
    my @cols = split("\t", $_);
    for( my $i = 0; $i <= $#cols; $i++ ) {
      $colnames{$cols[$i]} = $i;
      $colindex{$i} = $cols[$i];
    }
    my @samples;
    for (my $i = 9; $i < $#cols; $i++) {  #all samples
      #my $sample = $sampleRename{$colindex{$i}};
      my $sample = $cols[$i];
      if ($pvaf eq "TRUE") {
        $sample .= "mafc\t".$sample."refc\t".$sample."altc\t".$sample."d\t".$sample."mafa\t".$sample."ccf\t".$sample."ccfSD\t".$sample."time";
      }
      push(@samples, $sample);
    }
    my @header = qw(chr pos Change Gene);
    push(@header, @samples);
    printf("%s\n", join("\t", @header));
    if ($pvaf eq "FALSE"){
      printf STDERR ("%s\n", join("\t", @header));
    }
  } else {
    my @cols = split("\t", $_);
    my $chrom = $cols[$colnames{'CHROM'}];
    my $pos = $cols[$colnames{'POS'}];
    (my $change = $cols[$colnames{'ALT'}]) =~ s/[<>]//g;
    if ( $change =~ /chr/ ) {
      $change = 'TRA';
    }
    #$change = $cols[$colnames{'REF'}].'>A';
    #$change =~ s/N/A/;
    my $gene = $cols[$colnames{'ID'}];

    #IN: CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT
    my @formats = split(':', $cols[$colnames{'FORMAT'}]);
    my %formindex;
    for (my $f = 0; $f <= $#formats; $f++) {
      $formindex{$formats[$f]} = $f;
    }

    my @altc;
    my @depth;
    my @vaf;
    for (my $i = 9; $i < $#cols; $i++) {  #all samples
      my @sampleInfo = split(":", $cols[$i]);
      my $DR = $sampleInfo[$formindex{'DR'}];     #high-quality reference pairs
      my $DV = $sampleInfo[$formindex{'DV'}];     #high-quality variant pairs
      my $RR = $sampleInfo[$formindex{'RR'}];     #high-quality reference junction reads
      my $RV = $sampleInfo[$formindex{'RV'}];     #high-quality variant junction reads
      my $altc = $DV+$RV;
      my $depth = $DV+$RV+$DR+$RR;
      my $vaf = ($depth > 0)? sprintf("%.4f", $altc/$depth) : 0;
      push(@altc, $altc);
      push(@depth, $depth);
      push(@vaf, $vaf, $depth-$altc, $altc, $depth, 0, 0, 0, 0);
    }
    if ($pvaf eq "TRUE") {
      printf("%s\n", join("\t", $chrom, $pos, $change, $gene, @vaf));
    } else {
      printf("%s\n", join("\t", $chrom, $pos, $change, $gene, @altc));
      printf STDERR ("%s\n", join("\t", $chrom, $pos, $change, $gene, @depth));
    }
  }
}
close IN;