feat: Add sanity check to indexed FASTA file

src/main/java/htsjdk/samtools/reference/AbstractIndexedFastaSequenceFile.java

@@ -214,7 +214,7 @@ public ReferenceSequence getSubsequenceAt(String contig, long start, long stop)
         final int bytesPerLine = indexEntry.getBytesPerLine();
         final int terminatorLength = bytesPerLine - basesPerLine;
 
-        long startOffset = ((start - 1) / basesPerLine) * bytesPerLine + (start - 1) % basesPerLine;
+        long startOffset = indexEntry.getOffset(start);
         // Cast to long so the second argument cannot overflow a signed integer.
         final long minBufferSize = Math.min(
                 (long) Defaults.NON_ZERO_BUFFER_SIZE, (long) (length / basesPerLine + 2) * (long) bytesPerLine);

src/main/java/htsjdk/samtools/reference/FastaSequenceIndex.java

@@ -33,8 +33,10 @@
 import java.io.PrintStream;
 import java.nio.file.Files;
 import java.nio.file.Path;
+import java.util.ArrayList;
+import java.util.HashMap;
 import java.util.Iterator;
-import java.util.LinkedHashMap;
+import java.util.List;
 import java.util.Map;
 import java.util.Objects;
 import java.util.Scanner;
@@ -44,11 +46,11 @@
  * Reads/writes a fasta index file (.fai), as generated by `samtools faidx`.
  */
 public class FastaSequenceIndex implements Iterable<FastaSequenceIndexEntry> {
-    /**
-     * Store the entries.  Use a LinkedHashMap for consistent iteration in insertion order.
-     */
-    private final Map<String, FastaSequenceIndexEntry> sequenceEntries =
-            new LinkedHashMap<String, FastaSequenceIndexEntry>();
+    /** Entries in the order they appear in the index file, for ordered iteration and O(1) positional access. */
+    private final List<FastaSequenceIndexEntry> entries = new ArrayList<>();
+
+    /** Entries keyed by contig name, for O(1) lookup by name. */
+    private final Map<String, FastaSequenceIndexEntry> entriesByContig = new HashMap<>();
 
     /**
      * Build a sequence index from the specified file.
@@ -82,10 +84,10 @@ protected FastaSequenceIndex() {}
      * @param indexEntry New index entry to add.
      */
     protected void add(FastaSequenceIndexEntry indexEntry) {
-        final FastaSequenceIndexEntry ret = sequenceEntries.put(indexEntry.getContig(), indexEntry);
-        if (ret != null) {
+        if (entriesByContig.putIfAbsent(indexEntry.getContig(), indexEntry) != null) {
             throw new SAMException("Contig '" + indexEntry.getContig() + "' already exists in fasta index.");
         }
+        entries.add(indexEntry);
     }
 
     /**
@@ -94,9 +96,15 @@ protected void add(FastaSequenceIndexEntry indexEntry) {
      * @param newName New name for the index entry.
      */
     protected void rename(FastaSequenceIndexEntry entry, String newName) {
-        sequenceEntries.remove(entry.getContig());
+        // Check for a collision before mutating anything, so a failed rename leaves the index unchanged.
+        // (Renaming an entry to its own current name is a no-op, not a collision.)
+        if (!newName.equals(entry.getContig()) && entriesByContig.containsKey(newName)) {
+            throw new SAMException("Contig '" + newName + "' already exists in fasta index.");
+        }
+        entriesByContig.remove(entry.getContig());
         entry.setContig(newName);
-        add(entry);
+        entriesByContig.put(newName, entry);
+        // `entry` keeps its position in `entries` (it is the same object), so index order is preserved.
     }
 
     /**
@@ -123,7 +131,7 @@ public boolean equals(Object other) {
 
     @Override
     public int hashCode() {
-        return Objects.hash(sequenceEntries);
+        return Objects.hash(entries);
     }
 
     /**
@@ -167,15 +175,13 @@ private void parseIndexFile(InputStream in) {
      */
     public void write(final Path indexFile) throws IOException {
         try (final PrintStream writer = new PrintStream(Files.newOutputStream(indexFile))) {
-            sequenceEntries
-                    .values()
-                    .forEach(se -> writer.println(String.join(
-                            "\t",
-                            se.getContig(),
-                            String.valueOf(se.getSize()),
-                            String.valueOf(se.getLocation()),
-                            String.valueOf(se.getBasesPerLine()),
-                            String.valueOf(se.getBytesPerLine()))));
+            entries.forEach(se -> writer.println(String.join(
+                    "\t",
+                    se.getContig(),
+                    String.valueOf(se.getSize()),
+                    String.valueOf(se.getLocation()),
+                    String.valueOf(se.getBasesPerLine()),
+                    String.valueOf(se.getBytesPerLine()))));
         }
     }
 
@@ -185,7 +191,7 @@ public void write(final Path indexFile) throws IOException {
      * @return True if contig name is present; false otherwise.
      */
     public boolean hasIndexEntry(String contigName) {
-        return sequenceEntries.containsKey(contigName);
+        return entriesByContig.containsKey(contigName);
     }
 
     /**
@@ -197,7 +203,7 @@ public boolean hasIndexEntry(String contigName) {
     public FastaSequenceIndexEntry getIndexEntry(String contigName) {
         if (!hasIndexEntry(contigName)) throw new SAMException("Unable to find entry for contig: " + contigName);
 
-        return sequenceEntries.get(contigName);
+        return entriesByContig.get(contigName);
     }
 
     /**
@@ -206,14 +212,27 @@ public FastaSequenceIndexEntry getIndexEntry(String contigName) {
      */
     @Override
     public Iterator<FastaSequenceIndexEntry> iterator() {
-        return sequenceEntries.values().iterator();
+        return entries.iterator();
     }
 
     /**
      * Returns the number of elements in the index.
      * @return Number of elements in the index.
      */
     public int size() {
-        return sequenceEntries.size();
+        return entries.size();
+    }
+
+    /**
+     * Returns the last entry in the index, i.e. the one for the contig that appears last in the index file.
+     *
+     * @return the last {@link FastaSequenceIndexEntry} in the index.
+     * @throws SAMException if the index is empty.
+     */
+    protected FastaSequenceIndexEntry getLastIndexEntry() {
+        if (entries.isEmpty()) {
+            throw new SAMException("Cannot determine the last entry of an empty fasta index.");
+        }
+        return entries.get(entries.size() - 1);
     }
 }

src/main/java/htsjdk/samtools/reference/FastaSequenceIndexEntry.java

@@ -108,6 +108,18 @@ public int getSequenceIndex() {
         return sequenceIndex;
     }
 
+    /** Return the offset to pos as determined by the number of bases and bytes per line
+     *
+     * @param pos the (1-based) position in the contig that is requested
+     * @return the offset (0-based) from 'location' where pos is located in the file.
+     */
+    public long getOffset(long pos) {
+        final int basesPerLine = this.getBasesPerLine();
+        final int bytesPerLine = this.getBytesPerLine();
+
+        return ((pos - 1) / basesPerLine) * bytesPerLine + ((pos - 1) % basesPerLine);
+    }
+
     /**
      * For debugging.  Emit the contents of each contig line.
      * @return A string representation of the contig line.

src/main/java/htsjdk/samtools/reference/IndexedFastaSequenceFile.java

@@ -61,6 +61,7 @@ public IndexedFastaSequenceFile(final Path path, final FastaSequenceIndex index)
                 throw new SAMException("Indexed block-compressed FASTA file cannot be handled: " + path);
             }
             this.channel = Files.newByteChannel(path);
+            sanityCheckFastaAgainstIndex(path, index);
         } catch (IOException e) {
             throw new SAMException("FASTA file should be readable but is not: " + path, e);
         }
@@ -123,6 +124,63 @@ public static boolean canCreateIndexedFastaReader(final Path fastaFile) {
         }
     }
 
+    /**
+     * Do some basic checking to make sure the fasta and its index are consistent with one another.
+     * <p>
+     * Verifies that the fasta file is at least as long as the last base position claimed by the index, and
+     * that any bytes beyond that last base are only whitespace. A mismatch usually indicates a corrupt fasta
+     * or a stale/incorrect index (which would otherwise silently return the wrong bases), and should be
+     * resolved by reindexing the fasta.
+     * <p>
+     * This is a cheap, necessary-but-not-sufficient check: it reads only the file length and the trailing
+     * bytes, so it catches truncation, wrong-file, and appended-content mismatches, but not an index whose
+     * internal offsets are wrong yet happen to end at the right place. It applies only to plain (non
+     * block-compressed) fasta files, where the on-disk length equals the uncompressed length; block-compressed
+     * references (see {@link BlockCompressedIndexedFastaSequenceFile}) are not checked.
+     *
+     * @param fastaFile Path to the fasta file
+     * @param fastaSequenceIndex the index to check against the fasta file
+     * @throws IOException if an io-error occurs while reading fastaFile
+     */
+    private void sanityCheckFastaAgainstIndex(final Path fastaFile, final FastaSequenceIndex fastaSequenceIndex)
+            throws IOException {
+
+        final FastaSequenceIndexEntry lastSequence = fastaSequenceIndex.getLastIndexEntry();
+        // 0-based byte offset of the last base of the last contig; that base occupies exactly this one byte.
+        final long lastBaseOffset = lastSequence.getLocation() + lastSequence.getOffset(lastSequence.getSize());
+        final long fastaLength = Files.size(fastaFile);
+
+        // The file must actually contain the last base byte, i.e. be strictly longer than its offset.
+        if (lastBaseOffset >= fastaLength) {
+            throw new IllegalArgumentException(
+                    ("The fasta file (%s) is shorter (%d bytes) than its index claims: the last base is expected at "
+                                    + "byte %d. Please reindex the fasta.")
+                            .formatted(fastaFile.toUri(), fastaLength, lastBaseOffset));
+        }
+
+        // Everything strictly after the last base (its line terminator plus any trailing blank lines) must be
+        // whitespace; a non-whitespace byte there means the index does not account for all of the fasta's content.
+        long position = lastBaseOffset + 1;
+        final ByteBuffer buffer = ByteBuffer.allocate(100);
+        while (position < fastaLength) {
+            buffer.clear();
+            final int bytesRead = readFromPosition(buffer, position);
+            if (bytesRead <= 0) {
+                break;
+            }
+            for (int i = 0; i < bytesRead; i++) {
+                final byte b = buffer.get(i);
+                if (!Character.isWhitespace((char) b)) {
+                    throw new IllegalArgumentException(
+                            ("The fasta file (%s) is longer than its index accounts for: found a non-whitespace "
+                                            + "character (%c) at byte %d, past the last base at byte %d. Please reindex the fasta.")
+                                    .formatted(fastaFile.toUri(), (char) b, position + i, lastBaseOffset));
+                }
+            }
+            position += bytesRead;
+        }
+    }
+
     /**
      * Reads a sequence of bytes from this channel into the given buffer,
      * starting at the given file position.

src/test/java/htsjdk/samtools/reference/AbstractIndexedFastaSequenceFileTest.java

@@ -58,6 +58,13 @@ public class AbstractIndexedFastaSequenceFileTest extends HtsjdkTest {
     private static final Path SEQUENCE_FILE_GZI = TEST_DATA_DIR.resolve("Homo_sapiens_assembly18.trimmed.fasta.gz.gzi");
     private static final Path SEQUENCE_FILE_NODICT =
             TEST_DATA_DIR.resolve("Homo_sapiens_assembly18.trimmed.nodict.fasta");
+    private static final Path HEADER_WITH_WHITESPACE = TEST_DATA_DIR.resolve("header_with_white_space.fasta");
+    private static final Path HEADER_WITH_EXTRA_WHITESPACE =
+            TEST_DATA_DIR.resolve("header_with_extra_white_space.fasta");
+    private static final Path CRLF = TEST_DATA_DIR.resolve("crlf.fasta");
+    private static final Path HEADER_WITH_WHITESPACE_INDEX =
+            AbstractIndexedFastaSequenceFile.findFastaIndex(HEADER_WITH_WHITESPACE);
+    private static final Path CRLF_INDEX = AbstractIndexedFastaSequenceFile.findFastaIndex(CRLF);
 
     private final String firstBasesOfChrM = "GATCACAGGTCTATCACCCT";
     private final String extendedBasesOfChrM =
@@ -66,6 +73,14 @@ public class AbstractIndexedFastaSequenceFileTest extends HtsjdkTest {
     private final String lastBasesOfChr20 = "ttgtctgatgctcatattgt";
     private final int CHR20_LENGTH = 1000000;
 
+    @DataProvider(name = "mismatched_indexes")
+    public Object[][] provideMismatchedIndexes() {
+        return new Object[][] {
+            {HEADER_WITH_WHITESPACE, CRLF_INDEX},
+            {CRLF, HEADER_WITH_WHITESPACE_INDEX}
+        };
+    }
+
     @DataProvider(name = "homosapiens")
     public Object[][] provideSequenceFile() throws FileNotFoundException {
         return new Object[][] {
@@ -518,4 +533,59 @@ public void testBlockCompressedIndexedFastaSequenceFileFromNio() throws IOExcept
             Assert.assertEquals(rs.getBaseString(), "CACAGGT");
         }
     }
+
+    @Test(expectedExceptions = IllegalArgumentException.class, dataProvider = "mismatched_indexes")
+    public void testWrongIndex(final Path fasta, final Path index) throws IOException {
+        // Opening a fasta with an index that does not match it should fail the sanity check.
+        try (IndexedFastaSequenceFile ignored = new IndexedFastaSequenceFile(fasta, new FastaSequenceIndex(index))) {
+            Assert.fail("Expected an IllegalArgumentException for a mismatched fasta/index.");
+        }
+    }
+
+    @Test
+    public void testExtraWhitespace() throws IOException {
+        // Trailing whitespace beyond the last base (here, extra blank lines) should be tolerated.
+        try (IndexedFastaSequenceFile ignored = new IndexedFastaSequenceFile(
+                HEADER_WITH_EXTRA_WHITESPACE, new FastaSequenceIndex(HEADER_WITH_WHITESPACE_INDEX))) {
+            // no exception expected
+        }
+    }
+
+    @Test
+    public void testSanityCheckAcceptsFastaWithoutTrailingNewline() throws IOException {
+        // A single-line fasta whose final line has no trailing newline (terminator length 0) is valid and
+        // must open without a false positive from the sanity check.
+        final Path dir = Files.createTempDirectory("fastaSanity");
+        final Path fasta = dir.resolve("noNewline.fasta");
+        try {
+            Files.write(fasta, ">c\nACGTACGT".getBytes(java.nio.charset.StandardCharsets.US_ASCII));
+            final FastaSequenceIndex index = FastaSequenceIndexCreator.buildFromFasta(fasta);
+            try (IndexedFastaSequenceFile ignored = new IndexedFastaSequenceFile(fasta, index)) {
+                // no exception expected
+            }
+        } finally {
+            Files.deleteIfExists(fasta);
+            Files.deleteIfExists(dir);
+        }
+    }
+
+    @Test(expectedExceptions = IllegalArgumentException.class)
+    public void testSanityCheckRejectsTruncatedFasta() throws IOException {
+        // Truncating the fasta so it ends exactly where the last base should start (the boundary case) must
+        // be detected against an index that still claims the base is present.
+        final Path dir = Files.createTempDirectory("fastaSanity");
+        final Path fasta = dir.resolve("trunc.fasta");
+        try {
+            Files.write(fasta, ">c\nACGTACGT\n".getBytes(java.nio.charset.StandardCharsets.US_ASCII));
+            final FastaSequenceIndex index = FastaSequenceIndexCreator.buildFromFasta(fasta); // index for the full file
+            Files.write(
+                    fasta, ">c\nACGTACG".getBytes(java.nio.charset.StandardCharsets.US_ASCII)); // drop the last base
+            try (IndexedFastaSequenceFile ignored = new IndexedFastaSequenceFile(fasta, index)) {
+                Assert.fail("Expected the sanity check to reject a truncated fasta.");
+            }
+        } finally {
+            Files.deleteIfExists(fasta);
+            Files.deleteIfExists(dir);
+        }
+    }
 }

src/test/java/htsjdk/samtools/reference/FastaSequenceIndexTest.java

@@ -309,4 +309,46 @@ public void testWrite() throws Exception {
         final FastaSequenceIndex writtenIndex = new FastaSequenceIndex(fileToWrite);
         Assert.assertEquals(writtenIndex, originalIndex);
     }
+
+    private static FastaSequenceIndex twoContigIndex() {
+        final FastaSequenceIndex index = new FastaSequenceIndex();
+        index.add(new FastaSequenceIndexEntry("chr1", 0, 100, 60, 61, 0));
+        index.add(new FastaSequenceIndexEntry("chr2", 200, 100, 60, 61, 1));
+        return index;
+    }
+
+    @Test
+    public void testRenameToNewName() {
+        final FastaSequenceIndex index = twoContigIndex();
+        index.rename(index.getIndexEntry("chr1"), "chrOne");
+
+        Assert.assertFalse(index.hasIndexEntry("chr1"));
+        Assert.assertTrue(index.hasIndexEntry("chrOne"));
+        Assert.assertEquals(index.getIndexEntry("chrOne").getContig(), "chrOne");
+        Assert.assertEquals(index.size(), 2);
+        // renaming preserves list position: the renamed entry is still first.
+        Assert.assertEquals(index.iterator().next().getContig(), "chrOne");
+    }
+
+    @Test
+    public void testRenameToOwnNameIsNoOp() {
+        final FastaSequenceIndex index = twoContigIndex();
+        index.rename(index.getIndexEntry("chr1"), "chr1");
+
+        Assert.assertTrue(index.hasIndexEntry("chr1"));
+        Assert.assertEquals(index.size(), 2);
+    }
+
+    @Test
+    public void testRenameToExistingNameThrowsAndLeavesIndexUnchanged() {
+        final FastaSequenceIndex index = twoContigIndex();
+
+        Assert.assertThrows(SAMException.class, () -> index.rename(index.getIndexEntry("chr1"), "chr2"));
+
+        // A failed rename must not corrupt the index (the name lookup and the ordered list must stay in sync).
+        Assert.assertTrue(index.hasIndexEntry("chr1"), "chr1 must remain looked-up-able after a failed rename");
+        Assert.assertTrue(index.hasIndexEntry("chr2"));
+        Assert.assertEquals(index.getIndexEntry("chr1").getContig(), "chr1");
+        Assert.assertEquals(index.size(), 2);
+    }
 }

src/test/resources/htsjdk/samtools/reference/header_with_extra_white_space.fasta

@@ -0,0 +1,7 @@
+>a test white space
+ACTG
+>b test whitespace
+ACTG
+
+
+