GATK VariantFiltration can't use CSI-indexed raw.vcf.gz -- breaks hard filtering for long-contig genomes

[tsackton]

Summary

rule variant_filtration (added in #334) cannot succeed for any genome with a contig longer than the TBI limit (2^29 - 1 bp, ~512Mb) when variant_calling.long_contig_mode is (or auto-resolves to) true. GATK's VariantFiltration never auto-discovers a .csi index for its -V input -- it only checks for .tbi -- so there is no index format that both (a) can represent such a contig's coordinates and (b) GATK will actually recognize.

Reproduction

Encountered while backfilling results/vcfs/filtered.vcf.gz for existing runs per #332. One run's reference (GCA_965641915.1) has a contig >512Mb, so variant_calling.long_contig_mode correctly resolves to true and raw.vcf.gz is indexed with bcftools index -c (.csi), matching what rule concat_interval_vcfs / rule compress_interval_raw_vcf in gatk_intervals.smk already do for long-contig mode.

Running the variant_filtration shell command as Snakemake itself would build it:

gatk VariantFiltration -R reference.fa.gz -V raw.vcf.gz --output filtered.vcf.gz ...

fails immediately with:

A USER ERROR has occurred: An index is required but was not found for file drivingVariantFile:.../raw.vcf.gz.
Support for unindexed block-compressed files has been temporarily disabled. Try running IndexFeatureFile on the input.

This reproduces even on a small, non-long-contig VCF that's been indexed with bcftools index -c instead of -t -- i.e. it's not a contig-size problem for GATK's reader, it's that GATK's automatic index discovery for -V inputs simply never looks for .csi. Confirmed in htsjdk source: AbstractFeatureReader.isTabix() only checks for a .tbi companion file when no index path is explicitly supplied:

public static boolean isTabix(String resourcePath, String indexPath) throws IOException {
    if (indexPath == null) {
        indexPath = ParsingUtils.appendToPath(resourcePath, FileExtensions.TABIX_INDEX);
    }
    return IOUtil.hasBlockCompressedExtension(resourcePath) && ParsingUtils.resourceExists(indexPath);
}

GATK's own suggested workaround, IndexFeatureFile, doesn't help either -- for this same contig it crashes trying to build a .tbi (which can't represent the contig at all):

java.lang.ArrayIndexOutOfBoundsException: Index 37451 out of bounds for length 37451
	at htsjdk.samtools.BinningIndexBuilder.processFeature(BinningIndexBuilder.java:102)
	at htsjdk.tribble.index.tabix.TabixIndexCreator.finalizeFeature(TabixIndexCreator.java:106)
	at htsjdk.tribble.index.tabix.TabixIndexCreator.addFeature(TabixIndexCreator.java:92)
	at htsjdk.tribble.index.IndexFactory.createIndex(IndexFactory.java:529)
	at htsjdk.tribble.index.IndexFactory.createTabixIndex(IndexFactory.java:476)
	at org.broadinstitute.hellbender.tools.IndexFeatureFile.createAppropriateIndexInMemory(IndexFeatureFile.java:109)

(GATK version in use: gatk4 4.6.2.0, pinned via workflow/envs/gatk.yaml as gatk4>=4.1.)

Impact

Any genome routed into long_contig_mode (auto or explicit) that also uses a GATK-lineage caller with generate_filtered_vcf: true (the default) will hit this at the variant_filtration step. Since #334 made FILTERED_VCF a default rule all target for GATK-lineage callers, this now breaks full runs on such genomes out of the box, not just optional/manual filtering.

Possible directions

• Split variant_filtration (and any other GATK step reading RAW_VCF post long-contig-mode compression) by contig/region so no single query ever needs an index at all, if that's feasible for how VariantFiltration traverses.
• Investigate whether a newer GATK/htsjdk release adds .csi auto-discovery for -V inputs (did not find one as of 4.6.2.0).
• Explicitly pass an index via FeatureInput URI syntax if GATK's argument parser supports specifying an index path alongside -V (didn't find a public CLI option for this on VariantFiltration --help, but worth double-checking against engine internals).
• Fall back to a bcftools-based equivalent of the same hard-filter expressions for long-contig-mode genomes specifically, since bcftools can query/stream .csi-indexed files without this limitation.

Happy to share the exact commands/genome used to reproduce if helpful.