Add the ability to emit the original quals in the OQ tag

protected/java/test/org/broadinstitute/sting/gatk/walkers/bqsr/BQSRIntegrationTest.java

@@ -107,7 +107,8 @@ public class BQSRIntegrationTest extends WalkerTest {
                 {new PRTest("", "d2d6ed8667cdba7e56f5db97d6262676")},
                 {new PRTest(" -qq -1", "b7053d3d67aba6d8892f0a60f0ded338")},
                 {new PRTest(" -qq 6", "bfbf0855185b2b70aa35237fb71e4487")},
-                {new PRTest(" -DIQ", "66aa65223f192ee39c1773aa187fd493")}
+                {new PRTest(" -DIQ", "66aa65223f192ee39c1773aa187fd493")},
+                {new PRTest(" -EOQ", "")}
         };
     }
 

public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java

@@ -198,8 +198,8 @@ public class GenomeAnalysisEngine {
     private BaseRecalibration baseRecalibration = null;
     public BaseRecalibration getBaseRecalibration() { return baseRecalibration; }
     public boolean hasBaseRecalibration() { return baseRecalibration != null; }
-    public void setBaseRecalibration(final File recalFile, final int quantizationLevels, final boolean disableIndelQuals, final int preserveQLessThan) {
+    public void setBaseRecalibration(final File recalFile, final int quantizationLevels, final boolean disableIndelQuals, final int preserveQLessThan, final boolean emitOriginalQuals) {
-        baseRecalibration = new BaseRecalibration(recalFile, quantizationLevels, disableIndelQuals, preserveQLessThan);
+        baseRecalibration = new BaseRecalibration(recalFile, quantizationLevels, disableIndelQuals, preserveQLessThan, emitOriginalQuals);
     }
 
     /**
@@ -239,7 +239,7 @@ public class GenomeAnalysisEngine {
 
         // if the use specified an input BQSR recalibration table then enable on the fly recalibration
         if (args.BQSR_RECAL_FILE != null)
-            setBaseRecalibration(args.BQSR_RECAL_FILE, args.quantizationLevels, args.disableIndelQuals, args.PRESERVE_QSCORES_LESS_THAN);
+            setBaseRecalibration(args.BQSR_RECAL_FILE, args.quantizationLevels, args.disableIndelQuals, args.PRESERVE_QSCORES_LESS_THAN, args.emitOriginalQuals);
 
         // Determine how the threads should be divided between CPU vs. IO.
         determineThreadAllocation();

public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java

@@ -220,6 +220,12 @@ public class GATKArgumentCollection {
     @Argument(fullName="disable_indel_quals", shortName = "DIQ", doc = "If true, disables printing of base insertion and base deletion tags (with -BQSR)", required=false)
     public boolean disableIndelQuals = false;
 
+    /**
+     * By default, the OQ tag in not emitted when using the -BQSR argument.
+     */
+    @Argument(fullName="emit_original_quals", shortName = "EOQ", doc = "If true, enables printing of the OQ tag with the original base qualities (with -BQSR)", required=false)
+    public boolean emitOriginalQuals = false;
+
     /**
      * Do not modify quality scores less than this value but rather just write them out unmodified in the recalibrated BAM file.
      * In general it's unsafe to change qualities scores below < 6, since base callers use these values to indicate random or bad bases.

public/java/src/org/broadinstitute/sting/utils/recalibration/BaseRecalibration.java

@@ -25,11 +25,14 @@
 
 package org.broadinstitute.sting.utils.recalibration;
 
+import net.sf.samtools.SAMTag;
+import net.sf.samtools.SAMUtils;
 import org.broadinstitute.sting.gatk.walkers.bqsr.*;
 import org.broadinstitute.sting.utils.MathUtils;
 import org.broadinstitute.sting.utils.QualityUtils;
 import org.broadinstitute.sting.utils.collections.NestedIntegerArray;
 import org.broadinstitute.sting.utils.collections.NestedHashMap;
+import org.broadinstitute.sting.utils.exceptions.UserException;
 import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
 
 import java.io.File;
@@ -51,6 +54,7 @@ public class BaseRecalibration {
 
     private final boolean disableIndelQuals;
     private final int preserveQLessThan;
+    private final boolean emitOriginalQuals;
 
     private static final NestedHashMap[] qualityScoreByFullCovariateKey = new NestedHashMap[EventType.values().length]; // Caches the result of performSequentialQualityCalculation(..) for all sets of covariate values.
     static {
@@ -66,7 +70,7 @@ public class BaseRecalibration {
      * @param disableIndelQuals  if true, do not emit base indel qualities
      * @param preserveQLessThan  preserve quality scores less than this value
      */
-    public BaseRecalibration(final File RECAL_FILE, final int quantizationLevels, final boolean disableIndelQuals, final int preserveQLessThan) {
+    public BaseRecalibration(final File RECAL_FILE, final int quantizationLevels, final boolean disableIndelQuals, final int preserveQLessThan, final boolean emitOriginalQuals) {
         RecalibrationReport recalibrationReport = new RecalibrationReport(RECAL_FILE);
 
         recalibrationTables = recalibrationReport.getRecalibrationTables();
@@ -80,6 +84,7 @@ public class BaseRecalibration {
         readCovariates = new ReadCovariates(MAXIMUM_RECALIBRATED_READ_LENGTH, requestedCovariates.length);
         this.disableIndelQuals = disableIndelQuals;
         this.preserveQLessThan = preserveQLessThan;
+        this.emitOriginalQuals = emitOriginalQuals;
     }
 
     /**
@@ -90,6 +95,14 @@ public class BaseRecalibration {
      * @param read the read to recalibrate
      */
     public void recalibrateRead(final GATKSAMRecord read) {
+        if (emitOriginalQuals && read.getAttribute(SAMTag.OQ.name()) == null) { // Save the old qualities if the tag isn't already taken in the read
+            try {
+                read.setAttribute(SAMTag.OQ.name(), SAMUtils.phredToFastq(read.getBaseQualities()));
+            } catch (IllegalArgumentException e) {
+                throw new UserException.MalformedBAM(read, "illegal base quality encountered; " + e.getMessage());
+            }
+        }
+
         RecalDataManager.computeCovariates(read, requestedCovariates, readCovariates);                                  // compute all covariates for the read
         for (final EventType errorModel : EventType.values()) {                                                         // recalibrate all three quality strings
             if (disableIndelQuals && errorModel != EventType.BASE_SUBSTITUTION) {