BQSR: avoid throwing an error when read group is missing in the recal table, and some refactoring. #9020

takutosato · 2024-10-28T13:24:47Z

Addresses #6242.

Current behavior: when all the reads in a read group are filtered in the base recalibration step, the read group is not logged in the recal table. Then ApplyBQSR encounters these reads, can't find the read group in the recal table, and throws an error.

New behavior: if --allow-read-group flag is set to true, then ApplyBQSR outputs the original quantities (after quantizing).

I avoided the alternative approach of collapsing (marginalizing) across the read groups, mostly because it would require a complete overhaul of the code. I also think that using recal data from other read groups might not be a good idea. In any case, using OQ should be good enough; I assume that these "missing" read groups are low enough quality to be filtered out and are likely to be thrown out by downstream tools.

I also refactored the BQSR code, mostly to update the variable and class names to be more accurate and descriptive. For instance:

ReadCovariates.java -> PerReadCovariateMatrix.java
EstimatedQReported -> ReportedQuality

davidbenjamin

Looks good, just needs another pass at some of your comments and TODOs.

davidbenjamin · 2024-11-07T15:57:47Z

build.gradle

@@ -188,7 +188,7 @@ configurations.all {
 }

 tasks.withType(JavaCompile) {
-  options.compilerArgs = ['-proc:none', '-Xlint:all', '-Werror', '-Xdiags:verbose']
+  options.compilerArgs = ['-proc:none', '-Xlint:all', '-Xdiags:verbose']


This is beyond my job description but I trust there's a reason?

davidbenjamin · 2024-11-07T16:04:43Z

src/main/java/org/broadinstitute/hellbender/tools/ApplyBQSRArgumentCollection.java

+     * If set to true, do not throw an error upon encountering a read with a read group that's not in the recalibration table.
+     * Instead, simply set the quantized original base qualities as the recalibrated base qualities.
+     */ // tsato: should this be in the *unique* ApplyBQSRArgumentCollection?
+    @Argument(fullName = ALLOW_MISSING_READ_GROUPS_LONG_NAME, doc = "", optional = true)


empty doc string?

davidbenjamin · 2024-11-14T17:11:21Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/bqsr/ApplyBQSR.java

@@ -90,7 +90,7 @@ public final class ApplyBQSR extends ReadWalker{
     */
    @Argument(fullName=StandardArgumentDefinitions.BQSR_TABLE_LONG_NAME, shortName=StandardArgumentDefinitions.BQSR_TABLE_SHORT_NAME, doc="Input recalibration table for BQSR")
    @WorkflowInput
-    public File BQSR_RECAL_FILE;
+    public File bqsrRecalFile; // tsato: change to GATKPath?


Address this TODO

davidbenjamin · 2024-11-14T17:11:50Z

src/main/java/org/broadinstitute/hellbender/transformers/BQSRReadTransformer.java

@@ -26,14 +28,14 @@

 public final class BQSRReadTransformer implements ReadTransformer {
    private static final long serialVersionUID = 1L;
-
+    // tsato: remove?


davidbenjamin · 2024-11-14T17:15:01Z

src/main/java/org/broadinstitute/hellbender/transformers/BQSRReadTransformer.java


-        //clear indel qualities
+        // clear indel qualities // tsato: why? Is this ok? Do we update them?


davidbenjamin · 2024-11-14T17:22:59Z

src/main/java/org/broadinstitute/hellbender/utils/recalibration/BaseRecalibrationEngine.java

        for ( final NestedIntegerArray.Leaf<RecalDatum> leaf : byQualTable.getAllLeaves() ) {
-            final int rgKey = leaf.keys[0];
-            final int eventIndex = leaf.keys[2];
+            final int rgKey = leaf.keys[readGroupIndex]; // tsato: ??? What are these indices?


davidbenjamin · 2024-11-14T17:24:39Z

src/main/java/org/broadinstitute/hellbender/utils/recalibration/RecalUtils.java

@@ -32,7 +32,7 @@
 public final class RecalUtils {
    public static final String ARGUMENT_REPORT_TABLE_TITLE = "Arguments";
    public static final String QUANTIZED_REPORT_TABLE_TITLE = "Quantized";
-    public static final String READGROUP_REPORT_TABLE_TITLE = "RecalTable0";
+    public static final String READGROUP_REPORT_TABLE_TITLE = "RecalTable0"; // TODO: make them more descriptive


davidbenjamin · 2024-11-14T17:24:58Z

src/main/java/org/broadinstitute/hellbender/utils/recalibration/RecalUtils.java

    private static boolean warnUserNullPlatform = false;

    private static final String SCRIPT_FILE = "BQSR.R";
    public static final int EMPIRICAL_QUAL_DECIMAL_PLACES = 4;
-    public static final int EMPIRICAL_Q_REPORTED_DECIMAL_PLACES = 4;
+    public static final int REPORTED_QUALITY_DECIMAL_PLACES = 4; // tsato: "estimated" q reported...we need to rename (DONE)


davidbenjamin · 2024-11-14T17:25:59Z

...main/java/org/broadinstitute/hellbender/utils/recalibration/covariates/ContextCovariate.java

 public final class ContextCovariate implements Covariate {
    private static final long serialVersionUID = 1L;
    private static final Logger logger = LogManager.getLogger(ContextCovariate.class);

-    private final int mismatchesContextSize;
+    private final int mismatchesContextSize; // TODO: rename "mismatch" here


davidbenjamin · 2024-11-14T17:30:45Z

src/main/java/org/broadinstitute/hellbender/utils/recalibration/RecalDatum.java

@@ -136,27 +152,30 @@ public double getEmpiricalErrorRate() {
            return doubleMismatches / doubleObservations;
        }
    }
-
-    public void setEmpiricalQuality(final double empiricalQuality) {
+    // tsato: int-double


clean up comment

lbergelson

@takutosato I took a look. Definitely want to put back the gradle build to how it was. If you need help with the warning we can help. I'm not sure I totally understand the motivation for changing the internal quality value from a double to an int? I think you changed if to being PHRED scaled, but it seems like that might change the precision of the calculation which could change the results? Or is it always output as a PHRED value before further calculations are made? It's not clear to me.

Thank you for the clarifying comments.

Lots of unaddressed TODOs. Some of them probably should be removed but if you think that some of them are useful to future readers than that's fine.

lbergelson · 2024-10-30T21:15:11Z

build.gradle

@@ -188,7 +188,7 @@ configurations.all {
 }

 tasks.withType(JavaCompile) {
-  options.compilerArgs = ['-proc:none', '-Xlint:all', '-Werror', '-Xdiags:verbose']
+  options.compilerArgs = ['-proc:none', '-Xlint:all', '-Xdiags:verbose']


What warning are we seeing that made you take this out?

Lets fix the warning or supress it instead of turning this off.

lbergelson · 2024-10-30T21:17:53Z

src/main/java/org/broadinstitute/hellbender/tools/ApplyBQSRArgumentCollection.java

+     * If set to true, do not throw an error upon encountering a read with a read group that's not in the recalibration table.
+     * Instead, simply set the quantized original base qualities as the recalibrated base qualities.
+     */ // tsato: should this be in the *unique* ApplyBQSRArgumentCollection?
+    @Argument(fullName = ALLOW_MISSING_READ_GROUPS_LONG_NAME, doc = "", optional = true)


Forgot the doc string.

Thanks for this note. I think this argument should be moved up. to ApplyBQSRUniqueArgumentCollection. It is only relevant to the apply stage and not to the recalibration stage so that makes sense for it to be there. Also, I think you need to update the method
ApplyBQSRUniqueArgumentCollection.toApplyBQSRArgumentCollection to pass this argument to the newly created collection.

lbergelson · 2024-10-30T21:19:31Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/bqsr/ApplyBQSR.java

@@ -90,7 +90,7 @@ public final class ApplyBQSR extends ReadWalker{
     */
    @Argument(fullName=StandardArgumentDefinitions.BQSR_TABLE_LONG_NAME, shortName=StandardArgumentDefinitions.BQSR_TABLE_SHORT_NAME, doc="Input recalibration table for BQSR")
    @WorkflowInput
-    public File BQSR_RECAL_FILE;
+    public File bqsrRecalFile; // tsato: change to GATKPath?


That would be a good thing to do but we could do it in a separate PR.

lbergelson · 2024-10-30T21:27:11Z

src/main/java/org/broadinstitute/hellbender/transformers/BQSRReadTransformer.java

@@ -26,14 +28,14 @@

 public final class BQSRReadTransformer implements ReadTransformer {
    private static final long serialVersionUID = 1L;
-
+    // tsato: remove?


I think it can be removed.

lbergelson · 2024-10-30T21:30:16Z

src/main/java/org/broadinstitute/hellbender/transformers/BQSRReadTransformer.java

    private byte[] staticQuantizedMapping;
    private final CovariateKeyCache keyCache;

+    // tsato: new flag, needs to be moved to apply BQSR argument
+    private boolean allowMissingReadGroup;


allowMissingReadGroup -> allowMissingReadGroups

make this final

I think the comment is resolved.

lbergelson · 2024-11-18T19:56:27Z

src/test/java/org/broadinstitute/hellbender/GATKBaseTest.java

@@ -92,6 +92,8 @@ public abstract class GATKBaseTest extends BaseTest {
    // Variants from a DBSNP 138 VCF form the first 65Mb of chr1
    public static final String dbsnp_138_b37_1_65M_vcf = largeFileTestDir + "dbsnp_138.b37.1.1-65M.vcf";

+    public static final String BQSR_TEST_RESOURCE_DIR = toolsTestDir + "BQSR/";


This seems like the wrong place for this constant. I'd put it in a BQSR test or maybe just use the CommandLineProgram.getTestDataDir() if that works.

lbergelson · 2024-11-18T19:57:18Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/bqsr/ApplyBQSRIntegrationTest.java

-                IllegalStateException.class);
-        spec.executeTest("testMissingReadGroup", this);
+    public void testMissingReadGroup() {
+        // TODO: there are two copies of this file (one under Validation and one under BQSR) in the repo; conslidate.


That's fine. Duplicated test files are alright. Git deduplicates them internally for transfer so they don't really take up meaningful amounts of room.

lbergelson · 2024-11-18T19:58:30Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/bqsr/ApplyBQSRIntegrationTest.java

+            final byte[] oldQuals = ReadUtils.getOriginalBaseQualities(originalRead);
+
+            if (ReadUtils.getPlatformUnit(originalRead, originalBamHeader).equals(readGroupToFilterOut)) {
+                // These are the read groups that are not in teh recal table.


lbergelson · 2024-11-18T19:58:53Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/bqsr/ApplyBQSRIntegrationTest.java

+
+            if (ReadUtils.getPlatformUnit(originalRead, originalBamHeader).equals(readGroupToFilterOut)) {
+                // These are the read groups that are not in teh recal table.
+                // final Random random = new Random();


Can this be deleted?

lbergelson · 2024-11-18T20:00:35Z

src/test/java/org/broadinstitute/hellbender/utils/recalibration/RecalDatumUnitTest.java

@@ -97,12 +97,12 @@ public void testRecalDatumCopyAndCombine(RecalDatumTestProvider cfg) {
    @Test(dataProvider = "RecalDatumTestProvider")
    public void testRecalDatumModification(RecalDatumTestProvider cfg) {
        RecalDatum datum = cfg.makeRecalDatum();
-        datum.setEmpiricalQuality(10.1);
-        Assert.assertEquals(datum.getEmpiricalQuality(), 10.1);
+        datum.setEmpiricalQuality(10); // tsato: 10.1 -> 10. Downstream data might be affected.


Does changing the quality from double -> int potentially change the results of the recalibration table by rounding earlier? I'm not sure what the advantage of that is?

takutosato · 2024-11-19T00:15:46Z

@lbergelson thanks for the review! I will address the changes as soon as I can.

The motivation for the double -> int change is in RecalDatum::bayesianEstimateOfEmpiricalQuality(). In the master branch now, this method always returns a (phred-scaled integer) / RESOLUTION_BINS_PER_QUAL, and RESOLUTION_BINS_PER_QUAL is a final constant set to 1. So this method really returns an integer, and my change attempts to make this clear.

gatk-bot · 2024-11-19T00:27:17Z

Github actions tests reported job failures from actions build 11903509453
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
cloud	17.0.6+10	11903509453.10	logs
unit	17.0.6+10	11903509453.12	logs
integration	17.0.6+10	11903509453.11	logs
unit	17.0.6+10	11903509453.1	logs
integration	17.0.6+10	11903509453.0	logs

takutosato added 8 commits August 9, 2024 08:11

temp

1a574ac

commit and push data before it's too late per Louis

409f928

temporarily disable werror but will need to figure it out at some point

3a213b5

first draft implementation

b0a6019

test passing

4660b00

proof reading 1

499e0e7

proofreading 2

4dd794b

finished proofreading

d372f01

davidbenjamin self-requested a review November 4, 2024 14:55

davidbenjamin self-assigned this Nov 4, 2024

davidbenjamin reviewed Nov 14, 2024

View reviewed changes

review 1

07f5383

lbergelson reviewed Nov 18, 2024

View reviewed changes

takutosato added 2 commits November 18, 2024 19:08

another review checkpoint

c90f3fc

revert build.gradle

109ee36

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BQSR: avoid throwing an error when read group is missing in the recal table, and some refactoring. #9020

BQSR: avoid throwing an error when read group is missing in the recal table, and some refactoring. #9020

takutosato commented Oct 28, 2024

davidbenjamin left a comment

davidbenjamin Nov 7, 2024

davidbenjamin Nov 7, 2024

takutosato Nov 18, 2024

davidbenjamin Nov 14, 2024

takutosato Nov 18, 2024

davidbenjamin Nov 14, 2024

takutosato Nov 18, 2024

davidbenjamin Nov 14, 2024

takutosato Nov 18, 2024

davidbenjamin Nov 14, 2024

takutosato Nov 18, 2024

davidbenjamin Nov 14, 2024

davidbenjamin Nov 14, 2024

davidbenjamin Nov 14, 2024

davidbenjamin Nov 14, 2024

lbergelson left a comment

lbergelson Oct 30, 2024

lbergelson Oct 30, 2024

lbergelson Oct 30, 2024

lbergelson Oct 30, 2024

lbergelson Oct 30, 2024

lbergelson Oct 30, 2024

lbergelson Oct 30, 2024

lbergelson Nov 18, 2024

lbergelson Nov 18, 2024

lbergelson Nov 18, 2024

lbergelson Nov 18, 2024

lbergelson Nov 18, 2024

takutosato commented Nov 19, 2024

gatk-bot commented Nov 19, 2024 •

edited

Loading


		//clear indel qualities
		// clear indel qualities // tsato: why? Is this ok? Do we update them?

BQSR: avoid throwing an error when read group is missing in the recal table, and some refactoring. #9020

Are you sure you want to change the base?

BQSR: avoid throwing an error when read group is missing in the recal table, and some refactoring. #9020

Conversation

takutosato commented Oct 28, 2024

davidbenjamin left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lbergelson left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

takutosato commented Nov 19, 2024

gatk-bot commented Nov 19, 2024 • edited Loading

gatk-bot commented Nov 19, 2024 •

edited

Loading