Skip to content

fix: tolerate duplicate VCF header definitions #17

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

Merged
tfenne merged 1 commit into from
Jun 15, 2026
Merged

fix: tolerate duplicate VCF header definitions #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
9 changes: 9 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

## [Unreleased]

### Fixed

- `simulate` and `methylate` now tolerate VCFs that redeclare a header ID
(e.g. `duplicate INFO ID: BREAKSIMLENGTH`). Such duplicates are common in
files from upstream tools and are accepted by bcftools; holodeck drops the
repeated definitions (keeping the first) instead of erroring. Compression is
now detected from the file's magic bytes rather than its extension, so the
variant reader also accepts plain-gzip and extension-less inputs.

## [0.3.0] - 2026-06-13

### Added
Expand Down
23 changes: 2 additions & 21 deletions src/vcf/methylation.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1237,25 +1237,6 @@ fn per_hap_bot_bit_string(
.join("|")
}

/// Open a VCF file as a buffered line reader, transparently decompressing
/// BGZF-/gzip-compressed inputs detected by the gzip magic bytes
/// `0x1f 0x8b` at the start of the file. Returns a boxed `BufRead` so
/// header-only probes and full-body readers can share the open/peek logic.
fn open_vcf_buf_reader(path: &std::path::Path) -> std::io::Result<Box<dyn std::io::BufRead>> {
use std::io::{BufReader, Read as _};
let mut peek_buf = [0u8; 2];
{
let mut f = std::fs::File::open(path)?;
f.read_exact(&mut peek_buf)?;
}
let file = std::fs::File::open(path)?;
if peek_buf == [0x1f, 0x8b] {
Ok(Box::new(BufReader::new(flate2::read::MultiGzDecoder::new(file))))
} else {
Ok(Box::new(BufReader::new(file)))
}
}

/// Check whether a VCF actually carries methylation truth: the header must
/// declare both `MT` and `MB` FORMAT fields **and** at least one data record
/// must list `MT` (or `MB`) in its FORMAT column.
Expand Down Expand Up @@ -1283,7 +1264,7 @@ fn open_vcf_buf_reader(path: &std::path::Path) -> std::io::Result<Box<dyn std::i
/// Returns an I/O error if the file cannot be opened or read.
pub fn vcf_has_mt_mb_records(path: &std::path::Path) -> std::io::Result<bool> {
use std::io::BufRead as _;
let mut reader = open_vcf_buf_reader(path)?;
let mut reader = crate::vcf::open_vcf_buf_reader(path)?;
let mut saw_top_strand_field = false;
let mut saw_bot_strand_field = false;
let mut line = String::new();
Expand Down Expand Up @@ -1361,7 +1342,7 @@ pub fn parse_methylation_vcf(path: &std::path::Path) -> std::io::Result<Methylat
// memory: the whole-file `String` and the per-contig `Vec<u8>`s held
// the body bytes twice. Streaming halves that to roughly one body
// worth (the `HashMap` values) plus a single `read_line` buffer.
let mut reader = open_vcf_buf_reader(path)?;
let mut reader = crate::vcf::open_vcf_buf_reader(path)?;

let mut saw_top_strand_field = false;
let mut saw_bot_strand_field = false;
Expand Down
202 changes: 191 additions & 11 deletions src/vcf/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@ pub mod genotype;
/// the `methylate` subcommand and the methylation-aware `simulate` path.
pub mod methylation;

use std::collections::HashMap;
use std::collections::{HashMap, HashSet};
use std::io::{BufRead, Read};
use std::path::Path;

use anyhow::{Context, Result, bail};
Expand Down Expand Up @@ -59,7 +60,8 @@ pub struct ParsedVariants {
/// ploidy even on contigs (or whole samples) with no ALT calls.
///
/// # Arguments
/// * `path` — Path to the VCF file (plain or BGZF; noodles autodetects).
/// * `path` — Path to the VCF file (plain text, gzip, or BGZF; detected from
/// the leading magic bytes).
/// * `sample_name` — Sample to resolve genotypes for, or `None` to use the
/// only sample in a single-sample VCF.
/// * `_dict` — Sequence dictionary (reserved for future cross-validation
Expand All @@ -73,11 +75,7 @@ pub fn parse_variants_by_contig(
sample_name: Option<&str>,
_dict: &SequenceDictionary,
) -> Result<ParsedVariants> {
let mut reader = vcf::io::reader::Builder::default()
.build_from_path(path)
.with_context(|| format!("Failed to open VCF: {}", path.display()))?;

let header = reader.read_header()?;
let (header, mut reader) = open_lenient_vcf_reader(path)?;
let sample_index = resolve_sample_index(&header, sample_name)?;

let mut by_contig: HashMap<String, Vec<VariantRecord>> = HashMap::new();
Expand Down Expand Up @@ -153,10 +151,7 @@ pub fn parse_variants_by_contig(
/// Returns an error if the VCF cannot be read or the sample configuration
/// is invalid.
pub(crate) fn validate_vcf_sample(path: &Path, sample_name: Option<&str>) -> Result<String> {
let mut reader = vcf::io::reader::Builder::default()
.build_from_path(path)
.with_context(|| format!("Failed to open VCF: {}", path.display()))?;
let header = reader.read_header()?;
let (header, _reader) = open_lenient_vcf_reader(path)?;
let idx = resolve_sample_index(&header, sample_name)?;
Ok(header
.sample_names()
Expand Down Expand Up @@ -206,6 +201,112 @@ fn resolve_sample_index(header: &vcf::Header, sample_name: Option<&str>) -> Resu
}
}

/// Open a VCF and return its parsed header plus a record reader, tolerating
/// duplicate meta-information definitions that noodles would otherwise reject.
///
/// noodles' header parser errors on any redeclared `##INFO`/`##FORMAT`/
/// `##FILTER`/`##ALT`/`##contig` ID (e.g. `duplicate INFO ID: BREAKSIMLENGTH`).
/// Real VCFs from upstream tools (Manta, DELLY, hand-merged files) routinely
/// carry such duplicates, and bcftools reads them without complaint; holodeck
/// must not be stricter. This reads the raw header text first, drops every
/// repeated structured definition (keeping the first, matching bcftools), and
/// hands the cleaned header back to noodles for parsing. The record stream is
/// untouched — only the header is rewritten.
///
/// Compression is detected from the leading magic bytes (plain text, gzip, or
/// BGZF), so the reader is also indifferent to the file's extension.
///
/// # Errors
/// Returns an error if the file cannot be opened or the (cleaned) header fails
/// to parse.
fn open_lenient_vcf_reader(
path: &Path,
) -> Result<(vcf::Header, vcf::io::Reader<Box<dyn BufRead>>)> {
let mut inner = open_vcf_buf_reader(path)
.with_context(|| format!("Failed to open VCF: {}", path.display()))?;

// Read raw header lines up to and including the `#CHROM` column header,
// dropping duplicate structured definitions. Everything after `#CHROM`
// stays unread in `inner` and is streamed as records below.
let mut header_text = String::new();
let mut seen: HashSet<(String, String)> = HashSet::new();
let mut line = String::new();
loop {
line.clear();
let n = inner
.read_line(&mut line)
.with_context(|| format!("Failed to read VCF header: {}", path.display()))?;
if n == 0 {
break; // EOF before `#CHROM`; let noodles report the malformed header
}
let is_meta = line.starts_with("##");
if is_meta
&& let Some(key) = structured_header_key(&line)
&& !seen.insert(key.clone())
{
log::warn!("dropping duplicate VCF header definition: ##{}=<ID={},...>", key.0, key.1);
Comment on lines +243 to +247

@coderabbitai coderabbitai Bot Jun 15, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟡 Minor | ⚡ Quick win

Constrain deduplication to the intended structured directives.

Current logic deduplicates any ##...=<ID=...> line, which can silently drop unrelated structured metadata. Limit dedup to the targeted directive set so non-target header semantics are preserved.

💡 Proposed fix 
-        if is_meta
-            && let Some(key) = structured_header_key(&line)
-            && !seen.insert(key.clone())
-        {
+        if is_meta
+            && let Some(key) = structured_header_key(&line)
+            && (matches!(key.0.as_str(), "INFO" | "FORMAT" | "FILTER" | "ALT")
+                || key.0.eq_ignore_ascii_case("contig"))
+            && !seen.insert(key.clone())
+        {
             log::warn!("dropping duplicate VCF header definition: ##{}=<ID={},...>", key.0, key.1);
             continue;
         }

Also applies to: 294-307

🤖 Prompt for AI Agents 
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.

In `@src/vcf/mod.rs` around lines 243 - 247, The deduplication logic in the
is_meta condition block is too broad and will drop any structured metadata line
matching the ##...=<ID=...> pattern, not just the intended structured
directives. Modify the logic to first check if the line contains one of the
targeted directive types before attempting deduplication, so that unrelated
structured header definitions are preserved. Apply this same constraint to all
affected deduplication checks (including the additional occurrences noted in the
review).

continue;
}
let is_column_header = !is_meta; // first non-`##` line is `#CHROM`
header_text.push_str(&line);
if is_column_header {
break;
}
}

// Replay the cleaned header in front of the still-unread record stream so
// noodles parses the header and reads records from one continuous source.
let chained: Box<dyn BufRead> =
Box::new(std::io::Cursor::new(header_text.into_bytes()).chain(inner));
let mut reader = vcf::io::Reader::new(chained);
let header = reader
.read_header()
.with_context(|| format!("Failed to parse VCF header: {}", path.display()))?;
Ok((header, reader))
}

/// Open a VCF for buffered reading, transparently decompressing gzip/BGZF.
///
/// Compression is detected from the leading two magic bytes (`1f 8b`) rather
/// than the file extension, and `MultiGzDecoder` handles both single-member
/// gzip and the multi-member BGZF blocks `bgzip`/`bcftools` produce.
///
/// # Errors
/// Returns an I/O error if the file cannot be opened or its first bytes read.
pub(crate) fn open_vcf_buf_reader(path: &Path) -> std::io::Result<Box<dyn BufRead>> {
use std::io::BufReader;

let mut peek = [0u8; 2];
// A file shorter than two bytes cannot be gzip; only treat the input as
// compressed when both magic bytes were actually read.
let bytes_read = {
let mut f = std::fs::File::open(path)?;
f.read(&mut peek)?
};
let file = std::fs::File::open(path)?;
if bytes_read == 2 && peek == [0x1f, 0x8b] {
Ok(Box::new(BufReader::new(flate2::read::MultiGzDecoder::new(file))))
} else {
Ok(Box::new(BufReader::new(file)))
}
}

/// Extract the `(directive, ID)` identity of a structured meta-information
/// line, e.g. `##INFO=<ID=DP,...>` → `("INFO", "DP")`.
///
/// Returns `None` for lines without a leading `ID=` field — `##fileformat`,
/// free-form `##key=value`, or the rare case where `ID` is not the first
/// attribute. Such lines are never deduplicated, matching the categories
/// noodles enforces unique (`INFO`/`FORMAT`/`FILTER`/`ALT`/`contig`), whose
/// emitters all write `ID` first.
fn structured_header_key(line: &str) -> Option<(String, String)> {
let rest = line.trim_start_matches('#').trim_end_matches(['\n', '\r']);
let (directive, value) = rest.split_once('=')?;
let id_field = value.strip_prefix('<')?.strip_prefix("ID=")?;
let id = id_field.split([',', '>']).next()?;
Some((directive.to_string(), id.to_string()))
}

/// Extract the GT field as a string for a specific sample from a VCF record.
///
/// Returns `None` if the GT field is absent for this sample.
Expand Down Expand Up @@ -376,4 +477,83 @@ chr1\t10\t.\tA\tT\t.\t.\t.\tGT\t0|0|1
assert_eq!(parsed.sample_ploidy, 3);
assert_eq!(parsed.by_contig.get("chr1").map(Vec::len), Some(1));
}

/// A VCF that redeclares an `##INFO` ID must still parse. noodles' header
/// parser rejects duplicate IDs outright (`duplicate INFO ID: ...`), but
/// real-world VCFs from tools like Manta/DELLY carry repeated definitions
/// that bcftools tolerates; holodeck must not be stricter than bcftools.
#[test]
fn parse_variants_by_contig_tolerates_duplicate_info_headers() {
let vcf = "\
##fileformat=VCFv4.3
##contig=<ID=chr1,length=100>
##INFO=<ID=BREAKSIMLENGTH,Number=1,Type=Integer,Description=\"first\">
##INFO=<ID=BREAKSIMLENGTH,Number=1,Type=Integer,Description=\"second\">
##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSAMPLE
chr1\t10\t.\tA\tT\t.\tPASS\tBREAKSIMLENGTH=5\tGT\t0/1
";
let f = write_temp_vcf(vcf);
let dict = dict_for(&[("chr1", 100)]);
let parsed = parse_variants_by_contig(f.path(), None, &dict).unwrap();
assert_eq!(parsed.by_contig.get("chr1").map(Vec::len), Some(1));
assert_eq!(parsed.by_contig["chr1"][0].position, 9);
}

/// Duplicate `##FORMAT`/`##FILTER`/`##contig` declarations are tolerated
/// alongside `##INFO`; all of noodles' duplicate-ID categories are covered
/// by the same dedup pass.
#[test]
fn parse_variants_by_contig_tolerates_duplicate_format_filter_contig() {
let vcf = "\
##fileformat=VCFv4.3
##contig=<ID=chr1,length=100>
##contig=<ID=chr1,length=100>
##FILTER=<ID=q10,Description=\"first\">
##FILTER=<ID=q10,Description=\"second\">
##FORMAT=<ID=GT,Number=1,Type=String,Description=\"first\">
##FORMAT=<ID=GT,Number=1,Type=String,Description=\"second\">
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSAMPLE
chr1\t10\t.\tA\tT\t.\tPASS\t.\tGT\t0/1
";
let f = write_temp_vcf(vcf);
let dict = dict_for(&[("chr1", 100)]);
let parsed = parse_variants_by_contig(f.path(), None, &dict).unwrap();
assert_eq!(parsed.by_contig.get("chr1").map(Vec::len), Some(1));
}

/// `validate_vcf_sample` reads the header through the same lenient path, so
/// it too must accept a VCF with duplicate definitions.
#[test]
fn validate_vcf_sample_tolerates_duplicate_info_headers() {
let vcf = "\
##fileformat=VCFv4.3
##contig=<ID=chr1,length=100>
##INFO=<ID=DP,Number=1,Type=Integer,Description=\"first\">
##INFO=<ID=DP,Number=1,Type=Integer,Description=\"second\">
##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tNA12878
chr1\t10\t.\tA\tT\t.\tPASS\tDP=5\tGT\t0/1
";
let f = write_temp_vcf(vcf);
let sample = validate_vcf_sample(f.path(), None).unwrap();
assert_eq!(sample, "NA12878");
}

#[test]
fn structured_header_key_extracts_directive_and_id() {
assert_eq!(
structured_header_key("##INFO=<ID=DP,Number=1,Type=Integer>\n"),
Some(("INFO".to_string(), "DP".to_string()))
);
assert_eq!(
structured_header_key("##contig=<ID=chr1,length=100>"),
Some(("contig".to_string(), "chr1".to_string()))
);
// No ID field -> not a dedup candidate.
assert_eq!(structured_header_key("##fileformat=VCFv4.3"), None);
assert_eq!(structured_header_key("##source=holodeck-mutate"), None);
// ID not the first field -> conservatively not deduplicated.
assert_eq!(structured_header_key("##INFO=<Number=1,ID=DP>"), None);
}
}
Loading