Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix CRAMReferenceRegion updating. #1708

Merged
merged 9 commits into from
Jun 4, 2024
Merged

Fix CRAMReferenceRegion updating. #1708

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
2d5b20f
Fix CRAM reference region transitions.
cmnbroad Apr 9, 2024
969bdf1
Add roundtrip tests that fail without this change.
cmnbroad Apr 10, 2024
4f465ab
Add roundtrip tests with 2 and 3 containers aligned to position 1.
cmnbroad May 13, 2024
f276317
Add more roundtrip tests that validate bases.
cmnbroad May 13, 2024
64180f8
Add a large CRAM roundtrip test and samtools rountdtip test.
cmnbroad May 14, 2024
19be0c7
Always remember your index files.
cmnbroad May 14, 2024
669fd7e
Code review comments.
cmnbroad Jun 3, 2024
ca12df0
One more change based on review.
cmnbroad Jun 4, 2024
7c473de
One more update.
cmnbroad Jun 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion src/main/java/htsjdk/samtools/cram/build/CRAMReferenceRegion.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,9 +113,14 @@ public void fetchReferenceBases(final int referenceIndex) {

// Re-resolve the reference bases if we don't have a current region or if the region we have
// doesn't span the *entire* contig requested.
final SAMSequenceRecord newSequenceRecord = sequenceDictionary.getSequence(referenceIndex);
if (newSequenceRecord == null) {
throw new IllegalArgumentException(
String.format("Requested reference sequence index %d not found", referenceIndex));
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you call the existing getSAMSequenceRecord() method here? It does exactly the same thing.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh yeah, done.

if ((referenceIndex != this.referenceIndex) ||
regionStart != 0 ||
(regionLength < referenceBases.length)) {
(regionLength != newSequenceRecord.getSequenceLength())) {
setCurrentSequence(referenceIndex);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment that the setCurrentSequence() call needs to happen first in this if block (in particular, before the assignment to regionLength below). Or alternatively, you could have regionLength = newSequenceRecord.getSequenceLength(); below to eliminate the dependency on order of operations.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Latter is done.

referenceBases = referenceSource.getReferenceBases(sequenceRecord, true);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is still a potential order-of-operations issue here: since the call to getReferenceBases() uses the instance variable sequenceRecord, the call to setCurrentSequence() must happen first. If the order of the statements were reversed in some future refactoring, we might have another terrible bug. Either add a comment explaining that setCurrentSequence() must happen first, or eliminate the order dependency by using newSequenceRecord in the getReferenceBases() call.

if (referenceBases == null) {
Expand Down
79 changes: 65 additions & 14 deletions src/test/java/htsjdk/samtools/CRAMAllEncodingStrategiesTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,32 +21,72 @@ public class CRAMAllEncodingStrategiesTest extends HtsjdkTest {
private static final File TEST_DATA_DIR = new File("src/test/resources/htsjdk/samtools/cram");
private final CompressorCache compressorCache = new CompressorCache();

@DataProvider(name="roundTripTestFiles")
public Object[][] roundTripTestFiles() {
@DataProvider(name="defaultStrategyRoundTripTestFiles")
public Object[][] defaultStrategyRoundTripTestFiles() {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain why some of these test cases have to use lenient equality, while others don't?

Copy link
Collaborator Author

@cmnbroad cmnbroad Jun 3, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

CRAMs don't always roundtrip with complete fidelity, especially when roundtripping through samtools, which resurrects NM/MD tags whether they were present in the original or not, unless they are specifically excluded. So in some cases when roundtripping, you can't use full SAMRecord comparisons. In those cases we use lenient equality to test only a subset of the SAMRecord. I've added this as a comment.

return new Object[][] {
{ new File(TEST_DATA_DIR, "NA12878.20.21.1-100.100-SeqsPerSlice.500-unMapped.cram"),
new File(TEST_DATA_DIR, "human_g1k_v37.20.21.1-100.fasta") },
new File(TEST_DATA_DIR, "human_g1k_v37.20.21.1-100.fasta"),
false, false },
{ new File(TEST_DATA_DIR, "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.10m-10m100.cram"),
new File("src/test/resources/htsjdk/samtools/reference/human_g1k_v37.20.21.fasta.gz"),
false, false },
{ new File(TEST_DATA_DIR, "NA12878.unmapped.cram"),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add one-line comments explaining the provenance of these three crams, and what they are testing?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done with as much specificity as I have about these, since they've been in the repo for a while.

new File(TEST_DATA_DIR, "human_g1k_v37.20.21.1-100.fasta"),
false, false },
// generated with samtools 1.19 from the gatk bam file CEUTrio.HiSeq.WGS.b37.NA12878.20.21.bam
{ new File(TEST_DATA_DIR, "CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.0.samtools.cram"),
new File("src/test/resources/htsjdk/samtools/reference/human_g1k_v37.20.21.fasta.gz"),
true, false },
// these tests use lenient equality to only validate read names, bases and qual scores
{ new File(TEST_DATA_DIR, "mitoAlignmentStartTestGATKGen.cram"),
new File(TEST_DATA_DIR, "mitoAlignmentStartTest.fa"), true, false },
{ new File(TEST_DATA_DIR, "mitoAlignmentStartTest.cram"),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add comments explaining what these files are as well

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done with as much specificity as I have about these, since they've been in the repo for a while.

new File(TEST_DATA_DIR, "mitoAlignmentStartTest.fa"), true, false },
// files created by rewriting the htsjdk test file src/test/resources/htsjdk/samtools/cram/mitoAlignmentStartTest.cram
// using code that replicates the first read (which is aligned to position 1 of the mito contig) either
// 10,000 or 20,000 times, to create a file with 2 or 3 containers, respectively, that have reads aligned to
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you confirm via direct inspection of the file that it did in fact create 2 or 3 containers with reads mapped to position 1?

Copy link
Collaborator Author

@cmnbroad cmnbroad Jun 3, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. You can also see this in the GATK PR, which uses files that were created using the same code that created these files, to test the detector on files with multiple bad containers. Alternatively, you can run GATK PrintFileDiagnostics on these files and see the container alignments.

// position 1 of the contig
{ new File(TEST_DATA_DIR, "mitoAlignmentStartTest_2_containers_aligned_to_pos_1.cram"),
new File(TEST_DATA_DIR, "mitoAlignmentStartTest.fa"), true, false },
{ new File(TEST_DATA_DIR, "mitoAlignmentStartTest_3_containers_aligned_to_pos_1.cram"),
new File(TEST_DATA_DIR, "mitoAlignmentStartTest.fa"), true, false }
};
}

@Test(dataProvider = "roundTripTestFiles")
public final void testRoundTripDefaultEncodingStrategy(final File sourceFile, final File referenceFile) throws IOException {
@Test(dataProvider = "defaultStrategyRoundTripTestFiles")
public final void testRoundTripDefaultEncodingStrategy(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment before this test method describing the purpose of the test (since we're not just testing "encoding strategies" here....)

Copy link
Collaborator Author

@cmnbroad cmnbroad Jun 3, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is testing the default encoding strategy. Comment added.

final File sourceFile,
final File referenceFile,
final boolean lenientEquality,
final boolean emitDetail) throws IOException {
final CRAMEncodingStrategy testStrategy = new CRAMEncodingStrategy();
final File tempOutCRAM = File.createTempFile("testRoundTrip", ".cram");
tempOutCRAM.deleteOnExit();
CRAMTestUtils.writeToCRAMWithEncodingStrategy(testStrategy, sourceFile, tempOutCRAM, referenceFile);
assertRoundTripFidelity(sourceFile, tempOutCRAM, referenceFile, false);
assertRoundtripFidelityWithSamtools(tempOutCRAM, referenceFile);
assertRoundTripFidelity(sourceFile, tempOutCRAM, referenceFile, lenientEquality, emitDetail);
assertRoundtripFidelityWithSamtools(tempOutCRAM, referenceFile, lenientEquality, emitDetail);
}

@DataProvider(name="encodingStrategiesTestFiles")
public Object[][] encodingStrategiesTestFiles() {
return new Object[][] {
{ new File(TEST_DATA_DIR, "NA12878.20.21.1-100.100-SeqsPerSlice.500-unMapped.cram"),
new File(TEST_DATA_DIR, "human_g1k_v37.20.21.1-100.fasta"), false, false },
};
}

@Test(dataProvider = "roundTripTestFiles")
public final void testAllEncodingStrategyCombinations(final File cramSourceFile, final File referenceFile) throws IOException {
@Test(dataProvider = "encodingStrategiesTestFiles")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would there be value in having this testAllEncodingStrategyCombinations test use the defaultStrategyRoundTripTestFiles DataProvider as well, with its many additional test cases?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps, but there are currently 81 encoding strategies, so using all of those files would make this test take a really long time, and I would guess, dominate the CI test time, especially now that we have the large CRAM test file in that list. So instead I separated them this way.

public final void testAllEncodingStrategyCombinations(
final File cramSourceFile,
final File referenceFile,
final boolean lenientEquality,
final boolean emitDetail) throws IOException {
for (final Tuple<String, CRAMEncodingStrategy> testStrategy : getAllEncodingStrategies()) {
final File tempOutCRAM = File.createTempFile("allEncodingStrategyCombinations", ".cram");
tempOutCRAM.deleteOnExit();
CRAMTestUtils.writeToCRAMWithEncodingStrategy(testStrategy.b, cramSourceFile, tempOutCRAM, referenceFile);
assertRoundTripFidelity(cramSourceFile, tempOutCRAM, referenceFile, false);
assertRoundtripFidelityWithSamtools(tempOutCRAM, referenceFile);
assertRoundTripFidelity(cramSourceFile, tempOutCRAM, referenceFile, lenientEquality, emitDetail);
assertRoundtripFidelityWithSamtools(tempOutCRAM, referenceFile, lenientEquality, emitDetail);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment explaining why we also do a test with samtools

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Its basically testing interoperability with samtools for all of the different encodings, to make sure both tool sets agree on the data. Added a comment saying its for interop testing.

}
}

Expand Down Expand Up @@ -82,6 +122,7 @@ public void assertRoundTripFidelity(
final File sourceFile,
final File targetCRAMFile,
final File referenceFile,
final boolean lenientEquality,
final boolean emitDetail) throws IOException {
try (final SamReader sourceReader = SamReaderFactory.makeDefault()
.referenceSequence(referenceFile)
Expand All @@ -91,7 +132,13 @@ public void assertRoundTripFidelity(
final SAMRecordIterator sourceIterator = sourceReader.iterator();
final SAMRecordIterator targetIterator = copyReader.getIterator();
while (sourceIterator.hasNext() && targetIterator.hasNext()) {
if (emitDetail) {
if (lenientEquality) {
final SAMRecord sourceRec = sourceIterator.next();
final SAMRecord targetRec = targetIterator.next();
Assert.assertEquals(targetRec.getReadName(), sourceRec.getReadName());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can lenient mode also check the start/end positions of the read?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, added.

Assert.assertEquals(targetRec.getReadBases(), sourceRec.getReadBases());
Assert.assertEquals(targetRec.getBaseQualities(), sourceRec.getBaseQualities());
} else if (emitDetail) {
final SAMRecord sourceRec = sourceIterator.next();
final SAMRecord targetRec = targetIterator.next();
if (!sourceRec.equals(targetRec)) {
Expand All @@ -108,13 +155,17 @@ public void assertRoundTripFidelity(
}
}

private void assertRoundtripFidelityWithSamtools(final File sourceCRAM, final File referenceFile) throws IOException {
private void assertRoundtripFidelityWithSamtools(
final File sourceCRAM,
final File referenceFile,
final boolean lenientEquality,
final boolean emitDetail) throws IOException {
if (SamtoolsTestUtils.isSamtoolsAvailable()) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If samtools is not available, the test probably shouldn't just silently do nothing and succeed. Maybe throw a TestNG SkipException in that case?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense. Done.

final File samtoolsOutFile = SamtoolsTestUtils.convertToCRAM(
sourceCRAM,
referenceFile,
"--input-fmt-option decode_md=0 --output-fmt-option store_md=0 --output-fmt-option store_nm=0");
assertRoundTripFidelity(sourceCRAM, samtoolsOutFile, referenceFile, false);
assertRoundTripFidelity(sourceCRAM, samtoolsOutFile, referenceFile, lenientEquality, emitDetail);
}
}

Expand Down
34 changes: 34 additions & 0 deletions src/test/java/htsjdk/samtools/cram/ref/CRAMReferenceRegionTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import htsjdk.HtsjdkTest;
import htsjdk.samtools.SAMSequenceRecord;
import htsjdk.samtools.cram.build.CRAMReferenceRegion;
import htsjdk.samtools.cram.structure.AlignmentContext;
import htsjdk.samtools.cram.structure.CRAMStructureTestHelper;
import htsjdk.samtools.reference.InMemoryReferenceSequenceFile;
import org.testng.Assert;
Expand Down Expand Up @@ -167,6 +168,39 @@ public void testGetReferenceBasesByRegionPositive(
Assert.assertEquals(bases, Arrays.copyOfRange(fullContigBases, requestedOffset, requestedOffset + requestedLength));
}

// simulate the state transitions that occur when writing a CRAM file
@Test
public void testSerialStateTransitions() {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment here explaining and referencing the bug that prompted us to write this test.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

// start with an entire reference sequence
final CRAMReferenceRegion cramReferenceRegion = getAlternatingReferenceRegion();
cramReferenceRegion.fetchReferenceBases(CRAMStructureTestHelper.REFERENCE_SEQUENCE_ZERO);
final long fullRegionFragmentLength = cramReferenceRegion.getRegionLength();
Assert.assertEquals(fullRegionFragmentLength, CRAMStructureTestHelper.REFERENCE_CONTIG_LENGTH);

// transition to a shorter reference fragment using fetchReferenceBasesByRegion, then back to the full region
final int SHORT_FRAGMENT_LENGTH = 5;
Assert.assertTrue(SHORT_FRAGMENT_LENGTH < fullRegionFragmentLength);
cramReferenceRegion.fetchReferenceBasesByRegion(CRAMStructureTestHelper.REFERENCE_SEQUENCE_ZERO, 0, SHORT_FRAGMENT_LENGTH);
Assert.assertEquals(cramReferenceRegion.getRegionLength(), SHORT_FRAGMENT_LENGTH);

// now transition back to the full sequence
cramReferenceRegion.fetchReferenceBases(CRAMStructureTestHelper.REFERENCE_SEQUENCE_ZERO);
Assert.assertEquals(cramReferenceRegion.getRegionLength(), fullRegionFragmentLength);

// transition to a shorter region fragment length using fetchReferenceBasesByRegion(AlignmentContext), then back to the full region
Assert.assertTrue(SHORT_FRAGMENT_LENGTH < fullRegionFragmentLength);
cramReferenceRegion.fetchReferenceBasesByRegion(
new AlignmentContext(
new ReferenceContext(CRAMStructureTestHelper.REFERENCE_SEQUENCE_ZERO),
1,
SHORT_FRAGMENT_LENGTH));
Assert.assertEquals(cramReferenceRegion.getRegionLength(), SHORT_FRAGMENT_LENGTH);

// now transition back to the full sequence
cramReferenceRegion.fetchReferenceBases(CRAMStructureTestHelper.REFERENCE_SEQUENCE_ZERO);
Assert.assertEquals(cramReferenceRegion.getRegionLength(), fullRegionFragmentLength);
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is CRAMReferenceRegion ever used across multiple threads? Give how the fetch methods update the internal state of the object on every call, it's clearly not safe to use in a parallelized manner. Do we need to make these methods synchronized?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is currently not used anywhere across threads. As you noted, it's not thread-safe. I guess I could add synchronized, but that still wouldn't make it thread-safe, since the class is is stateful, and the usage pattern is fetch... followed by one or more getCurrent.. methods. So synchronized won't fix the statefulness.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added a comment saying it's not thread-safe.


@Test
public void testGetReferenceBasesByRegionExceedsContigLength() {
int regionStart = CRAMStructureTestHelper.REFERENCE_CONTIG_LENGTH / 2;
Expand Down
Binary file added ...est/resources/htsjdk/samtools/cram/CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.0.samtools.cram
View file
Open in desktop
Binary file not shown.
Binary file added ...esources/htsjdk/samtools/cram/CEUTrio.HiSeq.WGS.b37.NA12878.20.21.v3.0.samtools.cram.crai
View file
Open in desktop
Binary file not shown.
Binary file added .../resources/htsjdk/samtools/cram/mitoAlignmentStartTest_2_containers_aligned_to_pos_1.cram
View file
Open in desktop
Binary file not shown.
Binary file added .../resources/htsjdk/samtools/cram/mitoAlignmentStartTest_3_containers_aligned_to_pos_1.cram
View file
Open in desktop
Binary file not shown.
Loading