SelectVariants JEXL filter fixes and refactor #8092

takutosato · 2022-11-14T16:36:08Z

This PR is an attempt at improving SelectVariants by

Rewriting unclear code;
Adding documentation where needed; and
Adding tests.

The initial motivation for this code change (#7497) was to improve performance of SelectVariants by adding an option to do the "INFO-level filtering" before doing "genotype filtering." Our assumption was that this would improve performance because then we would avoid the expensive genotype "fully-decode" operation, which turns string format fields into appropriate object/types (int, array, etc.). This is (we think) done in VariantContext.fullyDecode().

This turned out not to be possible for the following reasons. First, there are roughly four types of genotype subsetting you could do:

a) By the sample names (--sample-name NA12878)
b) JEXL (--select GQ > 0)
c) JEXL by accessing the variant context object (--select vc.getGenotype('NA12878').getGQ() > 1)
d) Others (e.g. --remove-fraction-genotype)

a) does not need "fully-decode." It turns out b) was never supported (GATK currently removes all variants and succeed.) And from my experiments, c) does not seem to ever trigger calling VariantContext.fullyDecode(). In fact the only code path I can see that calls fullyDecode() is by setting the fully-decode SelectVariants argument, which seems to just call fullyDecode at the beginning just for the sake of calling it (or so it appears to me. The utility of this command line argument is highly dubious.)

It's possible that apache code does something similar to fully decoding that could affect performance.

All that is to say that we cannot achieve performance improvement with our original blueprint simply because this expensive "fullyDecode" operation seems to be a mythical operation that is never used in reality.

So while I could not speed up SelectVariants, I cleaned up the code and added the following new arguments:

--select-genotype: with this new genotype-specific JEXL argument, we support filtering by genotype fields like 'GQ > 0', where the behavior in the multi-sample case is 'GQ > 0' in at least one sample. I have not added the ability to do 'GQ > 0 for all samples' but it should be a simple (but not easy…) exercise in boolean operations.
applyJexlFiltersBeforeFilteringGenotypes: if set to true, we do the JEXL checking before we subset by samples. In my tests, performance improvement from this option was very modest. Subsetting a ~3k 1kg SV vcf to a single sample was about 30 seconds faster (out of ~20 min total run time) than the default. I kept it in the PR because I thought some user might find it useful, but I wouldn't be opposed to removing it.

Tests needed:

takutosato · 2022-11-14T16:40:33Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+        // START Chris/Louis/David
+        // Attempt 1, does not work
+//        final ArgumentsBuilder args1 = new ArgumentsBuilder()
+//                .add("V", testVcf.toString())
+//                .add("O", testOutput.getAbsolutePath())
+//                .add("--select", "'AC > 1'")
+//                .add("--select", "'AF > 0'");
+//        // .add("select", "'GQ > 0'"); // tsato: this does not work for some reason...
+//        runCommandLine(args1, SelectVariants.class.getSimpleName());
+
+        // Attempt 2 (remove single quotes), does not work
+//        final ArgumentsBuilder args2 = new ArgumentsBuilder()
+//                .add("V", testVcf.toString())
+//                .add("O", testOutput.getAbsolutePath())
+//                .add("--select", "AC > 1")
+//                .add("--select", "AF > 0");
+//        // .add("select", "'GQ > 0'"); // tsato: this does not work for some reason...
+//        runCommandLine(args2, SelectVariants.class.getSimpleName());
+
+        // Atempt 3 (remove spaces). commandline parser lets this through but not the JEXL parser
+//        final ArgumentsBuilder args3 = new ArgumentsBuilder()
+//                .add("V", testVcf.toString())
+//                .add("O", testOutput.getAbsolutePath())
+//                .addRaw("--select 'AC>1'")
+//                .addRaw("--select 'AF>0'");
+//                // .add("select", "'GQ > 0'"); // tsato: this does not work for some reason...
+//        runCommandLine(args3, SelectVariants.class.getSimpleName());
+
+        // Attempt 4 (give a list to runCommandLine), works
+//        runCommandLine(Arrays.asList("-V", testVcf.toString(),
+//                        "-O", testOutput.getAbsolutePath(),
+//                        "--select","AC > 1"),
+//                SelectVariants.class.getSimpleName());
+
+        // Attempt 5, works!
+//        runCommandLine(Arrays.asList("-V", testVcf.toString(),
+//                        "-O", testOutput.getAbsolutePath(),
+//                        "--select", "ALGORITHMS == 'depth'"),
+//                SelectVariants.class.getSimpleName());
+
+        // Attempt 6, works!
+        runCommandLine(Arrays.asList("-V", testVcf.toString(),
+                        "-O", testOutput.getAbsolutePath(),
+                        "--select", "SVTYPE == 'DEL'"),
+                SelectVariants.class.getSimpleName());
+        // END Chris/Louis/David


@droazen @lbergelson @cmnbroad this block of code has the command line parsing issue I mentioned. My work around is Attempt 5 and 6, where I feed runCommandLine a list of arguments rather than an argument builder. Is that good enough? Otherwise I suspect we need to update ArgumentsBuilder to support JEXL command lines.

gatk-bot · 2022-11-14T16:56:09Z

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

Test Type	JDK	Job ID	Logs
cloud	8	3463372954.10	logs
cloud	11	3463372954.11	logs
integration	11	3463372954.12	logs
integration	8	3463372954.0	logs

gatk-bot · 2022-11-19T22:07:16Z

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

Test Type	JDK	Job ID	Logs
cloud	11	3505269846.11	logs
unit	11	3505269846.13	logs
integration	11	3505269846.12	logs
integration	8	3505269846.0	logs

gatk-bot · 2022-11-28T18:40:36Z

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

Test Type	JDK	Job ID	Logs
unit	11	3567484153.13	logs
integration	11	3567484153.12	logs

gatk-bot · 2022-11-28T18:55:09Z

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

Test Type	JDK	Job ID	Logs
unit	11	3567534442.13	logs
integration	11	3567534442.12	logs

droazen

@takutosato Back to you with my comments. Changes look really good overall, and definitely an improvement over the existing tool! Go ahead and merge this after addressing my comments and tests pass. You'll need to rebase the branch to get rid of some unrelated test failures.

droazen · 2022-12-06T04:10:10Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

-     * specified samples are extracted and the INFO field annotations are updated.
+     * See example commands above for detailed usage examples. The expressions given to this argument
+     * should either refer to INFO fields, or access FORMAT field with the VariantContext object
+     * e.g. --select "vc.getGenotype('NA12878').getGQ() > 35"


It might be a good idea to include an example with vc.getGenotype().getGQ() in the main tool docs above, along with a comment comparing that syntax to the new -select-genotype approach.

Also, --select should be -select here

droazen · 2022-12-06T04:23:58Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

-    @Argument(shortName="select", doc="One or more criteria to use when selecting the data", optional=true)
+    public static final String SELECT_SHORT_NAME = "select";
+    public static final String SELECT_LONG_NAME = "select-expressions";
+    @Argument(fullName=SELECT_LONG_NAME, shortName=SELECT_SHORT_NAME, doc="", optional=true)


Doc string should not be empty for this argument

A bit worried about breaking some users' script by not allowing --select, so I might revert this new long name change.

droazen · 2022-12-06T04:33:32Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

+     */
+    public static final String GENOTYPE_SELECT_SHORT_NAME = "select-genotype";
+    public static final String GENOTYPE_SELECT_LONG_NAME = "select-genotype-expressions";
+    @Argument(fullName=GENOTYPE_SELECT_LONG_NAME, shortName=GENOTYPE_SELECT_SHORT_NAME, doc="", optional=true)


Doc string should not be empty

droazen · 2022-12-06T04:35:23Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

+     * on INFO fields only, which may improve speed when working with a large cohort vcf that contains genotypes for
+     * thousands of samples (format fields).
+     */
+    @Argument(fullName="jexl-first", shortName="jexl-first", doc="", optional=true)


Fill in doc string

Long argument name should perhaps be more descriptive here -- eg., --apply-jexl-filters-first

droazen · 2022-12-06T04:38:30Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

@@ -480,7 +523,7 @@ public void onTraversalStart() {

        final List<String> genotypeField = getHeaderForVariants().getGenotypeSamples();
        if(!ParsingUtils.isSorted(genotypeField)){
-            if(genotypeField.size() > 10) {
+            if(!genotypeField.isEmpty()) {


I believe that the previous threshold of 10 was deliberate, so that we don't warn unnecessarily about performance when there are only a few samples.

Sounds good, done.

droazen · 2022-12-06T05:41:36Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+        final int threshold = 1;
+        final Predicate<VariantContext> filter = vc -> vc.getGenotypes().stream().anyMatch(g -> g.getGQ() > threshold);
+        final Pair<VCFHeader, List<VariantContext>> preFilter = VariantContextTestUtils.readEntireVCFIntoMemory(svHGDBMultiSampleVcf.toString());
+        final int expectedNumVCs = (int) preFilter.getRight().stream().filter(filter).count();


In these new tests you added, have you confirmed in each case that the number of expected records is > 0? Otherwise the tests will not be as meaningful as they could be.

Added extra assert statements to check this.

droazen · 2022-12-06T05:48:02Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+                SelectVariants.class.getSimpleName());
+        final Pair<VCFHeader, List<VariantContext>> result2 = VariantContextTestUtils.readEntireVCFIntoMemory(testOutput2.getAbsolutePath());
+        final int actualNumVCs2 = result2.getRight().size();
+        Assert.assertEquals(expectedNumVCs2, actualNumVCs2);


This second test case should probably be a separate test method rather than embedded within testCombineInfoAndGenotypeJexl()

droazen · 2022-12-06T05:48:36Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+                SelectVariants.class.getSimpleName());
+        final Pair<VCFHeader, List<VariantContext>> result3 = VariantContextTestUtils.readEntireVCFIntoMemory(testOutput3.getAbsolutePath());
+        final int actualNumVCs3 = result3.getRight().size();
+        Assert.assertEquals(0, actualNumVCs3);


This should also be pulled out into a separate named test method of its own.

droazen · 2022-12-06T05:49:03Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+
+
+
+        // Assert.assertEquals(expectedNumVCs2, actualNumVCs2);D


Delete commented-out code

droazen · 2022-12-06T05:49:49Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+//        runCommandLine(Arrays.asList("-V", testVcf.toString(),
+//                        "-O", testOutput.getAbsolutePath(),
+//                        "--select", "ALGORITHMS == 'depth'"),
+//                SelectVariants.class.getSimpleName());


Delete commented-out code (but consult with @lbergelson on why ArgumentsBuilder didn't work here)

takutosato · 2023-01-05T16:46:23Z

.../org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariantsIntegrationTest.java

+        final Predicate<VariantContext> filter = vc -> vc.getAttributeAsInt("AC", -1) > 0 ||
+                vc.getAttributeAsString("SVTYPE", "").equals("DEL");
+
+        final Pair<VCFHeader, List<VariantContext>> preFilter = VariantContextTestUtils.readEntireVCFIntoMemory(svHGDBMultiSampleVcf.toString());
+        final int expectedCount = (int) preFilter.getRight().stream().filter(filter).count();
+
+        runCommandLine(Arrays.asList("-V", svHGDBMultiSampleVcf.toString(),
+                        "-O", testOutput.getAbsolutePath(),
+                        "--select", "SVTYPE == 'DEL'",
+                        "--select", "AC > 0"),


@cmnbroad what do you think? The issue here is that "--select" should no longer be a valid command line but runCommandLine does not complain about it.

takutosato · 2023-01-05T16:50:38Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

@@ -2,13 +2,7 @@

 import htsjdk.samtools.SAMSequenceDictionary;
 import htsjdk.tribble.util.ParsingUtils;
-import htsjdk.variant.variantcontext.Allele;


Put the explicit imports back

takutosato · 2023-01-05T17:03:00Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

@@ -603,6 +654,11 @@ public void apply(VariantContext vc, ReadsContext readsContext, ReferenceContext
            return;
        }

+        if (applyJexlFiltersBeforeFilteringGenotypes && !passesJexlFilters(vc)) {
+            return;


This should be tested.

codecov · 2023-01-07T19:21:00Z

Codecov Report

Merging #8092 (78c226d) into master (90bf668) will increase coverage by 0.023%.
The diff coverage is 92.237%.

❗ Current head 78c226d differs from pull request most recent head b0dcf74. Consider uploading reports for the commit b0dcf74 to get more accurate results

Additional details and impacted files

@@               Coverage Diff               @@
##              master     #8092       +/-   ##
===============================================
+ Coverage     86.636%   86.659%   +0.023%     
- Complexity     38944     39021       +77     
===============================================
  Files           2335      2337        +2     
  Lines         182675    182958      +283     
  Branches       20053     20080       +27     
===============================================
+ Hits          158262    158549      +287     
+ Misses         17373     17362       -11     
- Partials        7040      7047        +7

Impacted Files	Coverage Δ
...ender/tools/walkers/filters/VariantFiltration.java	`92.593% <ø> (ø)`
...der/tools/walkers/variantutils/SelectVariants.java	`83.911% <86.765%> (+2.516%)`	⬆️
...rs/variantutils/SelectVariantsIntegrationTest.java	`96.928% <94.702%> (-0.783%)`	⬇️
.../scalable/data/LabeledVariantAnnotationsDatum.java	`68.421% <0.000%> (-3.801%)`	⬇️
...rs/vqsr/scalable/TrainVariantAnnotationsModel.java	`77.778% <0.000%> (-2.991%)`	⬇️
...bender/utils/runtime/AsynchronousStreamWriter.java	`81.633% <0.000%> (-2.041%)`	⬇️
...ct/CreateSomaticPanelOfNormalsIntegrationTest.java	`96.396% <0.000%> (-1.305%)`	⬇️
...stitute/hellbender/utils/config/ConfigFactory.java	`73.750% <0.000%> (-1.250%)`	⬇️
...ools/walkers/annotator/VariantAnnotatorEngine.java	`86.260% <0.000%> (-0.999%)`	⬇️
...org/broadinstitute/hellbender/utils/MathUtils.java	`80.622% <0.000%> (-0.957%)`	⬇️
... and 80 more

takutosato added 9 commits October 24, 2022 19:39

select variants tests

7960a9d

clean up and add comment

5a651c2

just a commit

3f1cb46

cleaned up

8f00b6b

add

6e0bf02

documentation

298ebba

oops

8ae12bf

interim

027bc98

integratttion test issues

aac4af6

takutosato commented Nov 14, 2022

View reviewed changes

takutosato marked this pull request as draft November 19, 2022 20:15

all the tests

9da35dc

takutosato added 3 commits November 27, 2022 19:37

some changes

e93d9c2

proofreading

94ff401

ready

615ac1a

takutosato requested a review from droazen November 28, 2022 18:20

takutosato marked this pull request as ready for review November 28, 2022 18:21

droazen approved these changes Dec 6, 2022

View reviewed changes

droazen assigned takutosato Dec 6, 2022

takutosato commented Jan 5, 2023

View reviewed changes

takutosato added 2 commits January 7, 2023 12:40

David R review and optimizations from profiling

8d5db3e

rename

78c226d

last check

b0dcf74

takutosato merged commit 9f77b1f into master Jan 7, 2023

takutosato deleted the ts_select_variants branch January 7, 2023 20:03

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

SelectVariants JEXL filter fixes and refactor #8092

SelectVariants JEXL filter fixes and refactor #8092

takutosato commented Nov 14, 2022 •

edited

Loading

takutosato Nov 14, 2022

gatk-bot commented Nov 14, 2022 •

edited

Loading

gatk-bot commented Nov 19, 2022 •

edited

Loading

gatk-bot commented Nov 28, 2022 •

edited

Loading

gatk-bot commented Nov 28, 2022 •

edited

Loading

droazen left a comment

droazen Dec 6, 2022

droazen Dec 6, 2022

takutosato Jan 5, 2023

droazen Dec 6, 2022

takutosato Jan 5, 2023

takutosato Jan 5, 2023

droazen Dec 6, 2022

takutosato Jan 5, 2023

droazen Dec 6, 2022

droazen Dec 6, 2022

takutosato Jan 5, 2023

droazen Dec 6, 2022

takutosato Jan 5, 2023

droazen Dec 6, 2022

takutosato Jan 7, 2023

takutosato Jan 7, 2023

droazen Dec 6, 2022

takutosato Jan 7, 2023

droazen Dec 6, 2022

takutosato Jan 7, 2023

droazen Dec 6, 2022

takutosato Jan 7, 2023

droazen Dec 6, 2022

takutosato Jan 5, 2023

takutosato Jan 5, 2023

takutosato Jan 5, 2023

takutosato Jan 7, 2023

codecov bot commented Jan 7, 2023 •

edited

Loading

SelectVariants JEXL filter fixes and refactor #8092

SelectVariants JEXL filter fixes and refactor #8092

Conversation

takutosato commented Nov 14, 2022 • edited Loading

Choose a reason for hiding this comment

gatk-bot commented Nov 14, 2022 • edited Loading

gatk-bot commented Nov 19, 2022 • edited Loading

gatk-bot commented Nov 28, 2022 • edited Loading

gatk-bot commented Nov 28, 2022 • edited Loading

droazen 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

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

codecov bot commented Jan 7, 2023 • edited Loading

Codecov Report

takutosato commented Nov 14, 2022 •

edited

Loading

gatk-bot commented Nov 14, 2022 •

edited

Loading

gatk-bot commented Nov 19, 2022 •

edited

Loading

gatk-bot commented Nov 28, 2022 •

edited

Loading

gatk-bot commented Nov 28, 2022 •

edited

Loading

codecov bot commented Jan 7, 2023 •

edited

Loading