Bg pileup detection #7432

bhanugandham · 2021-08-20T21:34:24Z

This is the beta version of the pileup based haplotype generation and variant detection code.

bhanugandham · 2021-08-20T21:35:12Z

@jamesemery as promised here is the draft PR of the pileup code that works for M2 and HC

jamesemery

I have given this code a first pass. Overall it looks good but there are some thing that need to be cleaned up/leftover comments/little things to deal with. I think my biggest comments are that there are a few methods that really need some unit tests about their behavior and I would be happy to talk over how to go about writing those. Given that there is a lot of interest in this code i think it will end up paying off investing the time early to write tests for some of the tricky methods you added (namely the filtering of the discovered haplotypes and some of the code you are using to add haplotypes).

I also think there should be more parameters exposed for this code in case other projects in the future want to try custom configurations of the pileup detection that might be more appropriate for their use case.

jamesemery · 2021-08-23T17:25:47Z

src/main/java/META-INF/MANIFEST.MF

@@ -0,0 +1,3 @@
+Manifest-Version: 1.0
+Main-Class: org.broadinstitute.hellbender.Main


This looks like an accidental addition.

will remove

jamesemery · 2021-09-01T16:55:46Z

...roadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtilsUnitTest.java

+//    @Test
+//    public void testAddGivenAlleles() {
+//        final int assemblyRegionStart = 1;
+//        final int maxMnpDistance = 0;


Reenable these tests

jamesemery · 2021-09-01T17:04:09Z

src/main/java/org/broadinstitute/hellbender/engine/AssemblyRegion.java

@@ -114,6 +123,13 @@ public AssemblyRegion(final SimpleInterval activeSpan, final int padding, final
        this(activeSpan, true, padding, header);
    }

+    public List<AlignmentData> getAlignmentData() {


Comment these

jamesemery · 2021-09-01T17:11:40Z

src/main/java/org/broadinstitute/hellbender/engine/AssemblyRegion.java

+        addToReadCollection(read, reads);
+    }
+
+    private void addToReadCollection(final GATKRead read, final List<GATKRead> collection) {


"addToReadCollection" -> "addToReadCollectionAndValidate"

jamesemery · 2021-09-01T17:13:35Z

src/main/java/org/broadinstitute/hellbender/engine/AssemblyRegionIterator.java

@@ -113,12 +117,21 @@ public AssemblyRegion next() {
        return toReturn;
    }

+//    private void addAlignmentData(AlignmentData alignmentData){


commented code

jamesemery · 2021-09-01T18:59:53Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+                }
+
+                // now look for indels
+                if (element.isBeforeInsertion()) {


There should be an argument in the read threading assembler to enable/disable the indel checks at this stage.

@jamesemery unclear what you mean here.

There are some other projects that are interested in this code. They are not necessarily interested in the pileup indel calling code right now (since its complicated and they might have their own versions of this code) Ideally there should be a toggle that disables all of the indel pileup code. (or that disables putting them into the haplotypes).

jamesemery · 2021-09-01T19:00:06Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+            // TODO: AH & BG add an option to deal with multiple alt alleles
+            Optional<Map.Entry<Byte, Integer>> maxAlt = altCounts.entrySet().stream().max(Comparator.comparingInt(Map.Entry::getValue));
+            if (maxAlt.isPresent() && ((float)maxAlt.get().getValue() / (float)numOfBases) > 0.10 && numOfBases >= 5 ) {
+//            if (maxAlt.isPresent() && ((float)maxAlt.get().getValue() / (float)numOfBases) > 0.1 && maxAlt.get().getValue() > 10 ) {


more commented code.

jamesemery · 2021-09-01T19:01:07Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+            alleles.add(Allele.create(referenceContext.getBase(), true));
+            // TODO: AH & BG add an option to deal with multiple alt alleles
+            Optional<Map.Entry<Byte, Integer>> maxAlt = altCounts.entrySet().stream().max(Comparator.comparingInt(Map.Entry::getValue));
+            if (maxAlt.isPresent() && ((float)maxAlt.get().getValue() / (float)numOfBases) > 0.10 && numOfBases >= 5 ) {


I would pull this out as an explicit method filterAlleles() or something along those lines. I would also make static variables for these magic numbers.

I also think these should ulitmately be exposed for configuration down the road... There are a few ways to think about doing it but we should probably discuss how that works here.

jamesemery · 2021-09-01T19:02:32Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+            indelAlleles.add(Allele.create(referenceContext.getBase(), true));
+            Optional<Map.Entry<String, Integer>> maxIns = insertionCounts.entrySet().stream().max(Comparator.comparingInt(Map.Entry::getValue));
+            if (maxIns.isPresent() && ((float)maxIns.get().getValue() / (float)numOfBases) > 0.50 && numOfBases >= 5 ) {
+                indelAlleles.add(Allele.create(String.format("%c", referenceContext.getBase()) + maxIns.get().getKey()));


drop the String.Format() call here since this might get triggered a lot.

src/main/java/org/broadinstitute/hellbender/engine/AssemblyRegionIterator.java

yfarjoun · 2021-11-29T00:01:38Z

src/main/java/org/broadinstitute/hellbender/engine/AssemblyRegionIterator.java

@@ -86,6 +88,8 @@ public AssemblyRegionIterator(final MultiIntervalShard<GATKRead> readShard,
        this.readCachingIterator = new ReadCachingIterator(readShard.iterator());
        this.readCache = new ArrayDeque<>();
        this.activityProfile = new BandPassActivityProfile(assemblyRegionArgs.maxProbPropagationDistance, assemblyRegionArgs.activeProbThreshold, BandPassActivityProfile.MAX_FILTER_SIZE, BandPassActivityProfile.DEFAULT_SIGMA, readHeader);
+// TODO: AH & BG handle changing contig


@ahaessly can you elaborate what was missing WRT this TODO?

Perhaps a simple test case that would cover the circumstance you were thinking of?

This was something that james warned us about. I think the shards can be from different contigs so I think we needed to add code that checked for that case.

src/main/java/org/broadinstitute/hellbender/utils/haplotype/Haplotype.java

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

...oadinstitute/hellbender/tools/walkers/haplotypecaller/PileupDetectionArgumentCollection.java

yfarjoun · 2022-01-09T23:07:52Z

...oadinstitute/hellbender/tools/walkers/haplotypecaller/PileupDetectionArgumentCollection.java

+     * Percentage of reads required to supprot the alt for a variant to be considered
+     */
+    @Advanced
+    @Argument(fullName= PILEUP_DETECTION_SNP_THRESHOLD, doc = "Percentage of alt supporting reads in order to consider alt SNP", optional = true)


I would rename argument PILEUP_SNP_DETECTION_RATIO (also below for INDEL) also, in doc make it clear that this is regarding pileup mode

It would be good to note in the doc= of all these arguments that it pertains to the PILEUP detection only. otherwise the only hint to that it in the name of the argument itself.

I've changed all of the arguments to have a "Pileup Detection: " argument. Hopefully that will make it a little bit clearer? I've hidden these arguments anyway so maybe its all moot...

...oadinstitute/hellbender/tools/walkers/haplotypecaller/PileupDetectionArgumentCollection.java

yfarjoun · 2022-01-09T23:25:27Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+     * that are visible in the pileups but might be dropped from assembly for any number of reasons. The basic algorithm works
+     * as follows:
+     *  - iterate over every pileup and count alt bases
+     *      - (unfinished) detect insertions overlapping this site (CURRENTLY ONLY WORKS FOR INSERTIONS)


it's unclear what is unfinished...do you want this to work for indels, but it only works for insertions?

I consider the implementation of indels to be relatively untested compared to the SNP code. I'm going to change this to beta.

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

yfarjoun · 2022-01-09T23:36:37Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+        if (args.badReadProperPair && !read.isProperlyPaired()) {
+            return true;
+        }
+        if (args.badReadSecondary && (read.isSecondaryAlignment() || read.hasAttribute("SA"))) {


in the doc for badReadSecondary there's no mention of secondary alignments....and yet, here it is...please make sure we are doing it right....

Good catch. The DRAGEN heuristic attempts to capture both (on the assumption that those messy reads are going to cause false positives). I'll update the documentation.

yfarjoun · 2022-01-09T23:36:57Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+        }
+
+
+        //TODO this conversion is really unnecessary...


huh? so why is it done?

Its done because the method SequenceUtil.calculateSamNmTag() works on SAMRecords and its in HTSJDK so I can't easily change it... I didn't want to reimplement it but as you can see i have been forced to thread the header down to this class as a result.

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

yfarjoun · 2022-01-09T23:39:10Z

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

+        }
+
+        //TODO add threshold descibed by illumina about insert size compared to the average
+        if (args.templateLenghtStd > 0 && args.templateLengthMean > 0) {


insert-size is a much better term than template length.....i think we should use that instead.

We are already using samRecordForRead.getInferredInsertSize() or do you mean we should change the title/arguments for this to use "insert size" as a name?

src/main/java/org/broadinstitute/hellbender/utils/pileup/PileupBasedAlleles.java

yfarjoun · 2022-01-09T23:43:33Z

...es/org/broadinstitute/hellbender/tools/haplotypecaller/expected.pileupCallerDRAGEN.gatk4.vcf

+20	10098786	.	C	T	32.64	.	AC=1;AF=0.500;AN=2;BaseQRankSum=0.431;DP=9;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=54.65;MQRankSum=-2.200;QD=3.63;ReadPosRankSum=-1.383;SOR=0.132	GT:AD:DP:GQ:PL	0/1:7,2:9:40:40,0,217
+20	10099220	.	A	G	598.64	.	AC=1;AF=0.500;AN=2;BaseQRankSum=-1.664;DP=47;ExcessHet=3.0103;FS=1.137;MLEAC=1;MLEAF=0.500;MQ=57.74;MQRankSum=-2.596;QD=13.30;ReadPosRankSum=0.629;SOR=0.519	GT:AD:DP:GQ:PL	0/1:25,20:45:99:606,0,873
+20	10099535	.	G	A	1356.64	.	AC=1;AF=0.500;AN=2;BaseQRankSum=5.034;DP=68;ExcessHet=3.0103;FS=2.172;MLEAC=1;MLEAF=0.500;MQ=59.04;MQRankSum=-1.418;QD=20.87;ReadPosRankSum=-0.777;SOR=0.408	GT:AD:DP:GQ:PL	0/1:26,39:65:99:1364,0,686
+20	10099565	.	C	T	1367.64	.	AC=1;AF=0.500;AN=2;BaseQRankSum=4.595;DP=70;ExcessHet=3.0103;FS=8.736;MLEAC=1;MLEAF=0.500;MQ=59.40;MQRankSum=-1.266;QD=19.82;ReadPosRankSum=1.225;SOR=0.191	GT:AD:DP:GQ:PL	0/1:31,38:69:99:1375,0,994


there's a utuilty method that indexes vcfs for testing, perhaps better to use it than to include indexes in the repository?

Looking at that folder it looks like (with a few exceptions that look like they aren't actually used anywhere and should probably be removed themselves) every file also has the indexes checked in. I'm going to call it out of scope for this PR to change the HC tests over to index on the fly.

yfarjoun

Done with this pass. let me know that you see my comments...they were only put on your last commit.

jamesemery · 2022-01-11T16:41:50Z

@yfarjoun rebased (It was a bit messy its quite possible I broke something) and back to you.

yfarjoun

I still see many unresolved reviewer comments from previous passes, but the tests that you added are good.

...oadinstitute/hellbender/tools/walkers/haplotypecaller/PileupDetectionArgumentCollection.java

...institute/hellbender/tools/walkers/haplotypecaller/readthreading/ReadThreadingAssembler.java

.../broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java

jamesemery · 2022-02-04T16:51:28Z

@yfarjoun I think I have responded to everything in there and gotten rid of a bunch of the todos that were left (it turns out a bunch of them were vestigial from development). There is really only one of my comments I was unsure was worth fixing and I will probably ask around engine team for opinions about it.

yfarjoun

a few small nits remain. thanks for moving this along!

yfarjoun · 2022-02-12T22:30:47Z

...roadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtilsUnitTest.java

+        final Haplotype hapD = new Haplotype("ACCTGAA".getBytes());
+        final Haplotype hapF = new Haplotype("GAAGAAG".getBytes()); // testing repeated kmers
+
+        Map<Kmer, Integer> flatSupportAllKmers = new HashMap<Kmer, Integer>() { private static final long serialVersionUID = 0L; {


put the serialVersionUID into a line of its own please (here and below)

Sure... I'm still annoyed that i have to do this but hash-map is serializable so this anonymous subclass generates a compiler warning otherwise...

src/main/java/org/broadinstitute/hellbender/engine/AssemblyRegionIterator.java

...va/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AssemblyBasedCallerUtils.java

yfarjoun · 2022-02-13T02:12:54Z

.../java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngine.java

@@ -771,7 +771,7 @@ private boolean containsCalls(final CalledHaplotypes calledHaplotypes) {
            //TODO - why the activeRegion cannot manage its own one-time finalization and filtering?
            //TODO - perhaps we can remove the last parameter of this method and the three lines bellow?


are these todos still relevant and correct?

no but they are old enough now maybe they are load bearing? There is a bunch of this commentary leftover from the first port of HaplotypeCaller back in the day that can probably go...

...oadinstitute/hellbender/tools/walkers/haplotypecaller/PileupDetectionArgumentCollection.java

…ling tails

…e place it was before despite the GGA code having been moved). And fixed a broken integration test

jamesemery · 2022-03-29T19:03:00Z

@yfarjoun That was a messy rebase. Assuming these tests pass everything should be in order for you to take another look.

gatk-bot · 2022-03-29T20:02:25Z

Travis reported job failures from build 38310
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
variantcalling	openjdk8	38310.4	logs

yfarjoun

just a small nit in one of the comments....thanks!

Its my own review and I have made the changes necessary

…ore assembly that supplements the assembly variants with variants that show up in the pileups. (broadinstitute#7432)

bhanugandham requested review from jamesemery and davidbenjamin August 20, 2021 21:34

droazen assigned jamesemery Aug 26, 2021

jamesemery previously requested changes Sep 1, 2021

View reviewed changes

yfarjoun reviewed Nov 29, 2021

View reviewed changes

yfarjoun force-pushed the bg_pileup_detection branch from 58e59ea to 23f29c8 Compare November 29, 2021 02:25

jamesemery reviewed Dec 6, 2021

View reviewed changes

src/main/java/org/broadinstitute/hellbender/utils/haplotype/Haplotype.java Outdated Show resolved Hide resolved