Skip to content

Commit

Permalink
Fix CRAM reference region transitions.
Browse files Browse the repository at this point in the history
  • Loading branch information
cmnbroad committed Apr 10, 2024
1 parent 222d10b commit 2d5b20f
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 1 deletion.
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));
}
if ((referenceIndex != this.referenceIndex) ||
regionStart != 0 ||
(regionLength < referenceBases.length)) {
(regionLength != newSequenceRecord.getSequenceLength())) {
setCurrentSequence(referenceIndex);
referenceBases = referenceSource.getReferenceBases(sequenceRecord, true);
if (referenceBases == null) {
Expand Down
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() {
// 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);
}

@Test
public void testGetReferenceBasesByRegionExceedsContigLength() {
int regionStart = CRAMStructureTestHelper.REFERENCE_CONTIG_LENGTH / 2;
Expand Down

0 comments on commit 2d5b20f

Please sign in to comment.