Changed calculation of call rate to account for zeroed out SNPs. (#1711)

broadinstitute · Sep 7, 2021 · 1372326 · 1372326
1 parent 20ddad5
commit 1372326
Showing 1 changed file with 30 additions and 8 deletions.
diff --git a/src/main/java/picard/arrays/GtcToVcf.java b/src/main/java/picard/arrays/GtcToVcf.java
@@ -164,6 +164,8 @@ public class GtcToVcf extends CommandLineProgram {
         df.setGroupingSize(0);
     }
 
+    int NumZeroedOutAssays = 0;
+
     @Override
     protected boolean requiresReference() {
         return true;
@@ -192,6 +194,22 @@ protected int doWork() {
             // fill the sorting collection
             fillContexts(contexts, infiniumGTCFile, manifest, infiniumEGTFile);
 
+            final int numCalls = infiniumGTCFile.getNumCalls();
+            final int numSnps = infiniumGTCFile.getNumberOfSnps();
+            final double callRate = ((double) numCalls) / (numSnps - NumZeroedOutAssays);
+
+            log.debug("NumCalls:    " + numCalls);
+            log.debug("NumNoCalls:  " + infiniumGTCFile.getNumNoCalls());
+            log.debug("NumSnps:     " + numSnps);
+            log.debug("Zeroed:      " + NumZeroedOutAssays);
+            log.debug("IntOnly:     " + infiniumGTCFile.getNumIntensityOnly());
+            log.debug("GtcCallRate: " + infiniumGTCFile.getCallRate());
+            log.debug("CallRate:    " + callRate);
+
+            // The Call Rate in the Illumina GTC File does not take into account zeroed out SNPs.
+            // So we recalculate it (and ignore zeroed out assays - which also includes intensity only probes)
+            vcfHeader.addMetaDataLine(new VCFHeaderLine(InfiniumVcfFields.GTC_CALL_RATE, String.valueOf(callRate)));
+
             writeVcf(contexts, OUTPUT, refSeq.getSequenceDictionary(), vcfHeader);
 
             // not sure this is needed but it's definitely **cleaner**
@@ -290,9 +308,17 @@ private void fillContexts(final SortingCollection<VariantContext> contexts, fina
         while (iterator.hasNext()) {
             final Build37ExtendedIlluminaManifestRecord record = iterator.next();
 
+            Integer egtIndex = egtFile.rsNameToIndex.get(record.getName());
+            if (egtIndex == null) {
+                throw new PicardException("Found no record in cluster file for manifest entry '" + record.getName() + "'");
+            }
+            if (egtFile.totalScore[egtIndex] == 0.0) {
+                NumZeroedOutAssays++;
+            }
+
             if (!record.isFail()) {
                 InfiniumGTCRecord gtcRecord = gtcFile.getRecord(gtcIndex);
-                VariantContext context = makeVariantContext(record, gtcRecord, egtFile, progressLogger);
+                VariantContext context = makeVariantContext(record, gtcRecord, egtFile, egtIndex, progressLogger);
                 numVariantsWritten++;
                 contexts.add(context);
             }
@@ -301,23 +327,20 @@ private void fillContexts(final SortingCollection<VariantContext> contexts, fina
 
         log.info(numVariantsWritten + " Variants were written to file");
         log.info(gtcFile.getNumberOfSnps() + " Variants in the GTC file");
+        log.info(NumZeroedOutAssays + " Assays were zeroed out in the EGT file");
         log.info(manifest.getNumAssays() + " Variants on the " + manifest.getDescriptorFileName() + " genotyping array manifest file");
     }
 
     private VariantContext makeVariantContext(Build37ExtendedIlluminaManifestRecord record, final InfiniumGTCRecord gtcRecord,
-                                              final InfiniumEGTFile egtFile, final ProgressLogger progressLogger) {
+                                              final InfiniumEGTFile egtFile, final int egtIndex, final ProgressLogger progressLogger) {
         // If the record is not flagged as errant in the manifest we include it in the VCF
         Allele A = record.getAlleleA();
         Allele B = record.getAlleleB();
         Allele ref = record.getRefAllele();
 
         final String chr = record.getB37Chr();
         final Integer position = record.getB37Pos();
-        final Integer endPosition = position + ref.length() - 1;
-        Integer egtIndex = egtFile.rsNameToIndex.get(record.getName());
-        if (egtIndex == null) {
-            throw new PicardException("Found no record in cluster file for manifest entry '" + record.getName() + "'");
-        }
+        final int endPosition = position + ref.length() - 1;
 
         progressLogger.record(chr, position);
 
@@ -554,7 +577,6 @@ private VCFHeader createVCFHeader(final Build37ExtendedIlluminaManifest manifest
             lines.add(new VCFHeaderLine(controlInfo.getControl(), controlInfo.toString() + "|" + redIntensity + "|" + greenIntensity));
         }
         lines.add(new VCFHeaderLine(InfiniumVcfFields.FINGERPRINT_GENDER, fingerprintGender.name()));
-        lines.add(new VCFHeaderLine(InfiniumVcfFields.GTC_CALL_RATE, String.valueOf(gtcFile.getCallRate())));
         if (gtcGender != null) {
             lines.add(new VCFHeaderLine(InfiniumVcfFields.AUTOCALL_GENDER, gtcGender));
         } else {