|
0
|
1 #!/usr/bin/env perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use warnings;
|
|
|
5
|
|
|
6
|
|
|
7 my $quiet = 0;
|
|
|
8 if(@ARGV and $ARGV[0] =~ /^-q/){
|
|
|
9 $quiet = 1;
|
|
|
10 shift @ARGV;
|
|
|
11 }
|
|
|
12 # Calculates Cohen's kappa measure for inter-annotator agreement, here for two genotype callers
|
|
|
13 # SPECIAL ATTENTION MUST BE PAID to include or exclude calls using depth of coverage information, and only look at targeted regions (e.g. exome analysis)
|
|
|
14 # Also, the probability of a non-call must be incorporated, as there is going to be 99.9% agreement by default due to the low rate of
|
|
|
15 # SNPs in the genome in human, for example.
|
|
|
16 @ARGV == 9 or die "Usage: $0 [-q(uiet)] <exome targets.bed> <coding regions.gtf> <genotypes1.hgvs.txt> <bam1> <coverage cutoff1> <genotypes2.hgvs.txt> <bam2> <coverage cutoff2> <output report.txt>\n";
|
|
|
17
|
|
|
18 my $targets_file = shift @ARGV;
|
|
|
19 my $coding_file = shift @ARGV;
|
|
|
20 my $calls1_file = shift @ARGV;
|
|
|
21 my $bam_file1 = shift @ARGV;
|
|
|
22 my $calls1_min_depth = shift @ARGV;
|
|
|
23 my $calls2_file = shift @ARGV;
|
|
|
24 my $bam_file2 = shift @ARGV;
|
|
|
25 my $calls2_min_depth = shift @ARGV;
|
|
|
26 my $output_file = shift @ARGV;
|
|
|
27
|
|
|
28 my %reporting;
|
|
|
29 print STDERR "Reading in target regions list...\n" unless $quiet;
|
|
|
30 open(BED, $targets_file)
|
|
|
31 or die "Cannot open target BED file $targets_file for reading: $!\n";
|
|
|
32 while(<BED>){
|
|
|
33 next if /^(\s*#|track=)/;
|
|
|
34 chomp;
|
|
|
35 my @fields = split /\t/, $_;
|
|
|
36
|
|
|
37 if(not exists $reporting{$fields[0]}){
|
|
|
38 $reporting{$fields[0]} = [[$fields[1], $fields[2]]];
|
|
|
39 next;
|
|
|
40 }
|
|
|
41 push @{$reporting{$fields[0]}}, [$fields[1], $fields[2]];
|
|
|
42 }
|
|
|
43 close(BED);
|
|
|
44 for my $c (keys %reporting){
|
|
|
45 $reporting{$c} = [sort {$a->[0] <=> $b->[0]} @{$reporting{$c}}];
|
|
|
46 }
|
|
|
47
|
|
|
48 my %coding;
|
|
|
49 print STDERR "Reading in coding regions list...\n" unless $quiet;
|
|
|
50 open(GTF, $coding_file)
|
|
|
51 or die "Cannot open $coding_file for reading: $!\n";
|
|
|
52 while(<GTF>){
|
|
|
53 next if /^\s*#/;
|
|
|
54 my @fields = split /\t/, $_;
|
|
|
55
|
|
|
56 if($fields[2] eq "CDS"){
|
|
|
57 if(not exists $coding{$fields[0]}){
|
|
|
58 $coding{$fields[0]} = [[$fields[3], $fields[4]]];
|
|
|
59 next;
|
|
|
60 }
|
|
|
61 push @{$coding{$fields[0]}}, [$fields[3], $fields[4]];
|
|
|
62 }
|
|
|
63 }
|
|
|
64 close(GTF);
|
|
|
65 for my $c (keys %coding){
|
|
|
66 $coding{$c} = [sort {$a->[0] <=> $b->[0]} @{$coding{$c}}];
|
|
|
67 }
|
|
|
68
|
|
|
69 print STDERR "Reading in genotype calls...\n" unless $quiet;
|
|
|
70 open(OUT, ">$output_file")
|
|
|
71 or die "Cannot open $output_file for writing: $!\n";
|
|
|
72 my ($calls1, $subpar_calls1) = retrieve_calls($calls1_file, \%reporting, \%coding, $calls1_min_depth);
|
|
|
73 my ($calls2, $subpar_calls2) = retrieve_calls($calls2_file, \%reporting, \%coding, $calls2_min_depth);
|
|
|
74
|
|
|
75 my $total_calls1 = 0;
|
|
|
76 for my $chr (keys %$calls1){
|
|
|
77 $total_calls1 += scalar(keys %{$calls1->{$chr}});
|
|
|
78 }
|
|
|
79 my $total_calls2 = 0;
|
|
|
80 for my $chr (keys %$calls2){
|
|
|
81 $total_calls2 += scalar(keys %{$calls2->{$chr}});
|
|
|
82 }
|
|
|
83
|
|
|
84 print OUT "Total Method 1 calls\t$total_calls1\n";
|
|
|
85 print OUT "Total Method 2 calls\t$total_calls2\n";
|
|
|
86 print STDERR "Comparing genotype calls...\n" unless $quiet;
|
|
|
87 # Go through the BAM file in all targeted locations to see if there is a reference or variant call (HGVS file only shows variants)
|
|
|
88 my $bases_diff = 0;
|
|
|
89 my $bases_same = 0;
|
|
|
90 my $zygosities_diff = 0;
|
|
|
91 my $zygosities_same = 0;
|
|
|
92 for my $chr (keys %$calls1){
|
|
|
93 for my $pos (keys %{$calls1->{$chr}}){
|
|
|
94 if(exists $calls2->{$chr}->{$pos}){
|
|
|
95 # same variant base
|
|
|
96 if($calls2->{$chr}->{$pos}->[0] eq $calls1->{$chr}->{$pos}->[0]){
|
|
|
97 $bases_same++;
|
|
|
98 if($calls1->{$chr}->{$pos}->[1] eq $calls2->{$chr}->{$pos}->[1]){
|
|
|
99 $zygosities_same++;
|
|
|
100 }
|
|
|
101 else{
|
|
|
102 $zygosities_diff++;
|
|
|
103 }
|
|
|
104 }
|
|
|
105 # diff variant base
|
|
|
106 else{
|
|
|
107 $bases_diff++;
|
|
|
108 }
|
|
|
109 delete $calls1->{$chr}->{$pos}; # so in the end all that's left is method-1-only calls in the hash
|
|
|
110 delete $calls2->{$chr}->{$pos}; # so in the end all that's left is method-2-only calls in the hash
|
|
|
111 delete $subpar_calls1->{$chr}->{$pos} if exists $subpar_calls1->{$chr}; # so in the end all that's left is method-1-only calls in the hash
|
|
|
112 delete $subpar_calls2->{$chr}->{$pos} if exists $subpar_calls2->{$chr}; # so in the end all that's left is method-2-only calls in the hash
|
|
|
113 }
|
|
|
114 }
|
|
|
115 }
|
|
|
116
|
|
|
117 print OUT "Same calls\t$zygosities_same\n";
|
|
|
118 print OUT "Diff zygosity\t$zygosities_diff\n";
|
|
|
119 print OUT "Diff bases\t$bases_diff\n";
|
|
|
120
|
|
|
121 my $covbed = "$$.bed";
|
|
|
122 my $num_to_check = 0;
|
|
|
123 open(TMPCOV, ">$covbed")
|
|
|
124 or die "Cannot open temporary BED file $covbed for samtools: $!\n";
|
|
|
125 for my $chr (keys %$calls1){
|
|
|
126 for my $pos (keys %{$calls1->{$chr}}){
|
|
|
127 if(exists $subpar_calls2->{$chr} and exists $subpar_calls2->{$chr}->{$pos}){
|
|
|
128 next; #exclude counting positions with sufficent coverage in one analysis, but not the other
|
|
|
129 }
|
|
|
130 else{
|
|
|
131 $num_to_check++;
|
|
|
132 print TMPCOV "$chr\t$pos\n";
|
|
|
133 }
|
|
|
134 }
|
|
|
135 }
|
|
|
136 close(TMPCOV);
|
|
|
137
|
|
|
138 print STDERR "Checking $num_to_check read depths for missing calls in $bam_file2...\n" unless $quiet;
|
|
|
139 my $remaining_calls1_homo = 0;
|
|
|
140 my $remaining_calls1_het = 0;
|
|
|
141 my $remaining_calls1_novel = 0;
|
|
|
142 my $remaining_calls1_novel_homo = 0;
|
|
|
143 open(SAM, "samtools mpileup -l $covbed -I $bam_file2 2>/dev/null |")
|
|
|
144 or die "Cannot run samtools: $!\n";
|
|
|
145 while(<SAM>){
|
|
|
146 my @B = split /\t/, $_;
|
|
|
147 my $depth = @B > 1 ? $B[3] : 0;
|
|
|
148
|
|
|
149 # The two methods may have used different patch versions of the same genome. Check if there was a no-call due to different reference bases.
|
|
|
150 my $called_base = @B > 1 ? $B[2] : "N";
|
|
|
151 if($depth and $called_base eq $calls1->{$B[0]}->{$B[1]}->[0]){
|
|
|
152 next;
|
|
|
153 }
|
|
|
154
|
|
|
155 if($depth >= $calls2_min_depth){ # have good coverage in sample 2, but didn't make a SNP call like in sample 1
|
|
|
156 if($calls1->{$B[0]}->{$B[1]}->[1] =~ /het/){
|
|
|
157 $remaining_calls1_het++;
|
|
|
158 }
|
|
|
159 else{
|
|
|
160 $remaining_calls1_homo++;
|
|
|
161 }
|
|
|
162 if($calls1->{$B[0]}->{$B[1]}->[3] eq "NA" or $calls1->{$B[0]}->{$B[1]}->[3] < 0.001){
|
|
|
163 $remaining_calls1_novel++;
|
|
|
164 $remaining_calls1_novel_homo++ if $calls1->{$B[0]}->{$B[1]}->[1] !~ /het/;
|
|
|
165 }
|
|
|
166 else{
|
|
|
167 #print STDERR "Got method 1 only non-novel call: ", join("/", @{$calls1->{$B[0]}->{$B[1]}}), "\n";
|
|
|
168 }
|
|
|
169 }
|
|
|
170 }
|
|
|
171 print OUT "Method 1 only\t",($remaining_calls1_homo+$remaining_calls1_het),"\n";
|
|
|
172 printf OUT "Method 1 only, het %%\t%.1f\n",$remaining_calls1_het/($remaining_calls1_homo+$remaining_calls1_het)*100;
|
|
|
173 printf OUT "Method 1 only, novel %%\t%.1f\n",$remaining_calls1_novel/($remaining_calls1_homo+$remaining_calls1_het)*100;
|
|
|
174 printf OUT "Method 1 only, novel, homo %%\t%.1f\n",$remaining_calls1_novel_homo/$remaining_calls1_novel*100;
|
|
|
175
|
|
|
176 $num_to_check = 0;
|
|
|
177 open(TMPCOV, ">$covbed")
|
|
|
178 or die "Cannot open temporary BED file $covbed for samtools: $!\n";
|
|
|
179 for my $chr (keys %$calls2){
|
|
|
180 for my $pos (keys %{$calls2->{$chr}}){
|
|
|
181 if(exists $subpar_calls1->{$chr} and exists $subpar_calls1->{$chr}->{$pos}){
|
|
|
182 next; #exclude counting positions with sufficent coverage in one analysis, but not the other
|
|
|
183 }
|
|
|
184 else{
|
|
|
185 $num_to_check++;
|
|
|
186 print TMPCOV "$chr\t$pos\n";
|
|
|
187 }
|
|
|
188 }
|
|
|
189 }
|
|
|
190 close(TMPCOV);
|
|
|
191
|
|
|
192 print STDERR "Checking $num_to_check read depths for missing calls in $bam_file1...\n" unless $quiet;
|
|
|
193 my $remaining_calls2_het = 0;
|
|
|
194 my $remaining_calls2_homo = 0;
|
|
|
195 my $remaining_calls2_novel = 0;
|
|
|
196 my $remaining_calls2_novel_homo = 0;
|
|
|
197 open(SAM, "samtools mpileup -l $covbed -I $bam_file1 2>/dev/null |")
|
|
|
198 or die "Cannot run samtools: $!\n";
|
|
|
199 while(<SAM>){
|
|
|
200 my @B = split /\t/, $_;
|
|
|
201 my $depth = @B > 1 ? $B[3] : 0;
|
|
|
202
|
|
|
203 # The two methods may have used different patch versions of the same genome. Check if there was a no-call due to different reference bases.
|
|
|
204 my $called_base = @B > 1 ? $B[2] : "N";
|
|
|
205 if($depth and $called_base eq $calls2->{$B[0]}->{$B[1]}->[0]){
|
|
|
206 next;
|
|
|
207 }
|
|
|
208
|
|
|
209 if($depth >= $calls1_min_depth){ # have good coverage in sample 1, but didn't make a SNP call like in sample 2
|
|
|
210 if($calls2->{$B[0]}->{$B[1]}->[1] =~ /het/){
|
|
|
211 $remaining_calls2_het++;
|
|
|
212 }
|
|
|
213 else{
|
|
|
214 $remaining_calls2_homo++;
|
|
|
215 }
|
|
|
216 if($calls2->{$B[0]}->{$B[1]}->[3] eq "NA" or $calls2->{$B[0]}->{$B[1]}->[3] < 0.001){
|
|
|
217 $remaining_calls2_novel++;
|
|
|
218 $remaining_calls2_novel_homo++ if $calls2->{$B[0]}->{$B[1]}->[1] !~ /het/;
|
|
|
219 }
|
|
|
220 }
|
|
|
221 }
|
|
|
222
|
|
|
223 print OUT "Method 2 only\t",($remaining_calls2_homo+$remaining_calls2_het),"\n";
|
|
|
224 printf OUT "Method 2 only, het %%\t%.1f\n",$remaining_calls2_het/($remaining_calls2_homo+$remaining_calls2_het)*100;
|
|
|
225 printf OUT "Method 2 only, novel %%\t%.1f\n",$remaining_calls2_novel/($remaining_calls2_homo+$remaining_calls2_het)*100;
|
|
|
226 printf OUT "Method 2 only, novel, homo %%\t%.1f\n",$remaining_calls2_novel_homo/$remaining_calls2_novel*100;
|
|
|
227 my $total_calls = $zygosities_same+$zygosities_diff+$bases_diff+$remaining_calls1_homo+$remaining_calls1_het+$remaining_calls2_homo+$remaining_calls2_het;
|
|
|
228 printf OUT "Total coding sequence genotype concordance%%\t%.1f\n", $zygosities_same/$total_calls*100;
|
|
|
229 printf OUT "Total coding sequence zygosity concordance%%\t%.1f\n", 100-$zygosities_diff/$total_calls*100;
|
|
|
230
|
|
|
231 unlink $covbed;
|
|
|
232
|
|
|
233 sub retrieve_calls{
|
|
|
234 my ($hgvs_file, $reporting, $coding, $coverage_cutoff) = @_;
|
|
|
235 my %calls;
|
|
|
236 my %bad_coverage;
|
|
|
237 open(HGVS, $hgvs_file)
|
|
|
238 or die "Cannot open the first genotypes file $hgvs_file for reading: $!\n";
|
|
|
239 <HGVS>; # header
|
|
|
240 while(<HGVS>){
|
|
|
241 my @F = split /\t/, $_;
|
|
|
242 my $HGVS = $F[2];
|
|
|
243 next unless $HGVS =~ />([ACGT])/; # only look at SNPs
|
|
|
244 my $new_base = $1;
|
|
|
245 $new_base =~ tr/ACGT/TGCA/ if $F[3] eq "-"; # rev comp cDNA HGVS
|
|
|
246 my $chr = $F[4];
|
|
|
247 my $pos = $F[5];
|
|
|
248 my $zygosity = $F[6];
|
|
|
249 my $coverage = $F[9];
|
|
|
250 my $maf = $F[14];
|
|
|
251 if($coverage_cutoff > $coverage or $F[16] ne ""){
|
|
|
252 $bad_coverage{$chr} = {} unless exists $bad_coverage{$chr};
|
|
|
253 $bad_coverage{$chr}->{$pos} = [$new_base, $zygosity, $coverage, $maf] unless exists $bad_coverage{$chr}->{$pos};
|
|
|
254 next;
|
|
|
255 }
|
|
|
256
|
|
|
257 if(exists $reporting->{$chr}){ # in a region reported in the annotation, and it's protein coding?
|
|
|
258 my $in_region = 0;
|
|
|
259 my $arrayref = $reporting->{$chr};
|
|
|
260 for(my $search_index = find_earliest_index($pos, $arrayref);$search_index <= $#{$arrayref}; $search_index++){
|
|
|
261 my $interval = $arrayref->[$search_index];
|
|
|
262 if($pos >= $interval->[0] and $pos <= $interval->[1]){
|
|
|
263 $in_region = 1; last;
|
|
|
264 }
|
|
|
265 last if $pos < $interval->[0];
|
|
|
266 }
|
|
|
267 next unless $in_region;
|
|
|
268 $in_region = 0;
|
|
|
269 $arrayref = $coding->{$chr};
|
|
|
270 for(my $search_index = find_earliest_index($pos, $arrayref);$search_index <= $#{$arrayref}; $search_index++){
|
|
|
271 my $interval = $arrayref->[$search_index];
|
|
|
272 if($pos >= $interval->[0] and $pos <= $interval->[1]){
|
|
|
273 $in_region = 1; last;
|
|
|
274 }
|
|
|
275 last if $pos < $interval->[0];
|
|
|
276 }
|
|
|
277 next unless $in_region;
|
|
|
278 }
|
|
|
279
|
|
|
280 if(not exists $calls{$chr}){
|
|
|
281 $calls{$chr} = {};
|
|
|
282 }
|
|
|
283 next if exists $calls{$chr}->{$pos};
|
|
|
284 $calls{$chr}->{$pos} = [$new_base, $zygosity, $coverage, $maf];
|
|
|
285 }
|
|
|
286 close(HGVS);
|
|
|
287 return (\%calls, \%bad_coverage);
|
|
|
288 }
|
|
|
289
|
|
|
290
|
|
|
291 sub find_earliest_index{
|
|
|
292 # employs a binary search to find the smallest index that must be the starting point of a search of [start,end] elements sorted in an array by start
|
|
|
293 my ($query, $array) = @_;
|
|
|
294
|
|
|
295 return 0 if $query < $array->[0]->[0];
|
|
|
296
|
|
|
297 my ($left, $right, $prevCenter) = (0, $#$array, -1);
|
|
|
298
|
|
|
299 while(1) {
|
|
|
300 my $center = int (($left + $right)/2);
|
|
|
301
|
|
|
302 my $cmp = $query <=> $array->[$center]->[0] || ($center == 0 || $query != $array->[$center-1]->[0] ? 0 : -1);
|
|
|
303
|
|
|
304 return $center if $cmp == 0;
|
|
|
305 if ($center == $prevCenter) {
|
|
|
306 return $left;
|
|
|
307 }
|
|
|
308 $right = $center if $cmp < 0;
|
|
|
309 $left = $center if $cmp > 0;
|
|
|
310 $prevCenter = $center;
|
|
|
311 }
|
|
|
312 }
|
|
|
313
|
