Update SV clustering classes #7243

mwalker174 · 2021-05-06T21:28:51Z

Reworks classes used by JointGermlineCNVSegmentationIntegration for SV clustering and defragmentation. The design of SVClusterEngine has been overhauled to enable the implementation of CNVDefragmenter and BinnedCNVDefragmenter subclasses. Logic for producing representative records from a collection of clustered SVs has been separated into an SVCollapser class, which provides enhanced functionality for handling genotypes for SVs more generally.

A number of bugs, particularly with max-clique clustering, have been fixed, as well as a parameter swap bug in JointGermlineCNVSegmentationIntegration.

This is the first of a series of PRs for an experimental Java-based implementation of some modules in gatk-sv pipeline, including SV vcf merging, clustering, evidence aggregation, and genotyping.

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java

ldgauthier

I can't find the code to do the defragmentation in bin-space that I need for exomes.

I also had a question regarding what the spec says about determining the ref allele.

Other than that, lots of little clarifications and suggestions you can take or leave. I didn't go over the tests with a fine toothed comb, but I can chip in if Chris doesn't have the time for that.

src/main/java/org/broadinstitute/hellbender/tools/sv/SVLocatable.java

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecordUtils.java

ldgauthier · 2021-06-09T17:42:48Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapser.java

+    protected Object collapseSingleVariantAttribute(final Set<Object> values) {
+        if (values.size() == 1) {
+            return values.iterator().next();
+        } else {


We don't want to add up DP or anything?

See comment above.

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapser.java

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/CNVDefragmenter.java

mwalker174 · 2021-06-30T20:40:40Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/BinnedCNVDefragmenter.java

+
+/**
+ * CNV defragmenter for when the intervals used for coverage collection are available. The intervals are used to adjust
+ * variant padding by rounding to the nearest bin boundary.


Is there something missing from getPaddedRecordInterval() below?

ldgauthier · 2021-06-09T19:47:57Z

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/CNVDefragmenterTest.java

+                call2End - call2Start + 1, Collections.singletonList(GATKSVVCFConstants.DEPTH_ALGORITHM),
+                Lists.newArrayList(Allele.REF_N, Allele.SV_SIMPLE_DEL),
+                Collections.emptyList(), Collections.emptyMap());
+        Assert.assertTrue(defragmenter.clusterTogether(call1, call2));


This is a nice improvement over the old checks.

cwhelan

This is great stuff! Looks really good. The tests are great, too. Feel free to take some of my comments as suggestions rather than requirements, or correct me if I'm misunderstood something as I'm commenting on it.

cwhelan · 2021-06-14T18:33:22Z

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java

+
+    // evidence metrics
+    public static final String COPY_NUMBER_LOG_POSTERIORS_KEY = "CNLP";
+    public static final String NEUTRAL_COPY_NUMBER_KEY = "NCN";


What's this one for? Is it to tag copy number neutral events? Is it meant to hold the reference ploidy?

It's for ploidy - future PR.

cwhelan · 2021-06-14T18:45:10Z

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java

+    public static final String CLUSTER_MEMBER_IDS_KEY = "MEMBERS";
+
+    // evidence metrics
+    public static final String COPY_NUMBER_LOG_POSTERIORS_KEY = "CNLP";


The new constants from this point on don't seem to be used anywhere in the GATK code base that I can see. Am I missing some way in which they are used? If not, do we want to add them now or hold off until they are used for something?

I tried to cherry-pick all the changes for this file, but these are actually for a later PR. I'll remove them for now.

cwhelan · 2021-06-14T18:45:54Z

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java

+    public static final String COPY_NUMBER_LOG_POSTERIORS_KEY = "CNLP";
+    public static final String NEUTRAL_COPY_NUMBER_KEY = "NCN";
+    public static final String DEPTH_OVERLAP_KEY = "RDOV";
+    public static final String DEPTH_P_HARDY_WEINBERG_LOSS_FIELD = "PHW_L";


I'm also not sure what some of these are. If you want to keep them, could you add some quick comments here or somewhere else to describe what they are?

cwhelan · 2021-06-14T18:49:20Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

-    @Override
-    public int getStart() {
-        return start;
+    public boolean isCNV() {


Some of the CPX types in our post 0506 VCFs are or contain CNVs. Maybe change this to something like isSimpleCNV?

cwhelan · 2021-06-14T18:56:45Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecordUtils.java

+     * @param attributes attributes to apply to all new genotypes
+     * @return genotypes augmented with missing samples
+     */
+    public static GenotypesContext fillMissingSamplesWithGenotypes(final GenotypesContext genotypes,


Just a suggestion, but I might change the name of this to populateGenotypesForMissingSamplesWithAlleles.

cwhelan · 2021-06-17T13:40:28Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapser.java

+        } else if (newType.equals(StructuralVariantType.INS)) {
+            // Median
+            final int[] lengths = items.stream().mapToInt(SVCallRecord::getLength).toArray();
+            final int midIndex = lengths.length / 2;


Would it make sense to use MathUtils.median()? You might have to update it or create another flavor of the method that would let you use the averaging EstimationType, but that might be more generally useful for GATK

I don't love the implementation but I suppose I should use the existing utility. If it ends up being slow we can propose a change there.

cwhelan · 2021-06-17T13:48:48Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapser.java

+        }
+    }
+
+    /**


A high-level description would be useful here.

cwhelan · 2021-06-17T13:51:28Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapser.java

+        final List<Integer> startPositions = items.stream().map(SVCallRecord::getPositionA).sorted().collect(Collectors.toList());
+        final List<Integer> endPositions = items.stream().map(SVCallRecord::getPositionB).sorted().collect(Collectors.toList());
+        //use the mid value of the sorted list so the start and end represent real breakpoint observations
+        final int medianStart = startPositions.get(startPositions.size() / 2);


Could use MathUtils.median

Done. Also replaced mean calculations below with (int)Math.round(MathUtils.sum(startPositions) / (double) startPositions.length); since it's now an array type.

cwhelan · 2021-06-17T13:55:08Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapser.java

+     * @param items
+     * @return (key, value) entry of (start, end)
+     */
+    protected Map.Entry<Integer,Integer> collapseInterval(final Collection<SVCallRecord> items) {


It might be clearer to use a org.apache.commons.lang3.tuple.Pair rather than a Map.Entry

Ah yes that is much better!

cwhelan · 2021-06-17T14:01:25Z

src/test/java/org/broadinstitute/hellbender/tools/sv/SVCallRecordUnitTest.java

+                {SVTestUtils.newCallRecordWithAlgorithms(Collections.emptyList()), false},
+                {SVTestUtils.newCallRecordWithAlgorithms(Collections.singletonList(GATKSVVCFConstants.DEPTH_ALGORITHM)), true},
+                {SVTestUtils.newCallRecordWithAlgorithms(Collections.singletonList("pesr")), false},
+                {SVTestUtils.newCallRecordWithAlgorithms(Lists.newArrayList(GATKSVVCFConstants.DEPTH_ALGORITHM, "pesr")), false}


MIght as well add "pesr" to GATKSVVCFConstants too

We usually use the caller name eg "manta", "melt", etc. for pesr and should allow for any caller name (that is not "depth"), whereas "depth" has a special meaning. Here "pesr" is just for testing.

mwalker174

Thanks @ldgauthier and @cwhelan for the reviews! I believe I've responded to all of them. There's been some major refactoring of the clustering classes as suggested by @cwhelan to make the clustering engine more flexible. I've also made some changes to the SVCallRecord class - in particular the ability to have null length/strand values.

mwalker174 · 2021-06-30T20:32:01Z

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java

+    public static final String CLUSTER_MEMBER_IDS_KEY = "MEMBERS";
+
+    // evidence metrics
+    public static final String COPY_NUMBER_LOG_POSTERIORS_KEY = "CNLP";


I tried to cherry-pick all the changes for this file, but these are actually for a later PR. I'll remove them for now.

mwalker174 · 2021-06-30T20:32:21Z

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java

+
+    // evidence metrics
+    public static final String COPY_NUMBER_LOG_POSTERIORS_KEY = "CNLP";
+    public static final String NEUTRAL_COPY_NUMBER_KEY = "NCN";


It's for ploidy - future PR.

src/main/java/org/broadinstitute/hellbender/tools/sv/SVLocatable.java

mwalker174 · 2021-06-30T20:40:40Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/BinnedCNVDefragmenter.java

+
+/**
+ * CNV defragmenter for when the intervals used for coverage collection are available. The intervals are used to adjust
+ * variant padding by rounding to the nearest bin boundary.


Is there something missing from getPaddedRecordInterval() below?

mwalker174 · 2021-06-30T21:51:47Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/BinnedCNVDefragmenter.java

+ */
+public class BinnedCNVDefragmenter extends CNVDefragmenter {
+
+    protected final TreeMap<GenomeLoc, Integer> genomicToBinMap;


Did you consider using SVIntervalTree? It could be faster because it's contig-aware and was written by @tedsharpe . The downside is converting from GenomeLoc to SVInterval which involves converting contig IDs to ints.

I'm hesitant to start mixing in SVIntervals mainly because of the confusion it can create with differences in coordinates.

Also I might make the name of this map a little more explicit, ie genomicCoordinatesToBinIndexMap

Done

src/test/java/org/broadinstitute/hellbender/tools/sv/SVTestUtils.java

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

mwalker174 · 2021-07-13T01:30:47Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

+                        final int positionB,
+                        final boolean strandB,
+                        final StructuralVariantType type,
+                        final int length,


I've made length and the strands into Integer and Boolean types, respectively, so that they can take on null values when appropriate (e.g. length for BNDs or strands for CNVs). This required quite a bit of refactoring and fixing of tests but I seem to have it working. There may be some performance hit from boxing/unboxing these values everywhere but I do agree that passing around -1 everywhere (and having non-sensical strand values) is problematic in some ways.

It is a little tricky to have constructors that allow for length/strand parameters that are potentially inconsistent, so I have checks to ensure the values are consistent and also sanitize attributes that could contain redundant information.

mwalker174 · 2021-07-13T01:31:28Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

-    @Override
-    public int getStart() {
-        return start;
+    public boolean isCNV() {


mwalker174 · 2021-07-13T01:32:18Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecordUtils.java

+        }
+        Utils.validateArg(record.isIntrachromosomal(), "Inversion is not intrachromosomal");
+        final SVCallRecord positiveBreakend = new SVCallRecord(record.getId(), record.getContigA(),
+                record.getPositionA(), true, record.getContigB(), record.getPositionB(), true, StructuralVariantType.BND, -1,


Fixed this - see comment in SVCallRecord.java.

cwhelan

This looks good to me now; sorry for the delay in re-review. I only have a couple of very minor stylistic comments; feel free to merge after deciding whether or not to make changes based on them.

cwhelan · 2021-08-06T15:38:41Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecordUtils.java

+        return newAttributes;
+    }
+
+    private static Integer getLength(final VariantContext variant, final StructuralVariantType type) {


I think it's a little confusing to have getLength and getStrands return null if the values are inferrable from the SV type -- in my view that'd be surprising behavior if I was reading through the SVCallRecord.create() method. I'd recommend putting the SVCallRecord can resolve these clauses of these if statements back into SVCallRecord.create() and having these methods just have the else clauses -- parsing the attribute values, etc.

cwhelan · 2021-08-13T16:00:01Z

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/CanonicalSVCollapser.java

+                if (i < tokens.length && subtype.equals(tokens[i])) {
+                    alleleSize = i + 1;
+                } else {
+                    break outerloop;


Labeled breaks always look a little weird to me (probably because I've worked on other codebases where they were explicitly against the style guide) and we don't have many of them in the GATK codebase. You could replace this with return Arrays.copyOf(firstTokens, alleleSize); since there's nothing else going on after the method, or refactor the loop labeled outerloop into a calculateAlleleSize method and just return from from it. But I guess they're not against current style guidelines and my aversion to them might be just a personal quirk, so feel free to keep this if you like it better.

Done, just replaced with return Arrays.copyOf(firstTokens, alleleSize);

ldgauthier

I still have to go over the tests, but here are the code comments. Mostly requests to rename, but there was one concern I had about dupes with no-call GTs.

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

ldgauthier · 2021-08-24T19:33:26Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

+        final int copyNumber = VariantContextGetters.getAttributeAsInt(genotype, COPY_NUMBER_FORMAT, 0);
+        if (altAllele.equals(Allele.SV_SIMPLE_DEL)) {
+            return copyNumber < ploidy;
+        } else {


I see that there are some nice checks in the calling method, but can you put javadoc to be explicit about the parameters?

Not sure if this is relevant, but I clearly documented the new isCarrier() method.

ldgauthier · 2021-08-24T19:36:57Z

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

-        for (final Genotype g : variant.getGenotypes()) {
-            if (g.isHomRef() || (g.isNoCall() && !g.hasExtendedAttribute(GATKSVVCFConstants.COPY_NUMBER_FORMAT))) {
-                continue;
+    private static Pair<Boolean, Boolean> inferStrands(final StructuralVariantType type,


I'll let Chris double check this logic.

src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/SVClusterEngineFactory.java

src/main/java/org/broadinstitute/hellbender/tools/sv/cluster/CanonicalSVCollapser.java

src/test/java/org/broadinstitute/hellbender/tools/sv/SVTestUtils.java

ldgauthier · 2021-08-25T18:43:25Z

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/CNVDefragmenterTest.java

+        Assert.assertFalse(defragmenter.areClusterable(deletionRAR, deletionAAR));
+        Assert.assertFalse(defragmenter.areClusterable(deletionRAR, deletionAAA));
+
+        Assert.assertTrue(defragmenter.areClusterable(deletionARR1, deletionARR1Copy));


Ooh good call. I have a branch where I ran JointGermlineCNVSegmentation on two multi-sample input VCFs, but never merged -- is that possible with this code now? If not then maybe I want to rebase before you merge...

ldgauthier

I think you might be missing a few test cases, but otherwise just tidying up that you can choose to do or not.

src/test/java/org/broadinstitute/hellbender/tools/sv/SVCallRecordUtilsUnitTest.java

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapserTest.java

ldgauthier · 2021-09-02T20:05:17Z

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapserTest.java

+                                Collections.singletonList(Allele.REF_N)
+                        ),
+                        Collections.singletonList(Allele.SV_SIMPLE_DEL),
+                        Lists.newArrayList(Allele.REF_N, Allele.SV_SIMPLE_DEL)


Can you add a comment about preferring higher ploidy?

ldgauthier · 2021-09-02T20:06:59Z

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapserTest.java

+                        Lists.newArrayList(Allele.SV_SIMPLE_DEL, Allele.SV_SIMPLE_DUP),
+                        Lists.newArrayList(Allele.SV_SIMPLE_DUP, Allele.SV_SIMPLE_DEL)
+                },
+        };


The javadoc for the method says "potentially augmented with additional ref alleles to match the ploidy" -- is there a test for that? Like DEL, REF/REF ?

Good catch, added a test for that and additional CNV cases such as haploid DUPs where we can infer GT exactly.

ldgauthier · 2021-09-03T14:31:54Z

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapserTest.java

+                {
+                        new int[]{-1},
+                        new String[]{"chr2"},
+                        new StructuralVariantType[]{StructuralVariantType.BND},


Can you add another BND here or another case with multiple?

ldgauthier · 2021-09-03T14:38:44Z

src/test/java/org/broadinstitute/hellbender/tools/sv/cluster/SVCollapserTest.java

+    public Object[][] collapseAlgorithmsTestData() {
+        return new Object[][]{
+                {Collections.singletonList(Collections.singletonList(GATKSVVCFConstants.DEPTH_ALGORITHM)), Collections.singletonList(GATKSVVCFConstants.DEPTH_ALGORITHM)},
+                {Collections.singletonList(Collections.singletonList("pesr")), Collections.singletonList("pesr")},


Are there cases where events can be supported by depth and pesr evidence?

Yes, added those cases now

mwalker174 · 2021-10-05T19:44:58Z

Back to you @ldgauthier

ldgauthier · 2021-10-06T16:33:24Z

This must be the most discussed tool in the GATK. I am relieved to close the book on this PR. Very nice work @mwalker174 -- thanks for all the refactoring!

mwalker174 requested review from ldgauthier and cwhelan May 6, 2021 21:28

ldgauthier reviewed May 13, 2021

View reviewed changes

src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java Outdated Show resolved Hide resolved

ldgauthier requested changes Jun 9, 2021

View reviewed changes

cwhelan requested changes Jun 17, 2021

View reviewed changes

mwalker174 force-pushed the mw_sv_cluster_engine branch from 75dd534 to fa907fd Compare July 13, 2021 01:33

mwalker174 commented Jul 13, 2021

View reviewed changes

cwhelan approved these changes Aug 13, 2021

View reviewed changes

mwalker174 added 7 commits August 18, 2021 15:34

Update SV clustering backend

910ca37

Fix compiler warning

f01628a

Start addressing reviewer comments

480df65

Fix SVTestUtils makeRecord

6edd62c

Set ref allele to start position base; fix tests

f6c9f96

Address remaining comments

41c691d

Fix compiler warnings and address final comments

a65420d

mwalker174 force-pushed the mw_sv_cluster_engine branch from fa907fd to a65420d Compare August 18, 2021 19:48

ldgauthier requested changes Aug 25, 2021

View reviewed changes

ldgauthier requested changes Sep 3, 2021

View reviewed changes

mwalker174 added 3 commits September 13, 2021 15:36

Start refactoring

b95753f

Finish refactoring; fix tests

f2fce13

Address comments

2453a21

ldgauthier approved these changes Oct 6, 2021

View reviewed changes

mwalker174 merged commit 95852a1 into master Oct 6, 2021

mwalker174 deleted the mw_sv_cluster_engine branch October 6, 2021 18:58

gileshall pushed a commit that referenced this pull request Oct 6, 2021

Update SV clustering classes (#7243)

f76fc90

mwalker174 mentioned this pull request Nov 1, 2021

SVCluster tool #7541

Merged

+                      }
+                  }
+                  /**

Update SV clustering classes #7243

Update SV clustering classes #7243

Conversation

mwalker174 commented May 6, 2021

ldgauthier 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

cwhelan 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

mwalker174 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

cwhelan 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

ldgauthier 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

ldgauthier 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

mwalker174 commented Oct 5, 2021

ldgauthier commented Oct 6, 2021