fasta gz support in gatk #5120

lbergelson · 2018-08-16T15:02:17Z

this is a bit of a wip but probably good enough to take a look at for sanity.

codecov-io · 2018-08-22T20:05:13Z

Codecov Report

Merging #5120 into master will increase coverage by 0.002%.
The diff coverage is 78.713%.

@@               Coverage Diff               @@
##              master     #5120       +/-   ##
===============================================
+ Coverage     86.797%   86.799%   +0.002%     
- Complexity     29706     29711        +5     
===============================================
  Files           1820      1820               
  Lines         137384    137406       +22     
  Branches       15151     15160        +9     
===============================================
+ Hits          119245    119267       +22     
+ Misses         12651     12635       -16     
- Partials        5488      5504       +16

Impacted Files	Coverage Δ	Complexity Δ
...oadinstitute/hellbender/engine/AssemblyRegion.java	`84.496% <ø> (ø)`	`48 <0> (ø)`	⬇️
...er/tools/walkers/rnaseq/OverhangFixingManager.java	`94.079% <ø> (ø)`	`67 <0> (ø)`	⬇️
.../walkers/rnaseq/OverhangFixingManagerUnitTest.java	`98.83% <ø> (ø)`	`31 <0> (ø)`	⬇️
...plotypecaller/ReadLikelihoodCalculationEngine.java	`100% <ø> (ø)`	`0 <0> (ø)`	⬇️
...haplotypecaller/HaplotypeCallerEngineUnitTest.java	`96.296% <ø> (ø)`	`6 <0> (ø)`	⬇️
...institute/hellbender/exceptions/UserException.java	`67.606% <0%> (-2.467%)`	`3 <0> (ø)`
...lbender/engine/AssemblyRegionIteratorUnitTest.java	`84.444% <100%> (+0.175%)`	`11 <6> (ø)`	⬇️
...llbender/engine/filters/VariantFilterUnitTest.java	`90.244% <100%> (ø)`	`6 <0> (ø)`	⬇️
...tute/hellbender/utils/GenomeLocParserUnitTest.java	`88.435% <100%> (-0.263%)`	`46 <2> (ø)`
.../readthreading/ReadThreadingAssemblerUnitTest.java	`98.718% <100%> (+0.011%)`	`38 <1> (+1)`	⬆️
... and 19 more

droazen

Review complete, back to @lbergelson. Found one bug plus a couple of other issues, but looks good overall.

droazen · 2018-08-28T20:30:48Z

src/main/java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFile.java

@@ -94,26 +77,21 @@ public CachingIndexedFastaSequenceFile(final Path fasta, final FastaSequenceInde
     * @param cacheSize the size of the cache to use in this CachingIndexedFastaReader, must be >= 0
     * @param preserveCase If true, we will keep the case of the underlying bases in the FASTA, otherwise everything is converted to upper case


Needs javadoc for the preserveIUPAC argument.

droazen · 2018-08-28T20:36:05Z

src/main/java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFile.java

-        this(fasta, index, cacheSize, false, false);
+    private static ReferenceSequenceFile requireIndex(Path fasta, ReferenceSequenceFile referenceSequenceFile) {
+        if( !referenceSequenceFile.isIndexed()) {
+            throw new UserException("Couldn't create indexed fasta: " + fasta.toAbsolutePath());


"Fasta index not found, but required" would be a better error message, I think.

The error reporting here is kind of bad, there's no easy way to tell what the problem was. I added a gross block of error diagnosing code, let me know what you think...

Actually, I punted here and moved the error checking all to be in the initial checks. This throws gatk error now if it ends up here without being indexed since we should be preventing that up front.

droazen · 2018-08-28T20:37:48Z

src/main/java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFile.java

+        if( !referenceSequenceFile.isIndexed()) {
+            throw new UserException("Couldn't create indexed fasta: " + fasta.toAbsolutePath());
+        }
+        return referenceSequenceFile;
    }

    /**


Why keep any of these constructors public? Why not force all clients to go through the checkAndCreate() factory methods?

I rewrote it to be only constructors and to make them all perform the necessary safety checks. It's much less confusing now I think.

droazen · 2018-08-28T20:39:39Z

src/main/java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFile.java

@@ -136,7 +114,7 @@ public CachingIndexedFastaSequenceFile(final Path fasta) throws FileNotFoundExce
     * @param fasta The file to open.
     * @param preserveCase If true, we will keep the case of the underlying bases in the FASTA, otherwise everything is converted to upper case


Constructor docs should note that IUPAC bases will not be preserved.

droazen · 2018-08-28T20:41:55Z

src/main/java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFile.java

@@ -221,7 +194,7 @@ public static CachingIndexedFastaSequenceFile checkAndCreate(final Path fastaPat
     * @param fasta The file to open.
     * @param cacheSize the size of the cache to use in this CachingIndexedFastaReader, must be >= 0
     */
-    public CachingIndexedFastaSequenceFile(final Path fasta, final long cacheSize ) throws FileNotFoundException {
+    public CachingIndexedFastaSequenceFile(final Path fasta, final long cacheSize ) {


Should be encapsulated, clients should go through checkAndCreate() instead.

Also add overload of checkAndCreate() that lets the client specify the cache size, and any other overloads necessary to replicate the functionality of the constructors.

I did this the other way around, getting rid of check and create and making all the constructors safe.

removed the unused constructor to set cache size, but exposing the full argument constructor as public so people can do what they want

droazen · 2018-08-28T21:41:56Z

.../java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFileUnitTest.java

+            // times.
+            final int minExpectedHits = (int) Math.floor(
+                    (Math.min(caching.getCacheSize(), contig.getSequenceLength()) - querySize * 2.0) / STEP_SIZE);
+            caching.printEfficiency(Level.WARN);


Again, WARN seems inappropriate here.

droazen · 2018-08-28T21:47:33Z

.../java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFileUnitTest.java

+            int middleStart = (contig.getSequenceLength() - querySize) / 2;
+            int middleStop = middleStart + querySize;
+
+            logger.warn(String.format(


Here too I question the use of warn()

droazen · 2018-08-28T21:50:39Z

.../java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFileUnitTest.java

+                if (stop <= contig.getSequenceLength()) {
+                    ReferenceSequence grabMiddle = caching.getSubsequenceAt(contig.getSequenceName(), middleStart,
+                                                                            middleStop);
+                    ReferenceSequence cachedVal = caching.getSubsequenceAt(contig.getSequenceName(), start, stop);


Should this test be asserting that you get a high cache hit rate, given that you query a nearby interval first?

droazen · 2018-08-28T21:54:09Z

.../java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFileUnitTest.java

-    private int testCases(final IndexedFastaSequenceFile original,
-                          final IndexedFastaSequenceFile casePreserving,
-                          final IndexedFastaSequenceFile allUpper,
+    private int testCases(final ReferenceSequenceFile original,


This method name is terrible...how about mixedCasesTestHelper() or something that ties it to the testMixedCasesInExample() test above?

droazen · 2018-08-28T21:56:21Z

.../java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFileUnitTest.java

-                          final IndexedFastaSequenceFile allUpper,
+    private int testCases(final ReferenceSequenceFile original,
+                          final ReferenceSequenceFile casePreserving,
+                          final ReferenceSequenceFile allUpper,
                          final String contig, final int start, final int stop ) {
        final String orig = fetchBaseString(original, contig, start, stop);
        final String keptCase = fetchBaseString(casePreserving, contig, start, stop);


Um, why does this test call toUpperCase() explicitly on allUpper (see line below this one)? Shouldn't it be testing that allUpper returns uppercase bases on its own?

yeah, that's really weird and wrong. Removed that and the test is fine so I guess it was a no-op.

sooheelee · 2018-09-07T22:30:55Z

Looking forward to this feature!

droazen · 2018-09-12T15:41:47Z

@lbergelson Ping me once this is ready for a final pass

lbergelson · 2018-09-12T19:51:35Z

@droazen Could you take another look at this? I had to hit a lot of files to do the constructor update and close the readers where they weren't being closed already.

droazen

Final review complete, back to @lbergelson. Merge after addressing remaining comments.

src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCaller.java

droazen · 2018-09-14T18:56:53Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCaller.java

@@ -257,5 +250,13 @@ public void closeTool() {
        if ( hcEngine != null ) {
            hcEngine.shutdown();
        }
+
+        if ( referenceReader != null){
+            try {


If you take my suggestion above, you can get rid of this try/catch.

Actually, it makes more sense to put the close() call in HaplotypeCallerEngine.shutdown(), which is already called from closeTool(), and get rid of the tool-level referenceReader field completely. This would also be more consistent with Mutect2

It seems wrong to have HaplotypeCallerEngine close the reader since it's passed in externally. It should be the responsibility of the caller to close it. In HaplotypeCallerSpark we currently broadcast an instance of ReferenceMultiSource which is reused between calls and should not be closed.

Yes, but it's gross having the reader as a new field in the tool itself, and inconsistent with existing patterns. I think it's better to transfer responsibility for the reference reader to the HaplotypeCallerEngine at construction time, and let the shutdown() call there govern when it gets closed. The Spark version can always refrain from calling shutdown() until it's ready to do so, can't it?

I guess so. It looks like it never bothers calling shutdown anyway... so it won't be effected... Maybe we should fix that also.

Actually, on further thought, I don't think the spark version will ever be able to properly clean up if we make this change, not without changing how it uses readers at least. Since it should be calling shutdown whenever it's done with a haplotypecaller engine, but the readers persist through multiple engines.

I changed it but I'm not sure if it's ideal.

src/main/java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFile.java

src/test/java/org/broadinstitute/hellbender/engine/AssemblyRegionUnitTest.java

...g/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngineUnitTest.java

src/test/java/org/broadinstitute/hellbender/utils/read/ReadUtilsUnitTest.java

.../java/org/broadinstitute/hellbender/utils/fasta/CachingIndexedFastaSequenceFileUnitTest.java

lbergelson · 2018-09-20T00:41:34Z

I think I'm cursed. This is failing the oracle build again.

* GATK can now read fasta files that are bgzipped This requires both a .fai and .gzi index A gzi index can be created using the bgzip program the same .fai can be used for the zipped and unzipped version of a file * replacing mixed use of CachingIndexFastaFile constructors and static factory methods with a single set of constructors * closing instances of CachingIndexFastaFile

lbergelson · 2018-09-20T02:22:32Z

It was cursed with a bad cache.

responded to comments

lbergelson requested a review from droazen August 16, 2018 15:02

lbergelson assigned droazen Aug 16, 2018

lbergelson force-pushed the lb_fasta_gz branch 2 times, most recently from c6bf892 to 2870e36 Compare August 22, 2018 19:01

droazen requested changes Aug 28, 2018

View reviewed changes

droazen assigned lbergelson and unassigned droazen Aug 28, 2018

lbergelson force-pushed the lb_fasta_gz branch from 2870e36 to 66a5f6c Compare September 12, 2018 19:50

droazen assigned droazen and unassigned lbergelson Sep 13, 2018

droazen previously requested changes Sep 14, 2018

View reviewed changes

droazen assigned lbergelson and unassigned droazen Sep 14, 2018

lbergelson force-pushed the lb_fasta_gz branch from 66a5f6c to e1ebd3a Compare September 19, 2018 18:22

lbergelson added 3 commits September 19, 2018 20:41

more comments

5510225

removing kaboom

21203af

lbergelson force-pushed the lb_fasta_gz branch from e66322a to 21203af Compare September 20, 2018 00:42

lbergelson merged commit 6e352bb into master Sep 20, 2018

lbergelson deleted the lb_fasta_gz branch September 20, 2018 02:27

		@@ -94,26 +77,21 @@ public CachingIndexedFastaSequenceFile(final Path fasta, final FastaSequenceInde
		* @param cacheSize the size of the cache to use in this CachingIndexedFastaReader, must be >= 0
		* @param preserveCase If true, we will keep the case of the underlying bases in the FASTA, otherwise everything is converted to upper case

		@@ -136,7 +114,7 @@ public CachingIndexedFastaSequenceFile(final Path fasta) throws FileNotFoundExce
		* @param fasta The file to open.
		* @param preserveCase If true, we will keep the case of the underlying bases in the FASTA, otherwise everything is converted to upper case

fasta gz support in gatk #5120

fasta gz support in gatk #5120

Conversation

lbergelson commented Aug 16, 2018

codecov-io commented Aug 22, 2018 • edited

Codecov Report

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

droazen Aug 28, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sooheelee commented Sep 7, 2018

droazen commented Sep 12, 2018

lbergelson commented Sep 12, 2018

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

lbergelson commented Sep 20, 2018

lbergelson commented Sep 20, 2018

codecov-io commented Aug 22, 2018 •

edited

droazen Aug 28, 2018 •

edited