VCF write: 15+ min stalls on individual CDS exons (1KG chr2 MANE Select, gvl 0.24.1)

[bschilder]

Summary

gvl.write(...) from VCF input on chr2 of the 1KG NYGC 30x phased panel makes forward progress at the start (~40 region/s on the easy regions) but then enters a regime where single CDS exons (43-200 bp) take 15+ minutes each to process. Across two runs (different max_mem budgets, different host) the slowdown reproduces in the same chr2 region range. The new Variant reader resident size: 0 B; ... log line introduced in #163 suggests the index-aware accounting is reporting zero footprint for the VCF reader, so the budgeting fix doesn't apply to VCF inputs in practice.

Environment

genvarloader	0.24.1 (head of #163, commit 2d7d660)
genoray	>=2.4.0,<3 (resolved via PR's pin)
Python	3.11
OS	Ubuntu 24.04 (Runpod runpod/pytorch:1.0.2-cu1281-torch280-ubuntu2404)
Host	16 vCPU / 250 GB RAM (A100-SXM4-80GB; GPU unused — gvl.write is CPU-bound)
Input VCF	1KG NYGC 30x phased panel, chr2 — 1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz (~1 GB compressed, ~5M variants, 3202 samples)
Input BED	MANE Select chr2 CDS (14,818 rows, from Ensembl 113 GTF, ordered by chrom+start)
max_mem	tried both "2g" and "16g"; behaviour identical

Reproduction

# pip install "genvarloader @ git+https://github.com/mcvickerlab/GenVarLoader@2d7d660" "genoray>=2.4.0"
import polars as pl
import pooch
import genvarloader as gvl
from pathlib import Path

# 1. VCF (1KG NYGC 30x phased panel — public)
vcf_url = (
    "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/"
    "1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/"
    "1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
)
vcf_path = pooch.retrieve(vcf_url, known_hash=None, path="/tmp/gvl-repro")
pooch.retrieve(vcf_url + ".tbi", known_hash=None, path="/tmp/gvl-repro")

# 2. MANE Select CDS BED for chr2 (Ensembl 113)
gtf_url = "https://ftp.ensembl.org/pub/release-113/gtf/homo_sapiens/Homo_sapiens.GRCh38.113.chr.gtf.gz"
gtf_path = pooch.retrieve(gtf_url, known_hash=None, path="/tmp/gvl-repro")
gtf = pl.read_csv(
    gtf_path, separator="\t", comment_prefix="#", has_header=False,
    new_columns=["chrom","source","feature","start","end","score","strand","frame","attrs"],
    infer_schema_length=0,
).with_columns(pl.col("start").cast(pl.Int64), pl.col("end").cast(pl.Int64))
gtf = gtf.filter(pl.col("feature") == "CDS")
gtf = gtf.with_columns(
    pl.col("attrs").str.extract(r'transcript_id "([^"]+)"').alias("transcript_id"),
    pl.col("attrs").str.extract(r'exon_number "?([0-9]+)"?').cast(pl.Int64, strict=False).alias("exon_number"),
    pl.col("attrs").str.contains("MANE_Select").alias("is_mane_select"),
)
mane = gtf.filter(pl.col("is_mane_select") & (pl.col("chrom") == "2"))
bed = mane.select(
    ("chr" + pl.col("chrom").cast(pl.Utf8)).alias("chrom"),
    (pl.col("start") - 1).alias("chromStart"),
    pl.col("end").alias("chromEnd"),
    "transcript_id", "exon_number", "strand",
)
print(f"chr2 MANE Select CDS rows: {bed.height}")   # 14818

# 3. write — this is where it stalls
gvl.write(
    path="/tmp/gvl-repro/chr2.svar",
    bed=bed,
    variants=str(vcf_path),
    overwrite=True,
    max_mem="16g",
)

Observed behaviour

The tqdm progress bar shows a clean fast start, a slow patch around region 2004-2078 (rate dropped from ~47 region/s to ~1.3 region/s, then recovered to ~14 region/s by region 2284), and then a 15+ minute stall on a single 43-base exon at region 2683:

gvl 0.24.1
[chr2] gvl.write(/<...>/chr2, bed=14818 regions, vcf=1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz)
INFO     | genvarloader._dataset._write:write:126 - Writing dataset to /<...>/chr2
INFO     | genvarloader._dataset._write:write:198 - Using 3202 samples.
INFO     | genvarloader._dataset._write:write:203 - Writing genotypes.
INFO     | genvarloader._dataset._write:write:207 - Variant reader resident size: 0 B; max_mem budget: 16.00 GiB; available for chunking: 16.00 GiB
INFO     | genvarloader._dataset._write:_write_from_vcf:363 - VCF genoray index is invalid, writing
INFO     | genoray._vcf:_load_index:1096 - Loading genoray index.
WARNING  | genvarloader._dataset._write:_write_from_vcf:438 - A region has no variants for any sample. ...

Processing genotypes for 14818 regions on contig chr2:   0%|          | 0/14818 [00:00<?, ? region/s]
Processing genotypes for 14818 regions on contig chr2:   3%|▎         | 405/14818 [00:10, 43.26 region/s]
Processing genotypes for 14818 regions on contig chr2:   8%|▊         | 1114/14818 [00:30, 38.08 region/s]
Processing genotypes for 14818 regions on contig chr2:  13%|█▎        | 1941/14818 [00:45, 47.73 region/s]   ← cruising
Processing genotypes for 14818 regions on contig chr2:  14%|█▎        | 2004/14818 [00:59, 47.73 region/s]
Processing genotypes for 14818 regions on contig chr2:  14%|█▎        | 2042/14818 [02:11<1:14:47,  2.85 region/s]  ← first slowdown
Processing genotypes for 14818 regions on contig chr2:  14%|█▍        | 2078/14818 [03:00, 1.33 region/s]   ← bottom of first slowdown
Processing genotypes for 14818 regions on contig chr2:  15%|█▌        | 2284/14818 [03:05,  7.04 region/s]  ← recovered
Processing genotypes for 14818 regions on contig chr2:  17%|█▋        | 2544/14818 [03:10, 14.68 region/s]
Processing genotypes for 14818 regions on contig chr2:  18%|█▊        | 2682/14818 [03:29<13:46, 14.68 region/s]   ← was here at t=03:29
Processing genotypes for 14818 regions on contig chr2:  18%|█▊        | 2683/14818 [19:58<7:51:19,  2.33s/ region] ← 16:29 min spent on ONE 43-bp exon

Region 2683 is a CDS exon of E2F6 at chr2:11,446,480-11,446,523 (43 bp, strand -) — structurally unremarkable. The first slow patch (regions 2004-2078) similarly hits common ordinary genes: MYL1, ANAPC1, SGPP2, CRACDL, TANC1. None of these are textbook hard cases (HLA, KIR, immunoglobulin, low-complexity / repetitive regions).

Hypothesis

The slowness doesn't correlate with variant density of the BED region nor with chunk boundaries (a chunk flush would be deterministic, not random). My current guess is some kind of per-region quadratic scan or cache-miss pattern in the VCF reader: each region triggers a scan over a large data structure whose size grows with cumulative regions processed.

The "0 B" reported by the new #163 accounting log strongly suggests the VCF reader path is not surfacing its memory footprint via genoray's .nbytes — so #163's budgeting can't influence chunk-vs-index trade-off for VCF input.

Questions

1. Is Variant reader resident size: 0 B expected for VCF inputs, or is the genoray VCF reader missing the .nbytes property?
2. Have you observed per-region runtime variation of this magnitude (40× region/s → 0.4 region/s) on other VCF inputs? Could you profile a single-region call on the E2F6 exon above against the 1KG NYGC 30x chr2 VCF?
3. Would converting the VCF to PGen up front (so the PGen .nbytes is populated) likely sidestep this?

Happy to attach a private slice of the in-pod log (~6 KB after stripping internal paths) or run additional instrumentation if helpful.

---

(Reporting from external pipeline use of gvl.write; not opening as a regression on #163 since the underlying behaviour is consistent across max_mem values — #163 made the budget accounting visible but didn't introduce this.)

[bschilder]

Update: switched to gvl.get_splice_bed, stall still reproduces

Two corrections + a fresh data point.

First, my reprex was using the API wrong. The original repro hand-rolled a CDS BED from the Ensembl GTF (preserving GTF row order, which lists - strand CDS exons in 5'→3' transcript order = decreasing chromStart). That's not what gvl.write expects. The correct call is gvl.get_splice_bed(gtf, contigs=[...]), which sorts ascending by chromStart and applies the GVL-recommended filters (transcript_support_level='1', require_multiple_of_3=True). Apologies for the noise.

Second, the stall pattern still reproduces with the corrected pipeline. Same VCF, same host class, max_mem="16g", same GVL 0.24.1 at #163 head — only the BED-build call changed:

import genvarloader as gvl
splice_bed = gvl.get_splice_bed(gtf_path, contigs=["2"])
# (chrom column is bare "2"; my pipeline prefixes "chr" to match the VCF)
splice_bed = splice_bed.with_columns(("chr" + pl.col("chrom").cast(pl.Utf8)).alias("chrom"))
chrom_bed = splice_bed.filter(pl.col("chrom") == "chr2")
gvl.write(path="/tmp/out/chr2", bed=chrom_bed, variants=str(vcf_path), overwrite=True, max_mem="16g")

get_splice_bed returned 11,497 chr2 CDS rows (vs my 14,818 from the broader MANE Select filter). Throughput trajectory looks identical to the original report: a fast burst (regions 0→2028 in 49s, ~52 region/s), then tqdm goes silent for 10+ minutes with no new progress lines emitted. Pod is alive and stdout-buffering settings (PYTHONUNBUFFERED=1, TQDM_MININTERVAL=5) are correct.

Same Variant reader resident size: 0 B log line for VCF input, same chunk budget visible to GVL.

New log (sanitized, full output up to silence)

gvl 0.24.1
[chr2] gvl.write(/tmp/out/chr2, bed=11497 regions, vcf=1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz)
INFO     | genvarloader._dataset._write:write:126 - Writing dataset to /tmp/out/chr2
INFO     | genvarloader._dataset._write:write:198 - Using 3202 samples.
INFO     | genvarloader._dataset._write:write:203 - Writing genotypes.
INFO     | genvarloader._dataset._write:write:207 - Variant reader resident size: 0 B; max_mem budget: 16.00 GiB; available for chunking: 16.00 GiB
INFO     | genvarloader._dataset._write:_write_from_vcf:363 - VCF genoray index is invalid, writing
INFO     | genoray._vcf:_load_index:1096 - Loading genoray index.    ← ~3m 38s reindex
27541rows [00:00, 71798.20rows/s]
Processing genotypes for 11497 regions on contig chr2:   0%|          | 0/11497 [00:00<?, ? region/s]
WARNING  | genvarloader._dataset._write:_write_from_vcf:438 - A region has no variants for any sample. ...

Processing genotypes for 11497 regions on contig chr2:   1%|          | 131/11497 [00:05<07:14, 26.14 region/s]
Processing genotypes for 11497 regions on contig chr2:   3%|▎         | 394/11497 [00:10<04:26, 41.61 region/s]
Processing genotypes for 11497 regions on contig chr2:   6%|▌         | 693/11497 [00:15<03:36, 49.86 region/s]
Processing genotypes for 11497 regions on contig chr2:   9%|▊         | 1000/11497 [00:20<03:13, 54.36 region/s]
Processing genotypes for 11497 regions on contig chr2:  12%|█▏        | 1334/11497 [00:25<02:52, 58.81 region/s]
Processing genotypes for 11497 regions on contig chr2:  14%|█▍        | 1629/11497 [00:31<03:09, 52.10 region/s]
Processing genotypes for 11497 regions on contig chr2:  16%|█▋        | 1896/11497 [00:37<03:03, 52.19 region/s]
Processing genotypes for 11497 regions on contig chr2:  18%|█▊        | 2028/11497 [00:49<03:01, 52.19 region/s]
   ← last tqdm line, then 10+ min of silence; pod alive, heartbeat re-uploading same content

So the early throughput at fix-time is exactly the same as in the original report (the buggy hand-rolled BED also raced past the first ~2000 regions in ~49s on the easy portion of chr2). What follows is silent for many minutes regardless of BED ordering.

Suspicion (revised)

The first ~2000 regions are processed fast, then a long quiet period begins. Candidates:

1. max_mem-driven chunk flush. With max_mem="16g", a chunk fills at some region count and tqdm is silent during the flush (no per-region updates while ~16 GiB is serialized to disk). If that's the explanation, it would be enormously helpful for gvl.write to emit a single "chunk N flushed (M bytes in T seconds)" line per flush so users can distinguish "writing to disk" from "stuck."
2. Per-region quadratic in the VCF→genotype path. The original 0.24 introduced index-aware accounting; would expect a quadratic to be visible as monotonically slowing region/s rather than a hard cliff — but if it's a once-per-chunk-buffer fill, it would look like a cliff.
3. Both — fast burst, then chunk flush + quadratic accumulating each subsequent chunk.

Reproducible reprex (corrected)

# pip install "genvarloader @ git+https://github.com/mcvickerlab/GenVarLoader@2d7d660" "genoray>=2.4.0" polars pooch
import polars as pl
import pooch
import genvarloader as gvl

vcf_url = (
    "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/"
    "1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/"
    "1kGP_high_coverage_Illumina.chr2.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
)
vcf_path = pooch.retrieve(vcf_url, known_hash=None, path="/tmp/gvl-repro")
pooch.retrieve(vcf_url + ".tbi", known_hash=None, path="/tmp/gvl-repro")

gtf_url = "https://ftp.ensembl.org/pub/release-113/gtf/homo_sapiens/Homo_sapiens.GRCh38.113.chr.gtf.gz"
gtf_path = pooch.retrieve(gtf_url, known_hash=None, path="/tmp/gvl-repro")

splice_bed = gvl.get_splice_bed(str(gtf_path), contigs=["2"])
splice_bed = splice_bed.with_columns(("chr" + pl.col("chrom").cast(pl.Utf8)).alias("chrom"))
print(f"chr2 rows: {splice_bed.height}")   # 11497

gvl.write(
    path="/tmp/gvl-repro/chr2.svar",
    bed=splice_bed,
    variants=str(vcf_path),
    overwrite=True,
    max_mem="16g",
)

Expect: a fast initial pass (0 → ~2028 regions, ~52 region/s, ~50 s), then tqdm silence for 10+ min on a host with 16 vCPU / 250 GB.

Happy to try max_mem="2g" / max_mem="64g" next as a quick probe of the chunk-flush hypothesis, or to test PGen-converted input if you have a preferred conversion path. Leaving the issue open.

[bschilder]

Big correction: our reprex was using gvl.write wrong

I owe an apology — the reprex above was passing a raw VCF path to gvl.write(... variants=vcf_path), which routes through _write_from_vcf. The correct workflow (per the maintainer / API design) is:

1. One-time: convert VCF → SVAR via genoray.SparseVar.from_vcf(out="<chrom>.svar", vcf=VCF(vcf_path), max_mem="32g"). The output is a real SparseVar directory.
2. Then: gvl.write(path="<chrom>.gvl", bed=bed, variants=SparseVar("<chrom>.svar")) — this goes through _write_from_svar, which is the supported/recommended path.

The _write_from_vcf shortcut path is where the chr2 ~18% stall we were chasing actually lives. Switching our downstream pipeline to the two-step (SparseVar.from_vcf + gvl.write(variants=SparseVar)) is the right fix; no GVL bug to track here, this was user error.

I also realized PR #165 (the deduplicate_overlapping_cds parameter I proposed) is biologically wrong — exon sharing across transcripts of the same gene is normal alternative splicing, not duplication to drop. Closed.

You can probably close this issue. Apologies for the noise; thanks for the patient review.

[bschilder]

Resolved — fully end-to-end-verified

Closing this with a complete write-up so the thread isn't open with stale conclusions.

Root cause: gvl.write(variants=<raw VCF>) shortcut + un-normalized VCFs

Our pipeline was calling gvl.write(path, bed, variants=<raw_vcf_path>), which routes through GVL's _write_from_vcf shortcut. Two compounding issues:

1. Skipping the VCF→SparseVar split that the GVL docs recommend. The supported route is genoray.SparseVar.from_vcf(...) → gvl.write(variants=SparseVar(...)), which goes through _write_from_svar. We were taking the shortcut.
2. Un-normalized VCFs. The 1KG NYGC 30x phased panel ships without the bcftools-norm preprocessing GVL requires (left-aligned + bi-allelic + atomic). Feeding _write_from_vcf the raw VCF triggered a slow path on chr2 specifically.

Both PR #163 (max_mem accounting) and #165 (a deduplicate_overlapping_cds flag I proposed and then withdrew as biologically incorrect) were red herrings — neither addressed the actual issue.

Fix in our downstream pipeline (no GVL change needed)

Switched to the proper two-step workflow + bcftools-norm preprocessing:

# 1. One-time per cohort × chrom — bcftools norm with GVL-recommended flags
bcftools norm -f <ref> -a --atom-overlaps . -m -any --multi-overlaps . -O z -o <out>

# 2. One-time per cohort × chrom — build the SparseVar
SparseVar.from_vcf(out=<chrom>.svar, vcf=VCF(<normed>), max_mem="32g")

# 3. Per panel — build the GVL dataset over a BED
gvl.write(path=<chrom>.gvl, bed=<bed>, variants=SparseVar(<chrom>.svar), max_mem="2g")

End-to-end verification

Full rebuild of the 1KG NYGC 30x panel across all 22 autosomes (3202 samples, gvl.get_splice_bed(gtf, transcript_support_level=None) = ~1.9k–~79k CDS rows per chrom) on 22 parallel A100 pods. chr2 — the previously-cursed chromosome — landed cleanly:

• chr2 SVAR: 7.49 GB on R2
• chr2 GVL: 6.69 GB on R2
• chr2 gvl.write step: 153 seconds at ~50 region/s sustained, no per-region stall
• Total fleet: 88.3 GB SparseVar + 75.4 GB GVL on R2 in ~45 min wall clock

build_gvl over a pre-built SparseVar runs at ~400–950 region/s depending on chrom size. The original _write_from_vcf path was hitting < 1 region/s on chr2 at the stall point. ~1000× throughput improvement.

Closing this issue

Not a GVL bug — user-pipeline issue we resolved by going through the documented _write_from_svar path. Thanks for the review on the earlier (incorrect) PRs.

If it's helpful, two things that would make the trap harder to fall into:

1. A docstring tweak on gvl.write calling out that passing a raw VCF goes through the slower _write_from_vcf path and recommending pre-conversion to SparseVar for any production workload.
2. A gvl.normalize_vcf(...) (or similar) helper that wraps the bcftools-norm one-liner from the write doc page, so callers can do vcf_normed = gvl.normalize_vcf(vcf, ref) instead of needing to shell out.

Happy to PR either if useful.