0
+ − 1 #!/usr/bin/env perl
+ − 2
+ − 3 use strict;
+ − 4 use warnings;
+ − 5 use File::Basename;
+ − 6 # Assume the microarray file looks something like this...
+ − 7 ## comment lines...
+ − 8 #Probe Set ID 11-01451_(GenomeWideSNP_5).brlmm-p.chp Forward Strand Base Calls dbSNP RS ID Chromosome Chromosomal Position
+ − 9 #SNP_A-1780520 GG rs16994928 20 48440771
+ − 10 #SNP_A-1780985 AG rs3859360 18 34178190
+ − 11 #...end example input file
+ − 12 if(@ARGV == 1 and $ARGV[0] eq "-v"){
+ − 13 print "Version 1.0\n";
+ − 14 exit;
+ − 15 }
+ − 16
+ − 17 # configuration file stuff
+ − 18 my $dirname = dirname(__FILE__);
+ − 19 my %config;
+ − 20 my $tool_data = shift @ARGV;
+ − 21 if(not -e "$tool_data/microarray_report.loc"){
+ − 22 system("cp $dirname/tool-data/microarray_report.loc $tool_data/microarray_report.loc");
+ − 23 }
+ − 24 open CONFIG, '<', "$tool_data/microarray_report.loc";
+ − 25 while(<CONFIG>){
+ − 26 (my $key, my $value) = split(/\s+/, $_);
+ − 27 $config{$key} = $value;
+ − 28 }
+ − 29 my $dbs_dir = $config{"dbs_directory"};
+ − 30 close CONFIG;
+ − 31
+ − 32 my $quiet = 0;
+ − 33 if(@ARGV and $ARGV[0] =~ /^-q/){
+ − 34 $quiet = 1;
+ − 35 shift @ARGV;
+ − 36 }
+ − 37
+ − 38 @ARGV == 8 or @ARGV == 9 or die "Usage: $0 [-q(uiet)] <output discordant.bed> <output summary.txt> <input ngs hgvs.txt> <input microarray calls.txt> <input.bam> <reporting cds regions.gtf> <reference genome.fasta> <coverage cutoff for stats> [exome target ngs regions.bed]\n";
+ − 39
+ − 40 print STDERR "Reading in reference sequence...\n" unless $quiet;
+ − 41 my %seq;
+ − 42 open(FASTA, "$dbs_dir/$ARGV[6]" )
+ − 43 or die "Cannot open $dbs_dir/$ARGV[6] for reading: $!\n";
+ − 44 $/ = "\n>";
+ − 45 while(<FASTA>){
+ − 46 chomp;
+ − 47 my ($name) = /^>?(\S+)/;
+ − 48 s/^[^\n]+//;
+ − 49 tr/\r\n//d;
+ − 50 $seq{$name} = $_;
+ − 51 }
+ − 52 close(FASTA);
+ − 53 $/ = "\n";
+ − 54
+ − 55 my $cov_cutoff = $ARGV[7];
+ − 56 my %reporting;
+ − 57 print STDERR "Reading in reporting regions list...\n" unless $quiet;
+ − 58 open(GTF, $ARGV[5])
+ − 59 or die "Cannot open $ARGV[5] for reading: $!\n";
+ − 60 while(<GTF>){
+ − 61 next if /^\s*#/;
+ − 62 my @fields = split /\t/, $_;
+ − 63
+ − 64 if($fields[2] eq "exon"){
+ − 65 if(not exists $reporting{$fields[0]}){
+ − 66 $reporting{$fields[0]} = [[$fields[3], $fields[4]]];
+ − 67 next;
+ − 68 }
+ − 69 push @{$reporting{$fields[0]}}, [$fields[3], $fields[4]];
+ − 70 }
+ − 71 }
+ − 72 close(GTF);
+ − 73 for my $c (keys %reporting){
+ − 74 $reporting{$c} = [sort {$a->[0] <=> $b->[0]} @{$reporting{$c}}];
+ − 75 }
+ − 76
+ − 77 my %regions;
+ − 78 if(@ARGV == 9){
+ − 79 print STDERR "Reading in target regions list...\n" unless $quiet;
+ − 80 open(BED, $ARGV[8])
+ − 81 or die "Cannot open $ARGV[8] for reading: $!\n";
+ − 82 while(<BED>){
+ − 83 chomp;
+ − 84 my @F = split /\t/, $_;
+ − 85 next unless @F > 4;
+ − 86 if(not exists $regions{$F[0]}){
+ − 87 $regions{$F[0]} = [];
+ − 88 }
+ − 89 push @{$regions{$F[0]}}, [$F[1],$F[2]]; # assume they are in order
+ − 90 }
+ − 91 }
+ − 92
+ − 93 open(BED, ">$ARGV[0]")
+ − 94 or die "Cannot open $ARGV[0] for writing: $!\n";
+ − 95 print BED "track name=\"NGSvsMicroarrayDiscordance\" columns=\"chr pos pos microGenotype/ngsGenotype NGSReadDepth\" comment=\"asterisk after genotype indicates likely microarray false positive based on NGS coverage, lack of caveats and no mismatched NGS data\" useScore=\"colour gradient\"\n";
+ − 96
+ − 97 open(MICRO, $ARGV[3])
+ − 98 or die "Cannot open microarray calls $ARGV[3] for reading: $!\n";
+ − 99 my $micro_count = 1;
+ − 100 do {$_ = <MICRO>} while /^#/; #header lines
+ − 101 $micro_count++ while <MICRO>; # get a line count for progress meter later
+ − 102 close(MICRO);
+ − 103 $micro_count = int($micro_count/100);
+ − 104
+ − 105 print STDERR "Reading in NGS calls...\n" unless $quiet;
+ − 106 open(MICRO, $ARGV[3])
+ − 107 or die "Cannot open microarray calls $ARGV[3] for reading: $!\n";
+ − 108 open(HGVS, $ARGV[2])
+ − 109 or die "Cannot open HGVS genotype calls $ARGV[2] for reading: $!\n";
+ − 110 open(SUMMARY, ">$ARGV[1]")
+ − 111 or die "Cannot open $ARGV[1] for writing: $!\n";
+ − 112 my $bam_file = $ARGV[4];
+ − 113 die "Input BAM file does not exist\n" if not -e $bam_file;
+ − 114 die "Input BAM file is not readable\n" if not -r $bam_file;
+ − 115 die "Input BAM file is empty\n" if -z $bam_file;
+ − 116 <HGVS>; # header
+ − 117 my (%ngs_call, %ngs_caveats, %ngs_methods, %ngs_varreads, %ngs_depth, %key2pos, %ignore);
+ − 118 while(<HGVS>){
+ − 119 chomp;
+ − 120 my @F = split /\t/, $_;
+ − 121 my $pos_key = "$F[4]:$F[5]";
+ − 122 # ignore non-SNPs
+ − 123 if($F[2] =~ /c.[\-*]?(\d+)_(\d+)/){ #multi-base events ignored, as microarray probes are likely to be reporting wrong here anyway
+ − 124 my $modlength = $2-$1;
+ − 125 for my $offset (0..$modlength){
+ − 126 if($F[3] eq "-"){
+ − 127 $ignore{"$F[4]:".($F[5]-$offset)} = 1;
+ − 128 }
+ − 129 else{
+ − 130 $ignore{"$F[4]:".($F[5]+$offset)} = 1;
+ − 131 }
+ − 132 }
+ − 133 }
+ − 134 elsif($F[2] =~ /(?:del|ins|inv)/){
+ − 135 next;
+ − 136 }
+ − 137 next unless $F[2] =~ />([ACGT])$/; # only looking for SNPs, to do add MNPs
+ − 138 my $new_base = $1;
+ − 139 $new_base =~ tr/ACGT/TGCA/ if $F[3] eq "-"; # rev comp cDNA HGVS
+ − 140 my $rs_key = $F[11]; # use both position and rsID as patches to hg19 have shifted some positions
+ − 141 my $ngs_call = ($F[6] =~ /homozygote/ ? $new_base : $F[10]) . $new_base;
+ − 142 $ngs_call{$pos_key} = $ngs_call{$rs_key} = $ngs_call;
+ − 143 $ngs_caveats{$pos_key} = $ngs_caveats{$rs_key} = $F[16] if $F[16] =~ /\S/;
+ − 144 $ngs_methods{$pos_key} = $ngs_methods{$rs_key} = $F[18];
+ − 145 $ngs_varreads{$pos_key} = $ngs_varreads{$rs_key} = $F[8];
+ − 146 $ngs_depth{$pos_key} = $ngs_depth{$rs_key} = $F[9];
+ − 147 $key2pos{$rs_key} = $key2pos{$pos_key} = [$F[4],$F[5]];
+ − 148 }
+ − 149 close(HGVS);
+ − 150
+ − 151 print STDERR "Comparing to microarray calls...\n" unless $quiet;
+ − 152 print STDERR " 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%\n" unless $quiet;
+ − 153 my $num_ignored_snps = 0;
+ − 154 my $num_usable_snps = 0;
+ − 155 my $num_targeted_snps = 0;
+ − 156 my %discordant; # key => [zygosity, method, caveats, variant depth, total depth];
+ − 157 my %cov_missing; # List all sites with unknown coverage
+ − 158 do {$_ = <MICRO>} while /^#/; #header lines
+ − 159 while(<MICRO>){
+ − 160 if(not $quiet){
+ − 161 if($.%($micro_count*10) == 1){
+ − 162 print STDERR "|";
+ − 163 }
+ − 164 elsif($.%($micro_count*5) == 1){
+ − 165 print STDERR ":";
+ − 166 }
+ − 167 elsif($.%$micro_count == 1){
+ − 168 print STDERR ".";
+ − 169 }
+ − 170 }
+ − 171 chomp;
+ − 172 my @F = split /\t/, $_;
+ − 173 my $micro_call = $F[1];
+ − 174 next if $micro_call eq "---";
+ − 175 my $rsid = $F[2];
+ − 176 my $chr = $F[3];
+ − 177 next if $chr eq "---";
+ − 178 my $pos = $F[4];
+ − 179 if(exists $ignore{"chr$chr:$pos"}){ # multinucleotide events probably wrong on the microarray
+ − 180 $num_ignored_snps++;
+ − 181 next;
+ − 182 }
+ − 183 $num_usable_snps++;
+ − 184
+ − 185 if(exists $reporting{"chr".$chr}){ # in a region reported in the annotation?
+ − 186 my $in_region = 0;
+ − 187 my $arrayref = $reporting{"chr".$chr};
+ − 188 for(my $search_index = find_earliest_index($pos, $arrayref);$search_index <= $#{$arrayref}; $search_index++){
+ − 189 my $interval = $arrayref->[$search_index];
+ − 190 if($pos >= $interval->[0] and $pos <= $interval->[1]){
+ − 191 $in_region = 1; last;
+ − 192 }
+ − 193 last if $pos < $interval->[0];
+ − 194 }
+ − 195 next unless $in_region;
+ − 196 }
+ − 197
+ − 198 if(exists $regions{"chr".$chr}){ # in the exome target region?
+ − 199 my $in_region = 0;
+ − 200 my $arrayref = $regions{"chr".$chr};
+ − 201 for (my $search_index = find_earliest_index($pos, $arrayref);$search_index <= $#{$arrayref}; $search_index++){
+ − 202 my $interval = $arrayref->[$search_index];
+ − 203 if($pos >= $interval->[0] and $pos <= $interval->[1]){
+ − 204 $in_region = 1; last;
+ − 205 }
+ − 206 last if $pos < $interval->[0];
+ − 207 }
+ − 208
+ − 209 next unless $in_region;
+ − 210 }
+ − 211
+ − 212 $num_targeted_snps++;
+ − 213 my $m = substr($micro_call, 0, 1);
+ − 214 my $m2 = substr($micro_call, 1, 1);
+ − 215 my $key;
+ − 216 if(exists $ngs_call{$rsid}){
+ − 217 $key = $rsid;
+ − 218 }
+ − 219 elsif(exists $ngs_call{"chr$chr:$pos"}){
+ − 220 $key = "chr$chr:$pos";
+ − 221 }
+ − 222 else{
+ − 223 if(not exists $seq{"chr$chr"}){
+ − 224 warn "Skipping microarray call, no reference sequence was provided for chr$chr\n";
+ − 225 $num_targeted_snps--;
+ − 226 next;
+ − 227 }
+ − 228 # 4 situations could be e.g. ref AA, but micro says AA or AG or GG or GC
+ − 229 my $ref_base = uc(substr($seq{"chr$chr"}, $pos-1, 1));
+ − 230 if($m eq $m2){
+ − 231 if($micro_call eq $ref_base.$ref_base){
+ − 232 # concordant homo calls as reference
+ − 233 # undefs will be filled in later en masse in a single call to samtools for efficiency
+ − 234 $discordant{$rsid} = ["micro ref, ngs ref", "", "", undef, undef, "chr$chr", $pos, "$micro_call/$micro_call"];
+ − 235 }
+ − 236 else{
+ − 237 # called homo var by microarray, but homo ref by NGS
+ − 238 #print STDERR "chr$chr $pos micro $micro_call vs NGS homo ref $ref_base";
+ − 239 $discordant{$rsid} = ["micro homo, ngs ref", "", "", undef, undef, "chr$chr", $pos, "$micro_call/$ref_base$ref_base"];
+ − 240 }
+ − 241 }
+ − 242 else { #$m ne $m2
+ − 243 if($m eq $ref_base or $m2 eq $ref_base){
+ − 244 # called het var by microarray, but homo ref by NGS
+ − 245 $discordant{$rsid} = ["micro het, ngs ref", "", "", undef, undef, "chr$chr", $pos, "$micro_call/$ref_base$ref_base"];
+ − 246 }
+ − 247 else{
+ − 248 $discordant{$rsid} = ["micro compound het, ngs ref", "", "", undef, undef, "chr$chr", $pos, "$micro_call/$ref_base$ref_base"];
+ − 249 }
+ − 250 }
+ − 251 $cov_missing{"chr$chr:$pos"} = [$rsid, $m, $m2];
+ − 252 next;
+ − 253 }
+ − 254
+ − 255 next if $ngs_depth{$key} <= $cov_cutoff;
+ − 256 my $ngs_call = $ngs_call{$key};
+ − 257 # if we get to here, there are both NGS and micro calls
+ − 258 if($micro_call eq $ngs_call or $micro_call eq reverse($ngs_call)){
+ − 259 #print STDERR "Concordant NGS call for chr$chr:$pos\n";
+ − 260 if($m eq $m2){
+ − 261 $discordant{$key} = ["micro homo, ngs homo", $ngs_methods{$key}, $ngs_caveats{$key}, $ngs_varreads{$key}, $ngs_depth{$key}, "chr$chr", $pos, "$micro_call/$ngs_call"];
+ − 262 }
+ − 263 else{
+ − 264 $discordant{$key} = ["micro het, ngs het", $ngs_methods{$key}, $ngs_caveats{$key}, $ngs_varreads{$key}, $ngs_depth{$key}, "chr$chr", $pos, "$micro_call/$ngs_call"];
+ − 265 }
+ − 266 next;
+ − 267 }
+ − 268 else{
+ − 269 # discordant
+ − 270 #print STDERR "Discordant NGS call for chr$chr:$pos $micro_call vs. $ngs_call\n";
+ − 271 # is it just a zygosity difference?
+ − 272 my $n = substr($ngs_call, 0, 1);
+ − 273 my $n2 = substr($ngs_call, 1, 1);
+ − 274 my $discordance;
+ − 275 if($micro_call eq $n.$n or $micro_call eq $n2.$n2){ # micro is homo for variant, ngs is het
+ − 276 $discordance = "micro homo, ngs het";
+ − 277 }
+ − 278 elsif($m.$m eq $ngs_call or $m2.$m2 eq $ngs_call){ # micro is het for variant, ngs is homo
+ − 279 $discordance = "micro het, ngs homo";
+ − 280 }
+ − 281 elsif($m eq $m2){
+ − 282 if($n eq $n2){
+ − 283 $discordance = "diff micro homo, ngs homo";
+ − 284 }
+ − 285 else{
+ − 286 $discordance = "diff micro homo, ngs het";
+ − 287 }
+ − 288 }
+ − 289 else{
+ − 290 if($n eq $n2){
+ − 291 $discordance = "diff micro het, ngs homo";
+ − 292 }
+ − 293 else{
+ − 294 $discordance = "diff micro het, ngs het";
+ − 295 }
+ − 296 }
+ − 297 $discordant{$key} = [$discordance, $ngs_methods{$key}, $ngs_caveats{$key}, $ngs_varreads{$key}, $ngs_depth{$key}, "chr$chr", $pos, "$micro_call/$ngs_call"];
+ − 298 }
+ − 299 }
+ − 300 close(MICRO);
+ − 301 print STDERR "\n" unless $quiet;
+ − 302
+ − 303 print STDERR "Retrieving reference call depth of coverage stats...\n" unless $quiet;
+ − 304 my $covbed = "$$.bed";
+ − 305 open(TMPCOV, ">$covbed")
+ − 306 or die "Cannot open temporary BED file for samtools: $!\n";
+ − 307 my $num_covmissing = 0;
+ − 308 for(sort keys %cov_missing){
+ − 309 $num_covmissing++;
+ − 310 my ($chr, $pos) = split /:/, $_;
+ − 311 print TMPCOV "$chr\t$pos\n";
+ − 312 }
+ − 313 close(TMPCOV);
+ − 314
+ − 315 my $cov_count = int($num_covmissing/100);
+ − 316 print STDERR " 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%\n" unless $quiet;
+ − 317 #print STDERR "samtools mpileup -l $covbed -I $bam_file 2>/dev/null\n";
+ − 318 #<STDIN>;
+ − 319 open(SAM, "samtools mpileup -l $covbed -I $bam_file 2>/dev/null |")
+ − 320 or die "Cannot run samtools: $!\n";;
+ − 321 while(<SAM>){
+ − 322 #while($depth_data =~ /([^\n]+)/sg){
+ − 323 if(not $quiet){
+ − 324 if($.%($cov_count*10) == 0){
+ − 325 print STDERR "|";
+ − 326 }
+ − 327 elsif($.%($cov_count*5) == 0){
+ − 328 print STDERR ":";
+ − 329 }
+ − 330 elsif($.%$cov_count == 0){
+ − 331 print STDERR ".";
+ − 332 }
+ − 333 }
+ − 334
+ − 335 chomp;
+ − 336 my @B = split /\t/, $_;
+ − 337 #my @B = split /\t/, $1;
+ − 338 my ($rsid, $m, $m2) = @{$cov_missing{"$B[0]:$B[1]"}};
+ − 339 my $depth = @B > 1 ? $B[3] : 0;
+ − 340 #print STDERR "Depth is $depth for $B[0]:$B[1]\n";
+ − 341 my $read_calls = uc($B[4]);
+ − 342 my %base_tot;
+ − 343 for my $read_call (split //, $read_calls){
+ − 344 $base_tot{$read_call}++;
+ − 345 }
+ − 346 $base_tot{$m} = 0 if not exists $base_tot{$m};
+ − 347 $base_tot{$m2} = 0 if not exists $base_tot{$m2};
+ − 348 if($depth <= $cov_cutoff){
+ − 349 $discordant{$rsid}->[0] = "no ngs call";
+ − 350 $discordant{$rsid}->[1] = "insufficient coverage";
+ − 351 $discordant{$rsid}->[3] = $base_tot{$m};
+ − 352 $discordant{$rsid}->[4] = $depth;
+ − 353 }
+ − 354 elsif($discordant{$rsid}->[0] eq "micro ref, ngs ref" or $discordant{$rsid}->[0] eq "micro homo, ngs ref"){
+ − 355 # concordant homo calls as reference or called homo var by microarray, but homo ref by NGS
+ − 356 $discordant{$rsid}->[3] = $base_tot{$m};
+ − 357 $discordant{$rsid}->[4] = $depth;
+ − 358 }
+ − 359 elsif($discordant{$rsid}->[0] eq "micro het, ngs ref" or $discordant{$rsid}->[0] eq "micro compound het, ngs ref"){
+ − 360 # called het var by microarray, but homo ref by NGS
+ − 361 $discordant{$rsid}->[3] = $base_tot{$m}.",".$base_tot{$m2};
+ − 362 $discordant{$rsid}->[4] = $depth;
+ − 363 }
+ − 364 else{
+ − 365 print STDERR "Didn't know how to deal with $B[0]:$B[1] -> ", join(" ", @{$discordant{$rsid}}),"\n" unless $quiet;
+ − 366 }
+ − 367 }
+ − 368 unlink $covbed;
+ − 369
+ − 370 #Count false negatives
+ − 371 my $false_neg_homo_count = 0;
+ − 372 my $false_neg_het_count = 0;
+ − 373 my $missing_count = 0;
+ − 374 my $wrong_base_count = 0;
+ − 375 my $wrong_base_caveat_count = 0;
+ − 376 my $discordant_caveat_count = 0;
+ − 377 my $tot_caveat_count = 0;
+ − 378 my $false_homo_count = 0;
+ − 379 my $false_het_count = 0;
+ − 380 my $micro_fdr_estimate = 0;
+ − 381 my %calls_by_method;
+ − 382 my %concordance_by_method;
+ − 383 for my $key (keys %discordant){
+ − 384 my $rec = $discordant{$key};
+ − 385 my $name = $rec->[7];
+ − 386 $tot_caveat_count++ if $rec->[2];
+ − 387 if($rec->[0] eq "micro het, ngs ref" or $rec->[0] eq "micro homo, ngs ref"){
+ − 388 if($rec->[0] eq "micro homo, ngs ref"){
+ − 389 $false_neg_homo_count++;
+ − 390 }
+ − 391 else{
+ − 392 $false_neg_het_count++;
+ − 393 }
+ − 394 # probably false micro call if no caveats, high coverage, and no/one mismatches in NGS
+ − 395 next unless defined $rec and defined $rec->[3] and defined $rec->[4];
+ − 396 if(not $rec->[2] and $rec->[4] >= 50 and $rec->[3] =~ /(\d{2,})/ and $1+1 >= $rec->[4]){
+ − 397 $micro_fdr_estimate++;
+ − 398 $name .= "*";
+ − 399 }
+ − 400 }
+ − 401 elsif($rec->[0] eq "micro compound het, ngs ref"){
+ − 402 $false_neg_het_count+=2;
+ − 403 # probably false micro call if no caveats, high coverage, and no/one mismatches in NGS
+ − 404 next unless defined $rec and defined $rec->[3] and defined $rec->[4];
+ − 405 if(not $rec->[2] and $rec->[4] >= 50 and $rec->[3] =~ /(\d{2,})/ and $1+1 >= $rec->[4]){
+ − 406 $micro_fdr_estimate+=2;
+ − 407 $name .= "*";
+ − 408 }
+ − 409 }
+ − 410 elsif($rec->[0] eq "no ngs call"){
+ − 411 $missing_count++;
+ − 412 }
+ − 413 elsif($rec->[0] =~ /^diff/){
+ − 414 $wrong_base_count++;
+ − 415 if($rec->[3] =~ /\S/){
+ − 416 $wrong_base_caveat_count++;
+ − 417 }
+ − 418 }
+ − 419 elsif($rec->[0] eq "micro het, ngs homo"){
+ − 420 $false_homo_count++;
+ − 421 }
+ − 422 elsif($rec->[0] eq "micro homo, ngs het"){
+ − 423 $false_het_count++;
+ − 424 }
+ − 425 else{
+ − 426 # calls concordant
+ − 427 $concordance_by_method{$rec->[1]}++;
+ − 428 next;
+ − 429 }
+ − 430 if($rec->[4] > $cov_cutoff and $rec->[2] and $name !~ /\*$/){
+ − 431 $discordant_caveat_count++;
+ − 432 }
+ − 433 my $chr = $rec->[5];
+ − 434 my $pos = $rec->[6];
+ − 435 my $score = $rec->[4];
+ − 436 $score = 1000 if $score > 1000;
+ − 437 print BED "$chr\t$pos\t$pos\t$name\t$score\n";
+ − 438 }
+ − 439
+ − 440 print SUMMARY "# NGS genotypes in $ARGV[2] vs. SNP microarray in $ARGV[3], minimum NGS coverage of ",$cov_cutoff+1, "\n";
+ − 441 print SUMMARY "#Columns: Measure\tCount\tPercentage\n";
+ − 442 print SUMMARY "Total ignored SNP microarray calls due to NGS putative indels or MNPs\t$num_ignored_snps\n";
+ − 443 print SUMMARY "Total usable SNP microarray calls\t$num_usable_snps\n";
+ − 444 printf SUMMARY "Total targeted SNP microarray calls (based on target file %s)\t%d\t%.2f\n", $ARGV[5],$num_targeted_snps,$num_targeted_snps/$num_usable_snps*100 if keys %regions;
+ − 445 printf SUMMARY "Targeted SNPs with insufficient NGS coverage (<=$cov_cutoff)\t%d\t%.2f\n", $missing_count, $missing_count/$num_targeted_snps*100;
+ − 446 $num_targeted_snps -= $missing_count;
+ − 447 my $tot_bad = $wrong_base_count+$false_neg_homo_count+$false_neg_het_count+$false_homo_count+$false_het_count;
+ − 448 printf SUMMARY "Total discordance\t%d\t%.2f\n", $tot_bad, $tot_bad/$num_targeted_snps*100;
+ − 449 $tot_bad -= $discordant_caveat_count+$micro_fdr_estimate;
+ − 450 printf SUMMARY "Significant discordance (excludes NGS calls with caveats, microarray het FDR)\t%d\t%.2f\n", $tot_bad, $tot_bad/$num_targeted_snps*100;
+ − 451 printf SUMMARY "Caveats discordant\t%d\t%.2f\n", $discordant_caveat_count, $tot_caveat_count == 0 ? 0 : $discordant_caveat_count/$tot_caveat_count*100;
+ − 452 printf SUMMARY "Incorrect NGS base called\t%d\t%.2f\n", $wrong_base_count, $wrong_base_count/$num_targeted_snps*100;
+ − 453 printf SUMMARY "Incorrect NGS base called, subset with caveats\t%d\t%.2f\n", $wrong_base_caveat_count, $wrong_base_count == 0 ? 0 : $wrong_base_caveat_count/$wrong_base_count*100;
+ − 454 printf SUMMARY "False negative NGS homo\t%d\t%.2f\n", $false_neg_homo_count, $false_neg_homo_count/$num_targeted_snps*100;
+ − 455 printf SUMMARY "False negative NGS het\t%d\t%.2f\n", $false_neg_het_count, $false_neg_het_count/$num_targeted_snps*100;
+ − 456 printf SUMMARY "Microarray est. FDR het\t%d\t%.2f\n", $micro_fdr_estimate, $micro_fdr_estimate/$num_targeted_snps*100;
+ − 457 printf SUMMARY "Het called NGS homo\t%d\t%.2f\n", $false_homo_count, $false_homo_count/$num_targeted_snps*100;
+ − 458 printf SUMMARY "Homo called NGS het\t%d\t%.2f\n", $false_het_count, $false_het_count/$num_targeted_snps*100;
+ − 459 for(sort keys %concordance_by_method){
+ − 460 printf SUMMARY "%s true positives\t%d\t%.2f\n", (length($_) ? $_ : "Reference"), $concordance_by_method{$_}, $concordance_by_method{$_}/$num_targeted_snps*100;
+ − 461 }
+ − 462 sub find_earliest_index{
+ − 463 # 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
+ − 464 my ($query, $array) = @_;
+ − 465
+ − 466 return 0 if $query < $array->[0]->[0];
+ − 467
+ − 468 my ($left, $right, $prevCenter) = (0, $#$array, -1);
+ − 469
+ − 470 while(1) {
+ − 471 my $center = int (($left + $right)/2);
+ − 472
+ − 473 my $cmp = $query <=> $array->[$center]->[0] || ($center == 0 || $query != $array->[$center-1]->[0] ? 0 : -1);
+ − 474
+ − 475 return $center if $cmp == 0;
+ − 476 if ($center == $prevCenter) {
+ − 477 return $left;
+ − 478 }
+ − 479 $right = $center if $cmp < 0;
+ − 480 $left = $center if $cmp > 0;
+ − 481 $prevCenter = $center;
+ − 482 }
+ − 483 }
