|
0
|
1 #!/usr/bin/env perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use warnings;
|
|
|
5
|
|
|
6 if(@ARGV == 1 and $ARGV[0] eq "-v"){
|
|
|
7 print "Version 1.0\n";
|
|
|
8 exit;
|
|
|
9 }
|
|
|
10
|
|
|
11 my $quiet = 0;
|
|
|
12 if(@ARGV and $ARGV[0] =~ /^-q(?:uiet)?$/){
|
|
|
13 $quiet = 1;
|
|
|
14 shift @ARGV;
|
|
|
15 }
|
|
|
16
|
|
|
17 @ARGV >= 6 and not @ARGV%2 or die "Usage: $0 [-q(uiet)] <true|false> <output file> <method_name1> <hgvs_file_1.txt> <method_name2> <hgvs_file_2.txt> ...\n",
|
|
|
18 "Where true or false determines if the genotypes calls should all be reported (if false, they are collapsed to the unique set of values observed)\n";
|
|
|
19
|
|
|
20 my %zygosity;
|
|
|
21 my %quality;
|
|
|
22 my %var_count;
|
|
|
23 my %tot_count;
|
|
|
24 my %methods;
|
|
|
25 my %data;
|
|
|
26 my $chr_prefix_exists = 0; # do any of the tools report chr1, if so we need to ensure 1 get reported as chr1 later for consistency across tools
|
|
|
27
|
|
|
28 my $collapse = $ARGV[0] =~ /^true$/i;
|
|
|
29
|
|
|
30 open(OUT, ">$ARGV[1]")
|
|
|
31 or die "Cannot open $ARGV[1] for writing: $!\n";
|
|
|
32
|
|
|
33 my @headers;
|
|
|
34 my ($chr_column, $zygosity_column, $pvalue_column, $alt_cnt_column, $tot_cnt_column, $transcript_column, $cdna_hgvs_column, $pos_column, $to_column, $srcs_column);
|
|
|
35 for(my $i=2; $i<$#ARGV; $i+=2){
|
|
|
36 my $method_name = $ARGV[$i];
|
|
|
37 my $infile = $ARGV[$i+1];
|
|
|
38 print STDERR "Processing $infile...\n" unless $quiet;
|
|
|
39 open(IN, $infile)
|
|
|
40 or die "Cannot open $infile for reading: $!\n";
|
|
|
41 my $header = <IN>; # header
|
|
|
42 chomp $header;
|
|
|
43 if($i == 2){
|
|
|
44 @headers = split /\t/, $header;
|
|
|
45 for(my $j = 0; $j <= $#headers; $j++){
|
|
|
46 if($headers[$j] eq "Chr"){
|
|
|
47 $chr_column = $j;
|
|
|
48 }
|
|
|
49 elsif($headers[$j] eq "Zygosity"){
|
|
|
50 $zygosity_column = $j;
|
|
|
51 }
|
|
|
52 elsif($headers[$j] eq "P-value"){
|
|
|
53 $pvalue_column = $j;
|
|
|
54 }
|
|
|
55 elsif($headers[$j] eq "Variant Reads"){
|
|
|
56 $alt_cnt_column = $j;
|
|
|
57 }
|
|
|
58 elsif($headers[$j] eq "Total Reads"){
|
|
|
59 $tot_cnt_column = $j;
|
|
|
60 }
|
|
|
61 elsif($headers[$j] eq "Selected transcript"){
|
|
|
62 $transcript_column = $j;
|
|
|
63 }
|
|
|
64 elsif($headers[$j] eq "Transcript HGVS"){
|
|
|
65 $cdna_hgvs_column = $j;
|
|
|
66 }
|
|
|
67 elsif($headers[$j] eq "DNA From"){
|
|
|
68 $pos_column = $j;
|
|
|
69 }
|
|
|
70 elsif($headers[$j] eq "DNA To"){
|
|
|
71 $to_column = $j;
|
|
|
72 }
|
|
|
73 elsif($headers[$j] eq "Sources"){
|
|
|
74 $srcs_column = $j;
|
|
|
75 }
|
|
|
76 }
|
|
|
77 if(defined $srcs_column){
|
|
|
78 # print header as-is
|
|
|
79 print OUT "$header\n";
|
|
|
80 }
|
|
|
81 else{
|
|
|
82 # print header from the first file, with extra column for methods at the end
|
|
|
83 push @headers, "Sources";
|
|
|
84 $srcs_column = $#headers;
|
|
|
85 print OUT join("\t", @headers), "\n";
|
|
|
86 }
|
|
|
87 }
|
|
|
88 else{
|
|
|
89 # Make sure the columns are the same in the various files
|
|
|
90 my @latest_headers = split /\t/, $header;
|
|
|
91 for(my $j = 0; $j <= $#latest_headers && $j <= $#headers; $j++){
|
|
|
92 if($latest_headers[$j] ne $headers[$j]){
|
|
|
93 die "Header column ", $j+1, "($latest_headers[$j]) in $ARGV[$i] is different from the equivalent column label ($headers[$j]) in $ARGV[2]. Aborting, cannot merge the files.\n";
|
|
|
94 }
|
|
|
95 }
|
|
|
96 }
|
|
|
97 while(<IN>){
|
|
|
98 chomp;
|
|
|
99 my @F = split /\t/, $_, -1; # -1 to keep trailing blank fields
|
|
|
100 if(not $chr_prefix_exists and $F[$chr_column] =~ /^chr/){
|
|
|
101 $chr_prefix_exists = 1;
|
|
|
102 }
|
|
|
103 $F[$chr_column] =~ s/^chr//; # for consistency
|
|
|
104 my $key = "$F[$transcript_column]:$F[$cdna_hgvs_column]"; # transcript_id:cdna_hgvs is shared id for common variants amongst files
|
|
|
105 # record disagreement on zygosity if any
|
|
|
106 if(not defined $zygosity{$key}){
|
|
|
107 $zygosity{$key} = [];
|
|
|
108 $quality{$key} = [];
|
|
|
109 $var_count{$key} = [];
|
|
|
110 $tot_count{$key} = [];
|
|
|
111 $data{$key} = \@F;
|
|
|
112 #print STDERR "Missing fields (only have ", scalar(@F), " for input '$_'\n" if $#F < $#headers;
|
|
|
113 }
|
|
|
114 push @{$zygosity{$key}}, split(/; /,$F[$zygosity_column]);
|
|
|
115 push @{$quality{$key}}, $F[$pvalue_column];
|
|
|
116 push @{$var_count{$key}}, $F[$alt_cnt_column];
|
|
|
117 push @{$tot_count{$key}}, $F[$tot_cnt_column];
|
|
|
118 push @{$methods{$key}}, $method_name;
|
|
|
119 }
|
|
|
120 close(IN);
|
|
|
121 }
|
|
|
122
|
|
|
123 if($chr_prefix_exists){
|
|
|
124 for my $key (keys %data){
|
|
|
125 # assuming there is a mix of chr1 and 1, always add the chr for consistency
|
|
|
126 $data{$key}->[$chr_column] = "chr".$data{$key}->[$chr_column] unless $data{$key}->[$chr_column] =~ /^chr/;
|
|
|
127 }
|
|
|
128 }
|
|
|
129
|
|
|
130 # Sort by chr, then pos, then transcript_id
|
|
|
131 for my $key (sort {$data{$a}->[$chr_column] cmp $data{$b}->[$chr_column] or
|
|
|
132 $data{$a}->[$pos_column] <=> $data{$b}->[$pos_column] or
|
|
|
133 $data{$a}->[$to_column] <=> $data{$b}->[$to_column] or
|
|
|
134 $data{$a}->[$transcript_column] cmp $data{$b}->[$transcript_column]} keys %data){
|
|
|
135 # Before printing a line, merge the zygosity, variant quality, read count data if requested
|
|
|
136 my %seen;
|
|
|
137 if($collapse){
|
|
|
138 my @zygosities = grep {$_ ne "NA" and not $seen{$_}++} @{$zygosity{$key}};
|
|
|
139 if(@zygosities == 0){
|
|
|
140 # do nothing, existing NA in 7th column will be a fine answer as they are all like tha
|
|
|
141 }
|
|
|
142 elsif(@zygosities == 1){
|
|
|
143 # agreement by all methods, just print the one not NA
|
|
|
144 $data{$key}->[$zygosity_column] = $zygosities[0] if $data{$key}->[$zygosity_column] eq "NA";
|
|
|
145 }
|
|
|
146 else{
|
|
|
147 $data{$key}->[$zygosity_column] = join("/", sort keys %seen); # list all unique values that occur
|
|
|
148 }
|
|
|
149 }
|
|
|
150 else{
|
|
|
151 $data{$key}->[$zygosity_column] = join("; ", @{$zygosity{$key}});
|
|
|
152 }
|
|
|
153
|
|
|
154 if($collapse){
|
|
|
155 undef %seen;
|
|
|
156 my @qualities = grep {$_ ne "NA" and not $seen{$_}++} @{$quality{$key}};
|
|
|
157 if(@qualities == 0){
|
|
|
158 # do nothing, existing NA in 8th column will be a fine answer as they are all like that
|
|
|
159 }
|
|
|
160 elsif(@qualities == 1){
|
|
|
161 # agreement by all methods, just print the one not NA
|
|
|
162 $data{$key}->[$pvalue_column] = $qualities[0] if $data{$key}->[8] eq "NA";
|
|
|
163 }
|
|
|
164 else{
|
|
|
165 # calculate the median for the collapsed value
|
|
|
166 my @sorted_qualities = sort {$a <=> $b} grep({$_ ne "NA"} @{$quality{$key}});
|
|
|
167 my $median_quality = $sorted_qualities[int($#sorted_qualities/2)]; # rounds down (actually better score as this is a p-value)
|
|
|
168 $data{$key}->[$pvalue_column] = $median_quality;
|
|
|
169 }
|
|
|
170 }
|
|
|
171 else{
|
|
|
172 $data{$key}->[$pvalue_column] = join("; ", @{$quality{$key}});
|
|
|
173 }
|
|
|
174
|
|
|
175 if($collapse){
|
|
|
176 undef %seen;
|
|
|
177 my @var_counts = grep {$_ ne "NA" and not $seen{$_}++} @{$var_count{$key}};
|
|
|
178 undef %seen;
|
|
|
179 my @tot_counts = grep {$_ ne "NA" and not $seen{$_}++} @{$tot_count{$key}};
|
|
|
180 if(@var_counts == 0 and @tot_counts == 0){
|
|
|
181 # do nothing, existing NAs okay
|
|
|
182 }
|
|
|
183 elsif(@var_counts == 1 and @tot_counts == 1){
|
|
|
184 # agreement by all methods, just print the one in %data unless it's NA
|
|
|
185 $data{$key}->[$alt_cnt_column] = $var_counts[0] if $data{$key}->[$alt_cnt_column] eq "NA";
|
|
|
186 $data{$key}->[$tot_cnt_column] = $tot_counts[0] if $data{$key}->[$tot_cnt_column] eq "NA";
|
|
|
187 }
|
|
|
188 else{
|
|
|
189 # calculate the median for the collapsed value
|
|
|
190 my @sorted_var_counts = sort {$a <=> $b} grep({$_ ne "NA"} @{$var_count{$key}});
|
|
|
191 my $pivot = $#sorted_var_counts/2;
|
|
|
192 my $median_var_count = $pivot == int($pivot) ? $sorted_var_counts[$pivot] : # arithmetic mean
|
|
|
193 int(($sorted_var_counts[int($pivot)]+$sorted_var_counts[int($pivot)+1])/2);
|
|
|
194 $data{$key}->[$alt_cnt_column] = $median_var_count;
|
|
|
195 my @sorted_tot_counts = sort {$a <=> $b} grep({$_ ne "NA"} @{$tot_count{$key}});
|
|
|
196 $pivot = $#sorted_tot_counts/2;
|
|
|
197 my $median_tot_count = $pivot == int($pivot) ? $sorted_tot_counts[$pivot] : # arithmetic mean
|
|
|
198 int(($sorted_tot_counts[int($pivot)]+$sorted_tot_counts[int($pivot)+1])/2);
|
|
|
199 $data{$key}->[$tot_cnt_column] = $median_tot_count;
|
|
|
200 }
|
|
|
201 }
|
|
|
202 else{
|
|
|
203 # list the raw numbers
|
|
|
204 $data{$key}->[$alt_cnt_column] = join("; ", @{$var_count{$key}});
|
|
|
205 $data{$key}->[$tot_cnt_column] = join("; ", @{$tot_count{$key}});
|
|
|
206 }
|
|
|
207
|
|
|
208 # to facilitate multiple rounds of combining, slash the extra column from the last round if it exists
|
|
|
209 $data{$key}->[$srcs_column] = join("; ", @{$methods{$key}});
|
|
|
210 for(my $i = 0; $i <= $#{$data{$key}}; $i++){
|
|
|
211 $data{$key}->[$i] = "" if not defined $data{$key}->[$i];
|
|
|
212 }
|
|
|
213 print OUT join("\t", @{$data{$key}}), "\n";
|
|
|
214 }
|
|
|
215 close(OUT);
|
