Paper

Summary

Variable number tandem repeats (VNTRs) are genomic regions in which a sequence motif is repeated in tandem a variable number of times across individuals. Large VNTRs — those with repeat units spanning hundreds to thousands of base pairs — are a substantial source of structural variation in the human genome and have been linked to complex traits and disease [1]. Because their repeat units far exceed the length of standard short sequencing reads, large VNTRs cannot be genotyped by conventional variant-calling approaches, and haplotype-resolved copy number estimation from short-read whole-genome sequencing (WGS) remains a methodological gap. Mukamel et al. [1] identified 734 protein-coding VNTRs with elevated identity-by-descent (IBD) rates across the human genome, many of which are likely functional but lack accessible copy number estimation tools for cohort-scale short-read data.

A particularly well-characterized example is the Kringle IV type-2 (KIV-2) VNTR within the LPA gene, which is the primary genetic determinant of Lipoprotein(a) [Lp(a)] concentration [1, 2]. Lp(a) is a causal, independent, and highly heritable (\(\text{h}^2 >70\%\)) cardiovascular disease (CVD) risk factor elevated in an estimated 1.5 billion people globally, with concentrations differing substantially across ancestry groups [3, 4, 5]. The KIV-2 repeat unit spans ~5.5 kb, individuals carry 1–40 copies per haplotype, and copy number is inversely correlated with Lp(a) concentration, making accurate haplotype-resolved KIV-2 CNV estimation central to ancestry-aware polygenic risk scores and causal inference studies [1, 2, 6, 7].

GRiD (Genomic Repeat inference from Depth) is an open-source Python pipeline for haplotype-resolved large VNTR copy number estimation from short-read WGS data. GRiD combines read-depth normalization, read-depth profile-matched nearest-neighbor estimation, and IBD-based haplotype inference [8] into a single installable, multi-threaded, YAML-configured workflow. GRiD was developed and validated for the LPA KIV-2 locus, but its core framework is locus-agnostic: any large VNTR in a cohort with sufficient sample size can be targeted by supplying the relevant genomic coordinates in the configuration file. The pipeline is designed to be accessible to researchers without extensive bioinformatics expertise.

Statement of Need

Large VNTRs present a shared set of computational challenges that place them outside the scope of existing repeat genotyping software. Their repeat units are orders of magnitude larger than the short tandem repeats (STRs, 1–6 bp motifs) targeted by most tools, and in high-copy-number individuals the VNTR region can span hundreds of kilobases, making read-pair spanning approaches infeasible with standard 150 bp paired-end reads. Haplotype resolution adds a further challenge: phasing cannot rely on long reads — which remain unavailable for most large epidemiological cohorts — or on pre-existing phased variant data at the locus. GRiD is designed to address these challenges directly through three features: (1) a read-depth normalization and nearest-neighbor framework that is locus-agnostic and parameterized entirely through a configuration file, (2) an IBD-based iterative phasing algorithm [8] that decomposes diploid copy number into haplotype-specific estimates using cohort relatedness alone, and (3) a single installable pipeline that encapsulates the full workflow from short-read WGS to haplotype-resolved copy number in a reproducible, multi-threaded form.

For the specific case of LPA KIV-2 — GRiD’s primary development and validation target — several tools have addressed parts of this problem, but none provides the complete short-read-to-haplotype-resolved workflow in open-source, cohort-scalable form. KILDA [9] estimates diploid KIV-2 copy number from FASTQ files using an alignment-free k-mer approach but does not resolve haplotype-specific copy numbers. The Illumina DRAGEN LPA Caller provides KIV-2 estimation but requires the proprietary DRAGEN platform. vntr-calling-nf [10] offers a Nextflow-based pipeline for VNTR variant calling but targets coding variant detection rather than copy number quantification and has been primarily validated on whole-exome sequencing data. As a result, researchers studying Lp(a) at cohort scale with short-read WGS have largely relied on in-house scripts or unpublished internal pipelines, creating barriers to reproducibility and cross-cohort comparability.

State of the Field

Several mature tools exist for genotyping tandem repeats from short-read sequencing data. ExpansionHunter [11] and GangSTR [12] use probabilistic models of read-pair orientation and coverage to estimate repeat lengths, but are designed for STRs and short VNTRs where reads can partially or fully span the repeat unit. HipSTR [13] performs haplotype-level STR genotyping using phased reads but is similarly constrained to short repeat units accessible by individual reads. Tools such as TRTools [14] provide downstream analysis utilities for STR calls but do not extend to large-VNTR loci. adVNTR [15] uses hidden Markov models to genotype VNTRs genome-wide and supports both Illumina and PacBio data, but is designed for VNTRs of up to a few hundred base pairs where reads can overlap the repeat boundary, which is qualitatively different from the 5.6 kb KIV-2 repeat unit. More recently, long-read platforms have enabled direct assembly of complex repeat regions [16], and Kivvi [17] specifically targets KIV-2 from PacBio HiFi data; however, long-read data are not yet available for most large epidemiological cohorts, which remain dominated by short-read WGS.

For the LPA KIV-2 locus, KILDA [9] provides an alignment-free diploid copy number estimate from short-read FASTQ files using a k-mer counting strategy calibrated against long-read assemblies, demonstrating strong concordance with assembly-based diploid estimates. However, no existing tool addresses haplotype-level phasing for KIV-2 or for large VNTRs more broadly from short-read data. The foundational algorithmic framework underlying GRiD was established by Mukamel et al. [1], who demonstrated that population-normalized read depth is the most tractable approach for large VNTR copy number estimation from short-read cohort data, and by Hujoel et al. [8], who developed an IBD-based iterative algorithm for decomposing diploid copy number into haplotype-specific contributions without requiring long reads or pre-existing phased variant data. GRiD integrates and operationalizes both methodological contributions into a single accessible software package with read-depth profile-aware nearest-neighbor normalization, a full configuration system, and multi-cohort scalability.

Software Design

GRiD is implemented in Python (≥3.8) and follows a modular architecture organized into seven sequential pipeline steps, each implemented as an independent utility module under grid/utils/. The pipeline is orchestrated by grid/pipeline.py and driven by a user-supplied YAML configuration file, which specifies file paths, genomic coordinates, per-step parameters, and which steps to execute. The command-line interface is implemented with Click [18] and provides an entry point: grid wgs for whole-genome sequencing data.

The Whole Genome Sequencing pipeline proceeds as follows:

  1. Index verification / creation: Ensures CRAM/BAM index files (.crai/.bai) exist for all input samples, creating them with samtools [19] if needed.

  2. Read counting: Counts properly paired reads in the target VNTR region using pysam [20], filtering by mapping quality and SAM flag.

  3. Coverage estimation: Runs mosdepth [21] across the genome in binned mode to produce per-sample depth profiles.

  4. Coverage normalization: Normalizes coverage within individuals (by sample mean depth) and across individuals (by a z-score transformation), then filters to high-variance regions enriched for VNTR signal.

  5. Nearest-neighbor identification: Uses scikit-learn [22] to compute Euclidean distances in normalized depth space and identify the top-N genomically similar neighbors per individual, providing a cohort-matched reference for CNV normalization.

  6. Diploid copy number estimation: Computes per-exon diploid copy number for each individual by normalizing their read counts relative to their nearest neighbors’ read counts, accounting for sample-specific sequencing depth.

  7. Haplotype inference: Applies an IBD-based iterative phasing algorithm [8] that uses haplotype-matched neighbors to decompose diploid copy number estimates into haplotype-specific contributions.

All multi-sample steps are parallelized using Python’s concurrent.futures.ThreadPoolExecutor, and progress is reported via the Rich terminal library [23]. Configuration validation at pipeline startup provides early, informative error messages before any compute-intensive steps run.

Research Impact Statement

GRiD has been validated across multiple large-scale, multi-ancestry whole-genome sequencing cohorts representing tens of thousands of individuals spanning Hispanic/Latino, African American, and European American ancestry groups. Its nearest-neighbor normalization approach naturally accounts for ancestry and sequencing batch effects without requiring explicit ancestry labels, making it well-suited for population-diverse studies where CNV tools calibrated on European-ancestry reference panels may systematically underperform.

GRiD has attracted national and international interest from research groups across multiple conferences, with inquiries from groups seeking to apply large VNTR copy number estimation to their own short-read WGS cohorts. It is designed for accessibility, with comprehensive documentation of configuration, outputs, and algorithmic details to lower the barrier to entry for epidemiologists and clinicians without extensive bioinformatics training.

Mathematics

Coverage estimation. Weighted mean depth across the VNTR region from mosdepth bins is computed as:

\[ \text{Coverage} = \left\lfloor 100 \cdot \frac{ \sum_{i \in R} \left( \bar{C}_i \cdot \bigl( \min(e, r_{i,e}) - \max(s, r_{i,s}) \bigr) \right) }{ \sum_{i \in R} \bigl( \min(e, r_{i,e}) - \max(s, r_{i,s}) \bigr) } \right\rceil \]

where \(R\) is the set of bins overlapping the region \([s, e]\), \(\bar{C}_i\) is the mean depth of bin \(i\), and \(r_{i,s}\), \(r_{i,e}\) are the bin boundaries.

Normalization. Let \(D \in \mathbb{R}^{N \times M}\) be the depth matrix for \(N\) individuals across \(M\) genomic bins. Within-individual normalization gives \(D^{(1)}_{ij} = D_{ij} / \bar{D}_i\), where \(\bar{D}_i\) is the individual mean. Across-individual normalization then yields:

\[ D^{\mathrm{norm}}_{ij} = \frac{D^{(1)}_{ij} - \mu_j}{\sqrt{\mu_j}} \]

where \(\mu_j\) is the population mean of \(D^{(1)}_{\cdot j}\). Regions are filtered by a variance ratio \(\sigma_j^2 / \mu_j\), and the top fraction of high-variance regions are retained for neighbor computation.

Neighbor-normalized diploid copy number. For individual \(i\) with read count \(r_i\) and depth scale \(s_i\), the diploid copy number estimate is:

\[ \widehat{\text{dipCN}}_i = \frac{r_i / s_i}{\frac{1}{|N_i|} \sum_{j \in N_i} r_j / s_j} \]

where \(N_i\) is the set of nearest neighbors identified from the normalized depth space.

Haplotype inference. Following Hujoel et al. [8], haplotype-specific copy numbers \(h_{i,1}\) and \(h_{i,2}\) are iteratively estimated by updating each haplotype’s value proportionally to the mean copy number of IBD-matched haplotype neighbors:

\[ h_{i,k}^{(t+1)} = \widehat{\text{dipCN}}_i \cdot \frac{\bar{h}_{N_{i,k}}^{(t)}}{\bar{h}_{N_{i,1}}^{(t)} + \bar{h}_{N_{i,2}}^{(t)}} \]

where \(\bar{h}_{N_{i,k}}^{(t)}\) is the mean haplotype copy number of IBD neighbors on haplotype \(k\) at iteration \(t\).

AI Usage Disclosure

GitHub Copilot and Anthropic Claude were used to assist with code development and documentation writing during the preparation of this software and manuscript. All AI-generated content was reviewed line-by-line by the authors, with corrections and edits applied as necessary to ensure accuracy, correctness, and consistency with the underlying methodology. No AI-generated code or text was incorporated without direct author verification.

Acknowledgements

This work was supported by [GRANT INFORMATION].

References

[1] (1,2,3,4,5)

Ronen E. Mukamel, Robert E. Handsaker, Maxwell A. Sherman, and others. Protein-coding repeat polymorphisms strongly shape diverse human phenotypes. Science, 373(6562):1499–1505, 2021. doi:10.1126/science.abg8635.

[2] (1,2)

Gerd Utermann, Hans J. Menzel, Hans G. Kraft, Hans C. Duba, Hans G. Kemmler, and Cornelia Seitz. Lp(a) glycoprotein phenotypes: inheritance and relation to Lp(a)-lipoprotein concentrations in plasma. Journal of Clinical Investigation, 80(2):458–465, 1987. doi:10.1172/JCI113093.

[3]

Sotirios Tsimikas. A test in context: lipoprotein(a): diagnosis, prognosis, controversies, and emerging therapies. Journal of the American College of Cardiology, 69(6):692–711, 2017. doi:10.1016/j.jacc.2016.11.042.

[4]

Florian Kronenberg, Serena Mora, Erik S. G. Stroes, and others. Lipoprotein(a) in atherosclerotic cardiovascular disease and aortic stenosis: a European Atherosclerosis Society consensus statement. European Heart Journal, 43(39):3925–3946, 2022. doi:10.1093/eurheartj/ehac361.

[5]

Børge G. Nordestgaard, M. John Chapman, Kausik Ray, Jan Borén, Felicita Andreotti, Gerald F. Watts, Henry N. Ginsberg, Pierre Amarenco, Alberico Catapano, Olivier S. Descamps, Edward Fisher, Petri T. Kovanen, Jan Albert Kuivenhoven, Philippe Lesnik, Luis Masana, Zeljko Reiner, Marja-Riitta Taskinen, Lale Tokgözoğlu, and Anne Tybjærg-Hansen. Lipoprotein(a) as a cardiovascular risk factor: current status. European Heart Journal, 31(23):2844–2853, 2010. doi:10.1093/eurheartj/ehq386.

[6]

Robert Clarke, John F. Peden, Jemma C. Hopewell, and others. Genetic variants associated with Lp(a) lipoprotein level and coronary disease. New England Journal of Medicine, 361(26):2518–2528, 2009. doi:10.1056/NEJMoa0902604.

[7]

Pia R. Kamstrup, Anne Tybjærg-Hansen, Rolf Steffensen, and Børge G. Nordestgaard. Genetically elevated lipoprotein(a) and increased risk of myocardial infarction. JAMA, 301(22):2331–2339, 2009. doi:10.1001/jama.2009.801.

[8] (1,2,3,4,5)

Margaux LA Hujoel, Robert E Handsaker, David Tang, Nolan Kamitaki, Ronen E Mukamel, Simone Rubinacci, Pier Francesco Palamara, Steven A McCarroll, and Po-Ru Loh. Insights into dna repeat expansions among 900,000 biobank participants. Nature, pages 1–10, 2026.

[9] (1,2)

Lydia Molitor and others. KILDA: k-mer-based inference of LPA diploid copy number from short-read alignment-free data. NAR Genomics and Bioinformatics, 2025. doi:10.1093/nargab/lqaf070.

[10]

genepi. Vntr-calling-nf: a Nextflow pipeline for VNTR variant calling. 2024. URL: https://github.com/genepi/vntr-calling-nf.

[11]

Egor Dolzhenko and others. ExpansionHunter: a sequence-graph-based tool to analyze variation in short tandem repeat regions. Bioinformatics, 35(22):4754–4756, 2019. doi:10.1093/bioinformatics/btz431.

[12]

Nima Mousavi, Sharona Shleizer-Burko, Richard Yanicky, and Melissa Gymrek. GangSTR: template-based estimation of short tandem repeats. Nucleic Acids Research, 47(15):e90, 2019. doi:10.1093/nar/gkz501.

[13]

Thomas Willems, Or Zuk, Yaniv Erlich, and Melissa Gymrek. Genome-wide profiling of heritable and de novo STR variations. Nature Methods, 14:590–592, 2017. doi:10.1038/nmeth.4267.

[14]

Nima Mousavi, Jonathan Margoliash, Neha Pusarla, and others. TRTools: a toolkit for genome-wide analysis of tandem repeats. Bioinformatics, 37(18):3062–3064, 2021. doi:10.1093/bioinformatics/btab358.

[15]

Mehrdad Bakhtiari, Sharona Shleizer-Burko, Melissa Gymrek, Vikas Bansal, and Vineet Bafna. Targeted genotyping of variable number tandem repeats with advntr. Genome Research, 28(11):1709–1719, 2018. doi:10.1101/gr.235119.118.

[16]

Katja Didericksen and others. Long-read sequencing resolves complex repeat structures in the human genome. Nature Genetics, 2024. doi:VERIFY.

[17]

PacificBiosciences. Kivvi: KIV-2 and D4Z4 copy number estimation from PacBio HiFi whole-genome sequencing. 2026. URL: https://github.com/PacificBiosciences/Kivvi.

[18]

Armin Ronacher and others. Click: python composable command line interface toolkit. 2024. URL: https://palletsprojects.com/p/click/.

[19]

Heng Li, Bob Handsaker, Alec Wysoker, and others. The sequence alignment/map format and SAMtools. Bioinformatics, 25(16):2078–2079, 2009. doi:10.1093/bioinformatics/btp352.

[20]

Andreas Heger and others. Pysam: a python module for reading and manipulating sam/bam/cram files. 2024. URL: https://github.com/pysam-developers/pysam.

[21]

Brent S. Pedersen and Aaron R. Quinlan. Mosdepth: quick coverage calculation for genomes and exomes. Bioinformatics, 34(5):867–868, 2018. doi:10.1093/bioinformatics/btx699.

[22]

Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, and others. Scikit-learn: machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011. URL: https://jmlr.org/papers/v12/pedregosa11a.html.

[23]

Will Willison. Rich: rich text and beautiful formatting in the terminal. 2024. URL: https://github.com/Textualize/rich.