Implement SVConcordance tool #7977

mwalker174 · 2022-08-05T20:57:14Z

Analogous to GenotypeConcordance, but for SVs.

codecov · 2022-08-05T22:20:43Z

Codecov Report

Merging #7977 (27cf6ef) into master (d1a4f94) will decrease coverage by 0.006%.
The diff coverage is 92.234%.

Additional details and impacted files

@@               Coverage Diff               @@
##              master     #7977       +/-   ##
===============================================
- Coverage     86.658%   86.652%   -0.006%     
- Complexity     39024     39154      +130     
===============================================
  Files           2337      2346        +9     
  Lines         182961    183608      +647     
  Branches       20079     20147       +68     
===============================================
+ Hits          158550    159100      +550     
- Misses         17365     17447       +82     
- Partials        7046      7061       +15

Impacted Files	Coverage Δ
...roadinstitute/hellbender/tools/sv/SVLocatable.java	`25.000% <ø> (ø)`
...titute/hellbender/utils/GenotypeUtilsUnitTest.java	`94.624% <ø> (ø)`
...ellbender/tools/sv/cluster/CanonicalSVLinkage.java	`83.146% <62.963%> (-6.043%)`	⬇️
...ender/utils/variant/GATKSVVariantContextUtils.java	`93.939% <66.667%> (-2.835%)`	⬇️
...stitute/hellbender/tools/sv/SVCallRecordUtils.java	`79.167% <73.810%> (-4.518%)`	⬇️
...der/tools/walkers/sv/SVClusterIntegrationTest.java	`98.974% <77.778%> (-0.522%)`	⬇️
...oadinstitute/hellbender/tools/sv/SVCallRecord.java	`84.348% <79.167%> (+2.144%)`	⬆️
...r/tools/sv/concordance/SVConcordanceAnnotator.java	`83.465% <83.465%> (ø)`
...roadinstitute/hellbender/tools/sv/SVTestUtils.java	`90.274% <84.286%> (-6.426%)`	⬇️
...lbender/tools/sv/cluster/CanonicalSVCollapser.java	`80.635% <86.508%> (-6.836%)`	⬇️
... and 30 more

TedBrookings

This was a challenging review for me: there's lots of detailed logic (if this kind of SV with this property, then do X) that is difficult for me to check. I did try to follow it all though, and I didn't notice anything wrong. In a few places, I indicated that it would be ideal to add comments explaining the code logic. I also found a few places where the code could be altered to achieve better efficiency / conform more to THE JAVA WAY (as I understand it). I did not see anything that absolutely needs to be changed, so I am approving.

TedBrookings · 2022-10-05T15:24:08Z

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

+    /**
+     * Get allele frequencies (AF). Avoid calling this repeatedly on the same object, as computations are performed
+     * on the fly.
+     */


It seems like if you wanted to prevent these frequencies from being recalculated, you could make them an assigned class member and make the routine to compute them a private function called by the constructor (and then have a public getter).

I have no idea why I used a factory here… I've changed this to a simple class like you suggested.

TedBrookings · 2022-10-05T15:48:46Z

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

@@ -215,7 +222,9 @@ private static Pair<Boolean, Boolean> inferStrands(final StructuralVariantType t
            return Pair.of(Boolean.FALSE, Boolean.TRUE);
        } else {
            if (type.equals(StructuralVariantType.INV)) {
-            Utils.validateArg(inputStrandA != null && inputStrandB != null, "Cannot create variant of type " + type + " with null strands");
+                if (inputStrandA == null && inputStrandB == null) {
+                    return Pair.of(Boolean.TRUE, Boolean.TRUE);


Maybe add a comment explaining why it's okay to set the strands to True if they are input as null?

Thanks for catching this. I did some digging and I actually found our cleaned vcfs don't have STRANDS for INV types and I think this is fine. We do keep them there in intermediate vcfs since manta reports them that way, but we only retain ones with all 4 breakends (see manta manual about that). I've modified the check to allow null.

TedBrookings · 2022-10-06T14:10:42Z

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

+     * Asserts presence of {@link GATKSVVCFConstants#EXPECTED_COPY_NUMBER_FORMAT} and
+     * {@link GATKSVVCFConstants#COPY_NUMBER_FORMAT} attributes.
+     */
+    public static void assertHasCopyStateFields(final Genotype g) {


I feel like the word should be written out. In complicated SV code every little bit helps.

Suggested change

public static void assertHasCopyStateFields(final Genotype g) {

public static void assertHasCopyStateFields(final Genotype genotype) {

TedBrookings · 2022-10-06T14:11:15Z

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

+            return record;
+        }
+        final ArrayList<Genotype> newGenotypes = new ArrayList<>(record.getGenotypes().size());
+        for (final Genotype g : record.getGenotypes()) {


Same deal, g -> genotype

TedBrookings · 2022-10-06T14:18:02Z

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

@@ -312,7 +314,7 @@ protected List<Allele> collapseSampleGenotypeAlleles(final Collection<Genotype>
        final List<Allele> alleles;
        if (altAlleles.contains(Allele.SV_SIMPLE_DEL) || altAlleles.contains(Allele.SV_SIMPLE_DUP)) {
            // For CNVs, collapse genotype using copy number
-            alleles = getCNVGenotypeAllelesFromCopyNumber(altAlleles, refAllele, expectedCopyNumber, copyNumber);
+            alleles = getCNVGenotypeAllelesFromCopyNumber(altAlleles, refAllele, expectedCopyNumber, copyNumber, true);


Since sometimes you assume biallelic DUPs and sometimes you don't, it might be nice to include a brief comment explaining the choice in each location.

Biallelic dups are now assumed everywhere for simple DUP types

TedBrookings · 2022-10-06T15:31:46Z

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

+                for (final String sample : samples) {
+                    final Genotype genotypeA = genotypesA.get(sample);
+                    final Genotype genotypeB = genotypesB.get(sample);
+                    final int cnA = genotypeA == null ?


Maybe add a quick comment here explaining why genotypeB having ECN == CN counts as a match if genotypeA is null? I'm guessing it's because null Genotype means that not CNV was found for that sample, but ECN=CN means that it was explicitly called REF in the other?

Yes if a sample doesn't exist in one of the records, this assumes it's copy-neutral. I've added a comment to clarify.

TedBrookings · 2022-10-06T16:14:04Z

src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java

+/**
+ * Efficiently clusters a set of evaluation ("eval") SVs with their closest truth SVs.
+ *
+ * "Closest" is defined in the {@link #getClosestItem} method an selects first on total breakpoint distance


Suggested change

* "Closest" is defined in the {@link #getClosestItem} method an selects first on total breakpoint distance

* "Closest" is defined in the {@link #getClosestItem} method and selects first on total breakpoint distance

TedBrookings · 2022-10-06T16:45:05Z

src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java

+            final List<Long> finalizedRefItems = truthIdToItemMap.entrySet().stream()
+                    .filter(e -> linkage.getMaxClusterableStartingPosition(e.getValue()) < lastItemStart)
+                    .map(Map.Entry::getKey)
+                    .collect(Collectors.toList());
+            finalizedRefItems.forEach(truthIdToItemMap::remove);


You can avoid creating the list and doing a second iteration. Changes to values() are backed by the set, and removeIf does exactly what you want.
https://docs.oracle.com/javase/8/docs/api/java/util/Map.html#values--
https://docs.oracle.com/javase/8/docs/api/java/util/Collection.html#removeIf-java.util.function.Predicate-

Suggested change

final List<Long> finalizedRefItems = truthIdToItemMap.entrySet().stream()

.filter(e -> linkage.getMaxClusterableStartingPosition(e.getValue()) < lastItemStart)

.map(Map.Entry::getKey)

.collect(Collectors.toList());

finalizedRefItems.forEach(truthIdToItemMap::remove);

truthIdToItemMap.values()

.removeIf(cluster -> linkage.getMaxClusterableStartingPosition(cluster) < lastItemStart);

Nice! Never knew that values() could be modified like that

TedBrookings · 2022-10-06T19:04:26Z

src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java

+        final List<Map.Entry<Long, SVCallRecord>> linkedItems = candidates.stream()
+                .filter(other -> linkage.areClusterable(evalRecord, other.getValue()))
+                .collect(Collectors.toList());
+        final int[] distances = linkedItems.stream()
+                .mapToInt(e -> totalDistance(evalRecord, e.getValue()))
+                .toArray();
+        if (distances.length == 0) {
+            return null;
+        }
+        final int minDistance = MathUtils.arrayMin(distances);
+        int numMin = 0;
+        int minDistIndex = 0;
+        for (int i = 0; i < distances.length; i++) {
+            if (distances[i] == minDistance) {
+                numMin++;
+                minDistIndex = i;
+            }
+        }
+        if (numMin == 1) {
+            return linkedItems.get(minDistIndex);
+        } else {
+            return getClosestItemWithTiebreakers(evalRecord, linkedItems);
+        }


Not sure if performance / etc is worth the time of rewriting, but you can construct a Comparator that breaks ties by chaining with .thenComparing(). You could then have a single stream terminated with min() acting on that Comparator. This way you could have arbitrary comparisons without passing around intermediate lists that get winnowed down. If the comparator was constructed from a function (with input evalRecord), it could be used in update() below as well, decreasing the potential for bugs if the two code paths get out of sync.

Yep that's a much better way to do it. That cut down the code a lot

TedBrookings · 2022-10-06T20:39:02Z

src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/SVConcordanceAnnotator.java

+        } else if (g.isMixed()) {
+            return GenotypeConcordanceStates.CallState.IS_MIXED;
+        } else {
+            throw new IllegalArgumentException("Could not determine truth state for genotype: " + g);


Presumably a typo. Although I wonder if you could combine getTruthState and getEvalState by passing a boolean allowNull

Suggested change

throw new IllegalArgumentException("Could not determine truth state for genotype: " + g);

throw new IllegalArgumentException("Could not determine eval state for genotype: " + g);

…upport

…ations

mwalker174 mentioned this pull request Sep 21, 2022

Genotype concordance workflows broadinstitute/gatk-sv#414

Merged

mwalker174 requested a review from TedBrookings October 4, 2022 16:15

TedBrookings approved these changes Oct 7, 2022

View reviewed changes

mwalker174 force-pushed the mw_sv_concordance branch from c0241d1 to c6626ce Compare January 21, 2023 23:13

mwalker174 added 3 commits February 1, 2023 11:36

Implement SVConcordance tool

77fa2f8

Address reviewer comments; greatly simplify clustering; add ctx/cpx s…

52017bb

…upport

Treat missing truth samples as no-call; remove no-call from AF calcul…

69ec8ad

…ations

mwalker174 force-pushed the mw_sv_concordance branch from eaacae4 to 69ec8ad Compare February 1, 2023 16:36

Don't write CN fields if not in the header

27cf6ef

mwalker174 merged commit d3b2c5c into master Feb 9, 2023

mwalker174 deleted the mw_sv_concordance branch February 9, 2023 03:15

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement SVConcordance tool #7977

Implement SVConcordance tool #7977

mwalker174 commented Aug 5, 2022

codecov bot commented Aug 5, 2022 •

edited

Loading

TedBrookings left a comment

TedBrookings Oct 5, 2022

mwalker174 Dec 19, 2022

TedBrookings Oct 5, 2022

mwalker174 Dec 19, 2022 •

edited

Loading

TedBrookings Oct 6, 2022

mwalker174 Dec 19, 2022

TedBrookings Oct 6, 2022

TedBrookings Oct 6, 2022

mwalker174 Dec 20, 2022

TedBrookings Oct 6, 2022

mwalker174 Dec 20, 2022

TedBrookings Oct 6, 2022

TedBrookings Oct 6, 2022

mwalker174 Dec 20, 2022

TedBrookings Oct 6, 2022

mwalker174 Dec 20, 2022

TedBrookings Oct 6, 2022

	public static void assertHasCopyStateFields(final Genotype g) {
	public static void assertHasCopyStateFields(final Genotype genotype) {

	* "Closest" is defined in the {@link #getClosestItem} method an selects first on total breakpoint distance
	* "Closest" is defined in the {@link #getClosestItem} method and selects first on total breakpoint distance

	throw new IllegalArgumentException("Could not determine truth state for genotype: " + g);
	throw new IllegalArgumentException("Could not determine eval state for genotype: " + g);

Implement SVConcordance tool #7977

Implement SVConcordance tool #7977

Conversation

mwalker174 commented Aug 5, 2022

codecov bot commented Aug 5, 2022 • edited Loading

Codecov Report

TedBrookings 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

mwalker174 Dec 19, 2022 • edited Loading

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 Aug 5, 2022 •

edited

Loading

mwalker174 Dec 19, 2022 •

edited

Loading