{"paper_id":"9c4e5e12-13d9-4f0d-ace8-6c46857bed44","body_text":"The Scalable Variant Call Representation:\nEnabling Genetic Analysis Beyond One Million Genomes\nTimothy Poterba1,2,3,†, Christopher Vittal1,2,3, Daniel King 1,2,3,4,†,\nDaniel Goldstein1,2,3, Jacqueline I. Goldstein 1,2,3, Patrick Schultz 1,2,3,\nKonrad J. Karczewski 1,2,3,4, Cotton Seed 1,2,3,4,*, Benjamin M. Neale 1,2,3,4,*\n1Program in Medical and Population Genetics,\nBroad Institute of MIT and Harvard, Cambridge, Massachusetts 02142, USA\n2Analytic and Translational Genetics Unit,\nMassachusetts General Hospital, Boston, Massachusetts 02114, USA\n3Stanley Center for Psychiatric Research,\nBroad Institute of MIT and Harvard, Cambridge, Massachusetts 02142, USA\n4Novo Nordisk Foundation Center for Genomic Mechanisms of Disease,\nBroad Institute of MIT and Harvard, Cambridge, Massachusetts 02142, USA\n*Indicates equal contribution\n†Corresponding authors: tpoterba@gmail.com, dking@broadinstitute.org\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nAbstract\nThe Variant Call Format (VCF) is widely used in genome sequencing but scales poorly . For\ninstance, we estimate a 150,000 genome VCF would occupy 900 TiB, making it both costly and\ncomplicated to produce and analyze. The issue stems from VCF’s requirement to densely\nrepresent both reference-genotypes and allele-indexed arrays. These requirements lead to\nunnecessary data duplication and, ultimately , very large ﬁles.\nTo address these challenges, we introduce the Scalable Variant Call Representation (SVCR). This\nrepresentation reduces ﬁle sizes by ensuring they scale linearly with samples. SVCR achieves\nthis by adopting reference blocks from the Genomic Variant Call Format (GVCF) and employing\nlocal allele indices. SVCR is also lossless and mergeable, allowing for N+1 and N+K incremental\njoint-calling.\nWe present two implementations of SVCR: SVCR-VCF, which encodes SVCR in VCF format,\nand VDS, which uses Hail’s native format. Our experiments con ﬁrm the linear scalability of\nSVCR-VCF and VDS, in contrast to the super-linear growth seen with standard VCF ﬁles. We\nalso discuss the VDS Combiner, a scalable, open-source tool for producing a VDS from GVCFs\nand unique features of VDS which enable rapid data analysis. SVCR, and VDS in particular,\nensure the scientiﬁc community can generate, analyze, and disseminate genetics datasets with\nmillions of samples.\nIntroduction\nThe pipeline for high-throughput sequencing involves a series of datatypes and the\ntransformations between them:\n1. Sequencing DNA to generate unaligned reads, often stored in a FASTQ.\n2. Aligning to a reference genome to generate aligned reads, often stored in a BAM/CRAM.\n3. Variant calling to generate genotype calls and metadata, often stored in a GVCF or VCF.\n4. Joint calling to generate an analysis-ready genotype matrix, often stored in a PVCF.\nSequencing, aligning, and variant calling are straightforwardly sample-parallel, but joint calling is\nnot. The latter necessarily combines multiple sample-oriented ﬁles into a single variant-oriented\nmatrix. In this paper, we focus exclusively on joint calling and its data formats.\nAs cohort sizes grew so did the challenge of combining many sample-major GVCFs into one\nvariant-major PVCF. These challenges have motivated the development of new tools (e.g.\nDRAGEN gVCF Genotyper, GLNexus), new representations (e.g. msVCF, spvcf), new\ncompressors (e.g. popvcf), and new formats (e.g. Savvy , Genomic Variant Store, GLNexus’s\nkey-value store, GenomicsDB). These e ﬀorts draw from three approaches: sparsity of variant\ngenotypes, sparsity of alleles, and amenability to sample parallelism. In this paper, we describe\na framework which employs all three techniques to address the above challenges and yield\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nexactly linear scaling of size in samples, sample parallelism in generation, incremental N+K\nmergeability , and a plaintext VCF representation.\nFormats and representations. The common representation of genomic sequences, Project VCF\n(PVCF), is untenable at scale. To see why , consider the addition of a single sample to a\nsequencing dataset with N samples. PVCF requires, for every novel variant in the new sample,\nan additional row containing one variant genotype and N homozygous reference genotypes.\nFurthermore, each new allele that the N+1th sample introduces requires extension of existing\nﬁelds to account for the new allele. For instance, the ﬁrst N samples’ allele depth (AD) ﬁelds need\na new “0” (no reads) entry for the new allele. Worse, the Phred-scaled genotype likelihood (PL)\nﬁeld grows quadratically in the number of alleles. The result is a repetitive and bloated\nrepresentation, where the same uninformative value is stored over and over again, leading to\nineﬃcient use of storage and computational resources.\nThroughout this paper we distinguish between a representation and a format. A representation, or\ndata model, deﬁnes the expected set of ﬁelds, what they mean, and how they are related. A\nformat describes a concrete implementation in terms of bytes. We similarly cleave the VCF\nSpeciﬁcation into the PVCF representation for cohort-level variant data and the VCF format for\nstoring a variant-by-sample matrix of arbitrary data types in tab-delimited plaintext.\nSpeciﬁcally , PVCF deﬁnes the semantics of ﬁelds such as GT, AD, GP , PL, and, for list ﬁelds, the\nrelationship between their length and the number of alternate alleles. VCF, as a format,\ndescribes, for example, how a number or a list is rendered in plaintext. The Variant Call Format\nsupports representations other than PVCF: the single-sample Genomic Variant Call Format\n(GVCF) and the structural variant VCF (SV VCF). Analogously , the PVCF representation\nsupports diﬀerent formats: VCF and Binary Call Format (BCF).\nV ariants scale with samples.When discussing the size of a dataset, we focus on the number of\ngenotype records. The number of genotype records in a GVCF does not vary substantially\nacross samples. We use K to refer to the average number of non-reference genotype records per\nsample, which, at 3-6 million, is much smaller than the number of bases in the human genome\n(3 billion).\nIn a multi-sample dataset, the expected number of loci at which at least one sample has a\nnon-reference genotype is a function of the number of samples, M(N). We call such loci variant\nsites, reserving the term variant to mean a particular allele at a particular locus. We use the term\ngenotype record to refer to a genotype call and its associated quality metrics. By de ﬁnition, M(1) =\nK and M(N), for large N, approaches the number of loci in the reference genome. In contrast, the\nnumber of variants is unbounded due to the in ﬁnitude of possible insertion and deletion\nvariants.\nConsider M(2). The second sample will share many variant sites with the ﬁrst sample; however,\nboth samples will have many unshared singletons. As a result, this second sample substantially\nincreases the number of variant sites in the dataset. In contrast, the millionth sample will add\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nsigniﬁcantly fewer variant sites because the previous 999,999 samples are likely to share many\nof its sites. Thus, M(N) is non-linear, which we con ﬁrm below.\nPVCF represents a collection of sequences as a dense matrix, with one column per sequenced\nsample and one row for every variant site. PVCF permits both a multiallelic representation\n(wherein each locus appears in at most one row) and a biallelic representation (wherein loci are\nrepeated such that each row contains one alternate allele). The number of genotype records in a\nmultiallelic PVCF grows as\nO(N × M(N))\nthe number of genotype records in a biallelic PVCF grows as\nO(N × M(N) × A(N))\nwhere A(N) is the expected number of alternate alleles at a variant site in a dataset with N\nsamples.\nUntil variants are observed at all loci, each sample brings with it some number of new singleton\nsites. We empirically measured the relationship between the number of variant sites and the\nnumber of samples in subsets of the PVCF representation of the HGDP+1kG dataset. We found\nthat the number of variant sites grows roughly with the square-root of N. Notice that we can\nmeasure the exponential relationship between M and N using a linear model in log-space.\nlog10(M) = a + b log 10(N)\nM = 10 a Nb\nTo explore this relationship, we generated multiallelic PVCFs from subsets of the HGDP+1kG\ndataset (Koenig et al. 2023). Figure 1a shows, for HGDP+1kG, M ≈ N 0.46. Figure 1b shows the\neﬀect of this growth on the number of genotypes. The same experiment applied to subsets of\ngnomAD v2 found M ≈ N 0.57. We stress that this is only evidence of an exponential relationship\nwhen N is between one and several thousand. Furthermore, this relationship depends on the\nancestry mix of the dataset: diverse ancestries should increase the rate at which new variant\nsites are discovered. As the number of samples grows, the growth of variant sites must\neventually decelerate, but not before sequencing datasets become impractically large.\nFields scale with alleles. In multiallelic PVCFs, not only does the number of genotype records\ngrow super-linearly , but the size of each genotype record increases with the size of the cohort. In\nparticular, sequencing datasets store per-allele quality metadata. For example, each genotype\nrecord stores an allele depth (AD) ﬁeld containing one integer per allele observed at this locus.\nThis ﬁeld grows linearly in the number of alleles at the locus even though most genotypes\nobserve at most two alleles. Even worse, the phred-scaled genotype likelihood (PL) ﬁeld stores the\nlikelihood of all possible genotypes given the alleles at this locus. Within the autosomes, which\nare diploid, this ﬁeld grows quadratically in the number of alleles. In practice, particularly in\nrepetitive regions, a small fraction of loci have thousands to tens of thousands of alleles.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nIn summary , ifA(N) is the average number of alleles per site, then the sum-total number of allele\ndepths in a multiallelic PVCF ﬁle grows as,\nO(N × M(N) × A(N))\nand the sum-total number of phred-scaled genotype likelihoods grows as,\nO(N × M(N) × A(N) 2)\nIn practice, at a small number of sites, we observe thousands of alternate alleles. Sites like these\nare particularly challenging to analyze because of the quadratic growth of the PL ﬁeld. We\nconjecture this quadratic growth explains why many PVCFs with more than a few thousand\nsamples have dropped the PL ﬁeld entirely or dropped it from multiallelic loci.\nNotice that both multiallelic and biallelic datasets su ﬀer from two forms of non-linear growth\n(genotypes and alleles). They di ﬀer only in where the growth appears. In the multiallelic\nrepresentation, both genotype records and allele indexed ﬁelds grow super-linearly in the\nnumber of samples. In the biallelic representation, allele-indexed ﬁelds are bounded to length\nthree but the number of genotype records grows even faster because highly multiallelic variant\nsites duplicate reference calls across every allele.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nFigure 1. The number of variant sites (a) and total non-reference genotypes (b) scale with the number of\nsamples. Each new sample brings new variant sites (a) which causes the dense matrix to grow\nsuper-linearly (b) in the number of samples. Dashed line indicates linear ﬁt in log-log space, using all\npoints with at least 100 individuals.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nImplementation\nWe propose a new representation for sequenced cohorts: the Scalable Variant Call\nRepresentation or SVCR. This representation is a generalization of GVCF to more than one\nsample. SVCR scales linearly in the number of samples. One or more GVCFs can be losslessly\nconverted to a single SVCR dataset. Two SVCR datasets can be losslessly merged into a new\nSVCR dataset comprising samples from both. An SVCR dataset can be stored in many formats\nincluding VCF, Hail native format, and Google BigQuery .\nOverview. A dataset in the SVCR representation consists of zero or more variant sites. Variant\nsites are indexed and identi ﬁed by genomic locus. There is at most one variant site per locus.\nVariant sites are multi-allelic. A dataset has a set of samples. Every variant site has a genotype\nrecord for every sample. A genotype record is either the unique missing value or a set of ﬁelds\nand their values. Any ﬁeld value may be missing. Any value permitted for a VCF FORMAT\nﬁeld is permitted as a ﬁeld value in SVCR.\nColumn-sparsity. Within a single sample, SVCR stores adjacent reference genotype records as\nintervals, which are termed “reference blocks”. Each reference block has a locus interval, which\nmust be contained within a single chromosome. A reference block is stored in the variant site\nidentiﬁed by the interval’s ﬁrst, or left-hand, locus. The size of the reference block, in base pairs,\nis stored in a genotype ﬁeld named LEN. All reference block intervals for any particular sample\nmust be disjoint.\nReference blocks are a form of column-sparsity by run-length encoding. At each locus, the\ngenotype record of samples with a homozygous reference genotype is encoded in the\nunderlying format as a missing value. The actual value is implicitly given by that sample’s\noverlapping reference block. Preservation of the reference blocks enables the combination of\ntwo or more SVCR datasets into a new SVCR dataset.\nLocal alleles. Consider a highly polymorphic low-complexity insertion-deletion locus with\n1,000 alternate alleles in a PVCF with 100,000 whole genomes. The VCF spec mandates that the\nallelic depth ﬁeld AD has one value per allele including the reference (R-numbered), so AD\nrecords at this site must contain 1,000 elements denoting the depth of reads supporting each\nallele. For any particular sample, it is likely that at most one or two of these alternate alleles are\nobserved. In the PVCF representation, the AD ﬁeld for such a sample would contain between\n998 and 999 zeros. The GVCF representation, in contrast, contains at most three entries. As\ndescribed previously , this problem is signiﬁcantly worse for ﬁelds like Phred-scaled genotype\nlikelihoods (PL), which have one value for each possible genotype con ﬁguration (G-numbered).\nData generators have mitigated this problem with approaches such as truncating the list of\nalternate alleles at any given locus to a maximum number, or disseminating PVCFs with only\ngenotype (GT) ﬁelds and no quality metadata. Each of these approaches has drawbacks:\ntruncating alternate alleles removes potentially pathogenic information from the dataset and\nmasks information from short tandem repeats (STRs), and dropping all quality ﬁelds makes\nquality control much more di ﬃcult.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nSVCR introduces a new entry-level ﬁeld, local alleles (LA). At each locus, for each sample, an LA\nﬁeld indicates which alleles were observed in this sample. The global alleles refers to the set of\nalleles observed in any sample at this locus. Concretely , the LA ﬁeld is an injective function from\nlocal allele indices to global allele indices (Figure 2). Since the number of local allele indices is\nﬁnite (in a diploid genotype, there are two), we represent the LA ﬁeld as an array . For example,\nat a locus with 1,000 global alleles, an LA value of [0, 4, 5] indicates that the ﬁrst local allele\ncorresponds to the zeroth (or reference) global allele, the second local allele corresponds to the\nfourth global allele, and the third local allele corresponds to the ﬁfth global allele. For\nconvenience, we require the ﬁrst element of LA to always be zero.\nFor each genotype-level allele-indexed ﬁeld from the VCF spec, we de ﬁne a new ﬁeld whose\nname begins with “L”. The ﬁeld only contains the elements corresponding to this sample’s local\nalleles. Consider a variant site with ﬁve alleles, A, AA, AAA, AAAA, AAAAA, and a sample\nwith the genotype AA/AAAA with 14 reads supporting AA and 16 reads supporting AAAA.\nThe local alleles array will be [0, 1, 3], the allelic depth will be [0, 14, 0, 16, 0], and the local allelic\ndepth array will be [0, 14, 16].\nAlthough the GT ﬁeld is not problematic, we strongly recommend using an LGT ﬁeld in SVCR.\nIn particular, when removing a sample, an implementation of SVCR may wish to remove global\nalleles which no longer appear in any samples. Removing a global allele requires modifying\nevery local alleles entry and also updating any globally-indexed ﬁelds, such as the GT.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nFigure 2. Local alleles. The AD ﬁeld is an array ﬁeld containing one element per allele, including the\nreference. The AD ﬁeld entries indicate, for each allele, the number of informative reads. The arrows\nterminating at an allele indicate to which allele each value corresponds. The LAD ﬁeld is a locally indexed\narray ﬁeld storing the same data. Note that we need not store information for the unobserved C allele, and\nthat the two-step paths lead to the same alleles as the one-step path.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nSVCR-VCF .We encode SVCR in VCF by requiring two FORMAT ﬁelds: LA, an array of\nintegers, and LEN, an integer (described above). Reference blocks must have a non-missing\nLEN, should have a non-missing DP & GQ, may have a non-missing LGT, and should have a\nmissing LA, LPL, and LAD. When a sample has an overlapping reference block at a locus, all of\nthat sample’s FORMAT ﬁeld values must be missing. Non-reference genotypes must have a\nmissing LEN, must have a non-missing LA, and should have a non-missing DP, GQ, LGT, LAD,\nand LPL.\nSVCR introduces three new settings for the “Number” of a VCF list ﬁeld:\n● LOCAL-A numbered ﬁelds should have one list element for each allele present in LA\nminus one (to exclude the reference). For example: LEC.\n● LOCAL-R numbered ﬁelds should have one list element for each allele present in LA\n(including the reference). For example: LAD.\n● LOCAL-G numbered ﬁelds should have one list element for each possible genotype\ndepending on the number of alleles N LA and the ploidy of the genotype call. If diploid,\n(NLA * (NLA+1)) / 2. If haploid, N LA. For example: LPL.\nFigure 3 demonstrates the equivalent representation of twenty base pairs of two samples in\nGVCF, PVCF, and SVCR-VCF.\nVDS. We also encode SVCR in a Hail native format called Variant Dataset or VDS (unrelated to\nthe Hail 0.1 VariantDataset except in name). A VDS comprises two Hail matrix tables. One Hail\nmatrix table contains reference data the other contains variant data. The reference matrix table\nhas rows keyed by locus, columns keyed by sample identi ﬁer, and exactly three entry ﬁelds:\nLEN (int32), DP (int32), and GQ (int32). Entries which are not the start of a reference block are\nmissing values.\nThe variant matrix table has rows keyed by locus and alleles, columns keyed by sample id, and\nat least two entry ﬁelds, LA (array<int32>) and LGT (call), and usually also has LAD\n(array<int32>), LPL (array<int32>), DP (int32), GQ (int32), and gvcf_info (struct<...>). The latter\ncontains all the GVCF INFO ﬁelds for that sample.\nFor data interchange, we recommend a format that uses a single ﬁle. For analysis, we\nrecommend storing reference and variant data in two ﬁles for three reasons:\n1. The schemas for reference and variant records di ﬀer; therefore, a combined\nrepresentation pays the overhead of mixing these schemas. We observed a roughly 10%\nsize reduction in both VCF and VDS from a split format.\n2. The split representation makes it easier to ﬁlter variants without inadvertently removing\nreference blocks.\n3. Analytical and quality control methods that operate on the smaller variant-only table are\nfaster and cheaper because they read substantially less data. For example, the\nnewly-developed CHARR (Lu et al., 2023) achieves substantial speedups as it only needs\nto access variant data (homozygous genotypes).\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nSparsity in representation versus sparsity in format. It should be noted that while SVCR is\nsparse—only a small number of matrix entries contain information—both the GZIP-compressed\nVCF and Hail VDS realizations of SVCR contain high numbers of explicit “missing” sentinel\nvalues. Both formats use compression techniques which rapidly & e ﬃciently compress and\ndecompress these repeated values, reducing the gain from switching to an explicitly sparse\nformat. We anticipate substantial reductions in analysis cost and time from a system that stores\nand computes directly on the sparse representation.\nDensiﬁcation. The reference information for all samples at some locus L is not contained within\na single row of the matrix at L. Instead, this information is stored in overlapping reference\nblocks that appear at their start location, which is necessarily some previous locus. We call\nproducing a single row containing all reference and variant information for all samples at a\nparticular locus densiﬁcation. In general, densifying a locus requires searching all previous loci\nuntil an overlapping reference block for every sample has been identi ﬁed.\nMitigation by periodic reference checkpoints. We propose two strategies to bound the size of\neach reference block and thus bound the search distance. Reference blocks could either be split\nat some ﬁxed size limit or at ﬁxed intervals over the genome. Both strategies add reference\nblocks and therefore trade increased space for reduced densi ﬁcation time. Both strategies have a\nparameter K such that, at a particular variant site, for all samples, all overlapping reference\nblocks begin no further than K base pairs prior.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\n(3a)\n(3b)\n(3c)\nFigure 3. Representing sequences. In each panel, two di ﬀerent samples, 01 and 02, are shown in di ﬀerent\nrepresentations, including a) two GVCF ﬁles, b) a PVCF ﬁle, and c) an SVCR-VCF ﬁle.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nSplitting at a ﬁxed period. Split reference blocks at every K bases. Store reference blocks for all\nnon-variant samples at each split location. If this dataset is also partitioned into K-sized\npartitions, the information necessary to realize any partition’s dense sub-matrix of genotypes is\nfully contained in that partition. This approach in ﬂates the reference data by splitting every\nreference block that spans the K-th locus even if that reference block is already tiny . As shown\nbelow, this inﬂation relative to the next strategy is very modest.\nSplitting at a size threshold. Split reference blocks spanning more than K bases into two or more\nreference blocks spanning at most K bases. For instance, with a value of K=10,000 a reference\nblock starting at chr1:4500 with length 38000 would be split into reference blocks\nchr1:4500-14499, chr1:14500-24499, chr1:24500-34499, chr1:34500-42499. This strategy introduces\nthe minimum number of new reference blocks to ensure the K invariant; however, in partitioned\ndatasets, data from previous partitions may be necessary to realize the dense sub-matrix of\ngenotypes at any partition after the ﬁrst.\nK is a parameter of the dataset. We do not mandate or recommend any particular value of K\nbecause the length distribution of reference blocks substantially depends on experimental\ndesign, variant calling, and reference block bin resolution. For example, whole exome\nsequencing produces shorter reference blocks due to the ﬂuctuations in depth from partial\ncapture of the genome. In contrast, high coverage whole genome sequencing can produce long\nblocks of similar-quality reference calls in regions of constant depth. Table 1 and Table 2 present\nthe measured number of reference blocks after splitting the HGDP+1kG dataset using each\nstrategy .\nIn our experience, the granularity of reference block GQ bins has a large impact on the average\nlength of reference blocks in a given dataset. Whole genome sequencing data also tends to have\nlonger reference blocks than whole exome data, due to the absence of ﬁxed capture boundaries\nand more consistent coverage across the genome.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nK (fixed period) Billions of Reference\nBlocks after Splitting\nRatio of Split to Unsplit\nReference Blocks\n50,000 73.1 1.0000\n10,000 73.1 1.0002\n5,000 73.2 1.0015\n1,000 77.7 1.0634\n500 87.1 1.1925\n100 175.5 2.4013\nT able 1. Data size is inversely proportional to the period when splitting at a ﬁxed period. Using this\nstrategy, a period of 5,000 bases introduces a modest 0.15% increase in size. The increased sizes of other\nperiods are also shown. At a period of 50,000 bases, only a few hundred reference blocks are added.\nK (maximum size)\nBillions of Reference\nBlocks after Splitting\nRatio of Split to Unsplit\nReference Blocks\n50,000 73.1 1.0000\n10,000 73.1 1.0002\n5,000 73.2 1.0015\n1,000 77.7 1.0634\n500 87.1 1.1922\n100 175.3 2.3985\nT able 2. Data size is inversely proportional to the threshold when splitting at a size threshold. In practice,\ncompared to splitting at a ﬁxed period, this strategy is only mildly more parsimonious.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nStreaming densiﬁcation. The previous discussion and mitigations address reading one\nparticular locus. We call this a “point-query”. In contrast, applications that read all or most of\nthe dataset do not need these mitigations. Instead, streaming densi ﬁcation keeps one\noverlapping reference block per sample, swapping in new reference blocks as they appear in the\nstream of records. As stated above, this dense matrix is super-linear in size. Although a\nstreaming application does not reify the entire matrix in memory , if it performs computation on\neach dense entry , the compute-time, and thus cost, will grow super-linearly . For example,\ncomputing the allele frequency necessarily inspects averages over every call, even the\nhomozygous reference calls. A clever implementation would avoid densi ﬁcation entirely .\nGVCF is a single-sample SVCR. SVCR is a multi-sample generalization of GVCF. Indeed, an\nSVCR with one sample is equivalent to a GVCF. In a GVCF, the global alleles and the local\nalleles are the same because there is only one sample. In this setting locally-indexed and\nglobally-indexed ﬁelds are equivalent. Speci ﬁcally , in a GVCF, theAD and PL ﬁelds must be\nequal to the LAD and LPL ﬁelds. The lossless transformation between a GVCF and SVCR-VCF is\nas simple as (1) adding a ﬁeld LA to create an identity function of the alleles, (2) adding an “L”\npreﬁx to the name of allele-indexed ﬁelds, (3) adding a LEN ﬁeld to the FORMAT derived from\nthe position and the END ﬁeld.\nHierarchical joint calling. An SVCR may be produced in a tree-structured manner because it is\ncombinable with itself and with GVCFs. This allows some sample-wise parallelism in the\nproduction of an SVCR dataset from GVCFs. Suppose we have 100,000 GVCFs. We may , in\nparallel, combine 1,000 groups of 100 GVCFs each into 1,000 SVCR datasets. We may then, in\nparallel, combine ten groups of 100 SVCR datasets into ten SVCR datasets. Finally , we may\ncombine the ten SVCR datasets into one SVCR dataset. This entire process is also fully variant\nparallel.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nFigure 4. The graph of interconversion and combination operations. Notice that PVCF is a sink: crucial\ninformation is lost when producing a PVCF—missing data & homozygous reference calls can not be\nunambiguously resolved if the quality metrics at non-variant sites are discarded.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nResults\nJoint calling & analysis of nearly one million exomes. SVCR has been used since 2019 to\nmerge, store, and analyze large sequencing cohorts. The Hail library includes a “VDS\nCombiner” that can combine zero or more VDSes with zero or more GVCFs into a single VDS in\na parallel and fault-tolerant manner. Five large datasets produced by the VDS Combiner directly\nfrom GVCF ﬁles are:\n1. gnomAD v3.1. 153,000 genomes. 2019. (Chen, et al. 2020)\n2. Center for Common Disease Genomics (CCDG) exomes. 203,000 exomes. 2021.\n(Felsenfeld 2018)\n3. CCDG genomes. 136,000 genomes. 2021. (Felsenfeld 2018)\n4. gnomAD v4. 955,000 exomes. 2021 (Chao, et al. 2023).\n5. Blended Genome Exome (BGE) Wave 1. 82,000 exomes. 2023. (Howrigan 2022)\nThe All of Us research project also used the VDS Combiner, but started from tens of Avro ﬁles\neach containing 4,000 complete samples in a GVCF-like representation.\n6. All of Us, April 2023 data freeze. 245,350 genomes. (All of Us, 2022)\nIn addition to producing the VDS for each of these datasets, the Hail system was used to\nperform quality control and analyze each dataset.\nThe VDS Combiner is implemented in Hail Query , an open-source, partitioned,\nhorizontally-scalable, spot-tolerant, dataframe system and genomic analysis library with a\nPython API. Hail Query , and therefore the VDS Combiner, can run on any Apache Spark or Hail\nBatch cluster. A managed version of the former can be found in all the major clouds including\nAWS, Google, Azure, and Alibaba.\nHail Query is pervasively designed to support spot instances. At the time of writing, the cost of\na spot core as a fraction of the non-spot price of general purpose families in AWS was about\n40-50%; in GCP , about 21-40%; and in Azure, about 10-25%.\nUnlike other joint calling systems, Hail Query is also a general purpose analysis system\nsupporting relational algebra (e.g. ﬁlter, aggregate, group-by , order-by), distributed linear\nalgebra (e.g. PCA, matrix multiplication), and export to many formats (e.g. BGEN, PLINK, VCF,\nTSV). The Hail Query language, which is used to manipulate values within datasets, supports a\nwide variety of data types (e.g. integers, ﬂoating-point numbers, strings, genomic loci, genotype\ncalls, arrays, sets, dictionaries, tuples, ndarrays, and structs) and a wide variety of operations on\nthese types (e.g. random functions, statistical distributions, statistical tests, linear and logistic\nregression, linear algebra, collection iteration & aggregation).\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nFigure 5. Scaling of dataset size of sequenced whole genomes represented in Hail LZ4-compressed VDS,\ngzip-compressed SVCR-VCF and gzip-compressed PVCF . In (a) total size in Gibibytes is plotted on the y\naxis, and in (b), plotting the size per sample (total size divided by number of samples) provides a clear\nview of the super-linear scaling of PVCF .\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nSize in bytes. Consider subsets of the HGDP+1kG dataset in PVCF, SVCR-VCF, and VDS. The\ntotal size grows super-linearly with sample size for PVCF but linearly for SVCR-VCF and VDS\n(Figure 5a). Further con ﬁrming this, we see the size per sample increases for PVCF, while\napproaching an asymptote for SVCR-VCF and VDS (Figure 5b), as expected since each human\ngenome contains an approximately ﬁxed amount of information.\nIt is critical to note that exact ﬁgures of size per sample should not be used to compare\nrepresentations and formats without additional context that re ﬂects the granularity and type of\ninformation contained. For example, choices in resolution of reference blocks and the selection\nof included FORMAT ﬁelds can lead to a substantial size di ﬀerence in single-sample GVCFs for\nthe same sequencing experiments, and these ratios are preserved in SVCR.\nWe did not systematically generate PVCFs or SVCR-VCFs for larger sample sizes due to the\nhigh cost of these experiments, but we show that this asymptotic trend persists in a VDS at\nlarger cohort sizes with subsets of gnomAD (Figure 6), and list a few selected examples of PVCF\nand VDS size in Table 3.\nCost of joint calling. As of August 2023, using Hail 0.2.120, a VDS produced from GVCFs costs\napproximately 0.005 USD per exome, and in a large-scale application of the VDS combiner for\nwhole genomes in 2019, we observed a cost of 0.10 USD per genome. These numbers represent\nthe sum total cost of using a Google Dataproc cluster with n1-standard spot instances to\ngenerate a VDS.\nCost of storage. The VDS is stored as a hierarchy of ﬁles or blobs. As mentioned before, the size\nof a VDS depends signi ﬁcantly on the granularity of reference blocks. As a concrete example, at\n260 MiB per genome or 26 MiB per exome, a VDS costs about 0.0052 USD per genome per\nmonth or 0.00052 USD per exome per month to store in a cloud object store (at a typical cost of\n0.02 USD per GiB per month).\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nFigure 6. The growth of size in bytes per sample of subsets (a) gnomAD v3.1 exomes and (b) of gnomAD\nv4 exomes represented as VDS.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nDescription Data type Year Samples Size (TiB) Size per sample (MiB)\nPVCF GVCFs VDS PVCF GVCFs VDS\ngnomAD v2 genomes 2016 15,496 16.00 5.00 1,082.68 338.34\ngnomAD v3.1 genomes 2019 153,030 37.66 258.05\nCCDG exomes exomes 2021 203,664 3.36 17.30\nCCDG genomes genomes 2021 136,959 108.73 832.45\ngnomAD v4 exomes 2021 955,359 18.30 20.09\nBGE Wave 1 exomes 2023 82,000 8.76 2.07 112.02 26.47\nAll of Us April 2023 Data Freeze genomes 2023 245,350 292.61 20.95 1,250.56 89.54\nT able 3. Sizes and size per sample of selected large callsets. A blank PVCF cell indicates that a PVCF was\nnever generated. A blank GVCF cell indicates the sum-total size of the GVCFs is no longer available. The\nlarge variability in size per sample is due to variability in the granularity of reference blocks. In\nparticular, integer precision GQs produce very large ﬁles. GVCF sizes should not be compared directly to\nother representations to assess compression, because in many cases some ﬁelds are dropped during or after\ncombining GVCFs.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nConclusion\nA concise, portable, & analyzable representation. We have presented the Scalable Variant Call\nRepresentation, a concise representation of genome sequencing datasets that enables quality\ncontrol and analysis of nearly one million whole exome sequences. We rei ﬁed SVCR in the VCF\nformat as SVCR-VCF. We also presented the Hail VDS, a format generated and consumed by the\nopen-source Hail Query library .\nJoint calling. “Joint calling” has historically referred to combining samples into one dataset and\nalso possibly adjusting genotypes given information from other samples. In this paper, we use\nthe term “joint calling” to mirror conventional descriptions of these processes; however,\n“combining” may be more precise. The Hail VDS Combiner does not adjust genotypes nor is the\nSVCR format aﬀected by adjustment, although notably , GATK-based “joint calling” pipelines\nsuch as GenotypeGVCFs do not do so either. Instead, SVCR and the VDS Combiner provide a\nrepresentation and platform on which genotype adjustment could be implemented, but the\nprimary challenge at scale is the transposition of single-sample GVCF ﬁles into a variant-major\nlayout.\nStandardization. Local allele indexing has been under discussion for a few years at the\nhts-specs repository under issue 434 (Farjoun 2019). We hope the release of this paper motivates\nﬁnalizing those changes. Further work remains to standardize reference blocks.\nFuture directions. Sparse representations of genotype matrices o ﬀer an eﬃcient way to store\nlarge-scale sequenced cohorts, but many analyses still expect dense matrices as inputs.\nDeveloping new methods for querying and modeling genotype matrices that operate directly on\nsparse and/or compressed data will lead to substantial savings of compute time in addition to\nstorage. Further, representing missing genotypes in VDS or a split SVCR format using a third\ntable (in addition to reference and variant data) can support e ﬃcient analysis when reference\nquality information is not necessary by obviating the need to read the reference data; however,\ndetermining the best approach to do so in the context of dynamic ﬁltering makes this\nchallenging.\nRelated work.\nGLNexus is both a tool and a particular modi ﬁcation of the PVCF format. (Lin, et al. 2018).\nGLNexus combines GVCFs into PVCFs in three steps: (1) ﬁnd a uniﬁed set of loci, (2) “project”\ngenotype-level ﬁelds from each GVCF into the PVCF, and (3) adjust genotype calls.\nGLNexus proposes a new allele uni ﬁcation algorithm which introduces a new PVCF variant\ntype: <MONOALLELIC>. This algorithm and new variant are designed to generate variants\nthat are “completely non-overlapping, with mutually-exclusive alleles”.\nGLNexus stores genotype records in a key-value store. GLNexus splits the reference genome\ninto a disjoint set of 30 kilobase “bins”. Each genotype record is associated with a key pair of the\nbin containing this record and the sample identi ﬁer. GLNexus uses an LSM-tree based key-value\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nstore which permits the addition of new samples in amortized linear time and space. Genomic\nintervals may also be retrieved in amortized linear time.\n(Lin, et al. 2019) and (Yun, et al. 2021) report the runtime and PVCF size for several applications\nof GLNexus to produce a PVCF from GVCFs, including an exome cohort of nearly 250,000\nsamples and a genome cohort of nearly 23,000 samples. We summarize their results in Table 4.\nThe innovations of GLNexus and SVCR address di ﬀerent PVCF challenges. SVCR, as de ﬁned\nabove, requires one row per locus; however we believe GLNexus PVCFs are otherwise\ncompatible with local allele indexing and preservation of reference blocks. The GLNexus\ndatabase preserves reference information and e ﬀectively uses local allele indexing because it\nexactly preserves the GVCF records. The 30 kilobase reference bins resolve the challenge of\npoint-queries for a dense vector in a manner similar to our “split at ﬁxed periods” approach.\nThe Hail VDS Combiner does not adjust genotype calls; however we believe the GLNexus\nempirical frequency prior could be implemented in Hail Query .\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nDescription Reference Samples Time\n(core-hours / sample)\nSize\n(MB / sample)\nChromosome 2 of an\nexome,\nunspeciﬁed version\nLin, et. al 2019 50,000 0.005 Not reported.\nExome,\nstandalone GLNexus\nLin, et. al 2019 16,521 0.058 Not reported.\nExome,\ndistributed GLNexus\nLin, et. al 2019 243,953 0.23 (“thread”-hours) 28\nGenome,\ndistributed GLNexus\nLin, et. al 2019 22,609 7 Not reported.\nGenome,\nunspeciﬁed version\nYun, et. al 2021 2,504 Not reported. 387\nT able 4. GLNexus performance. We collect here the time to generate and size of PVCFs generated by\nGLNexus.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nGATK’s GenotypeGVCFs command produces a PVCF from either one or more GVCFs or a\nGenomicsDB workspace. The largest callset known to the authors to use GenotypeGVCFs\ndirectly on GVCFs is the 91,796 sample ExAC un ﬁltered callset (Lek et al. 2016). The largest\ncallset to use GenomicsDB is the 164,332 individuals with exome sequence data and 20,314\nindividuals with genome sequence data from the un ﬁltered gnomAD v2 callset (Karczewski et\nal. 2020). GenotypeGVCFs su ﬀers from both of the problems described above: for every variant\nsite, every sample must have a genotype and for every allele, every genotype must have quality\nmetrics.\nGenomicsDB (Datta 2017) and TileDB (Papadapoulos 2016) are projects with a shared origin\nwhich use similar techniques of sparse matrix representations to support better scaling. TileDB\nis an open-source multidimensional array storage manager written in C++ with support for\nlarge dense and sparse arrays. TileDB-VCF is an open-source C++ library with Python and\ncommand line interfaces which includes schemas and functionality for variant stores, including\nthe ability to ingest new samples and query slices. TileDB-VCF stores GVCF variant and\nreference block records directly without PL expansion in a sparse array , indicating that\nTileDB-VCF disk footprints will scale linearly in the number of samples. GenomicsDB is a\nseparately-maintained open-source project originally developed on TileDB with the same data\nmodel described above, and tighter integration with GATK.\nThe multi-sample VCF (msVCF) and DRAGEN gVCF genotyper were developed by Illumina for\ngenerating and representing a jointly called dataset (Illumina 2023) (Schulz-Triegla ﬀ 2023). The\nDRAGEN gVCF genotyper accepts either GVCFs or “multi-sample GVCFs” as input and\nproduces an msVCF as output. This process involves three steps:\n1. DRAGEN aggregates batches of GVCFs into census ﬁles and cohort ﬁles. A census ﬁle\nstores variant metadata and per-sample reference blocks for a batch. A cohort ﬁle is a\nGVCF-like format containing all the samples of a batch.\n2. DRAGEN aggregates all census ﬁles into a single global census ﬁle.\n3. DRAGEN uses the global census ﬁle to generate, per cohort ﬁle, an msVCF.\nThe ﬁrst step is parallelizable per batch of GVCFs. If desired, the second and third steps can be\ndone in a tree-like fashion which allows both scaling of samples as well as incremental addition\nof samples. The DRAGEN gVCF genotyper is available both on-sequencer and in the cloud.\nThe msVCF uses a slight variant of local allele indexing. Instead of “LA” they use “LAA” and\ninstead of requiring and including an entry for the reference allele, they always elide the\nreference allele. In contrast, this paper requires the ﬁrst entry of LA to always be 0. The msVCF\ndeﬁnes the LPL and LGT ﬁelds equivalently to this paper. Similarly to VDS, but unlike\nSVCR-VCF, DRAGEN stores the reference blocks in a separate ﬁle, the census ﬁle. We believe a\nsplit version of SVCF-VCF could serve the needs of DRAGEN while also being interoperable\nwith other tools.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nDRAGEN 4.0 has generated a 500,000 sample cohort (Schulz-Triegla ﬀ 2023). Illumina’s\ncloud-based DRAGEN platform, Illumina Connected Analytics, has been tested with cohorts of\nat least 100,000 samples. Illumina reports a cost of “0.3 iCredits per sample”. The cost of an\niCredit depends on the volume of iCredits purchased, but is roughly 1 USD as of 2023-11-20.\nDRAGEN also supports a “compact GVCF” representation which omits statistics only necessary\nfor de novo variant calling in pedigrees. This innovation is compatible with the SVCR and VDS\nCombiner.\nThe Genomic Variant Store, developed by the Broad Institute’s Data Sciences Platform (DSP), is a\nGoogle BigQuery based solution to GVCF aggregation and variant storage. The GVS\nrepresentation is a variant of SVCR. It consists of three tables: vet (variants), ref_ranges (reference\nblocks), and alt_allele (variants). The vet and ref_ranges table are indexed by genomic position.\nThe alt_allele table is indexed by sample. Much like GLNexus, it preserves GVCF records and\nthus eﬀectively uses local allele indexing. GVS stores reference and variant data in distinct\ntables; indeed, it inspired the Hail VDS split representation. Storing the variant information\ntwice, once row-major and once column-major appears unique to GVS and allows for rapid\naccess to a single sample’s sequence. GVS can export to both PVCF and Hail VDS. As of October\n2022, GVS can produce a 10,000 sample joint callset in less than half a day at a cost of 0.06 USD\nper genome (Degatano 2022). GVS, using Hail Query and BigQuery together, produced the All\nof Us April 2023 freeze callset, a 245,000 whole genome VDS.\nOther projects. Multiple projects propose better encodings and compression of PVCF data. Savvy\n(LeFaive 2021) is a storage layer and C++ query API for e ﬃcient storage and queries of variant\ndata. Savvy’s data model has the same scaling characteristics as PVCF. Lossless textual\nencodings of PVCF ﬁles have also been proposed. Lin et al. proposed spVCF which achieves\nO(N1.1) scaling of bytes in samples (Lin 2020). Eggertsson et al. propose popvcf which\ndemonstrates further improvements over spVCF but they do not report a scaling factor\n(Eggertsson 2022). Unlike Savvy , spVCF, and popvcf, the size in bytes of SVCR and VDS both\nscale as O(N) in the samples. Moreover, the VDS is integrated with the scalable Hail Query\ndataframe and linear algebra system.\nBCFTools supports local allele indexing, albeit using the LAA name and omitting the entry\nmapping the local reference allele to the global reference allele (Li 2023). This de ﬁnition of LAA\nmatches that of msVCF, described above.\nAcknowledgements\nWe thank Mike Wilson for help with HGDP data, as well as Katherine Chao, Grace Tiao, and the\ngnomAD consortium for making data available for methods development. We thank the rest of\nthe Hail team and its alumni for their development of the Hail system. We thank the Data\nSciences Platform (formerly Data Science and Data Engineering) broadly and especially Laura\nGauthier, Yossi Farjoun, Eric Banks, and Louis Bergelson for discussion and collaboration\naround early prototypes. We also thank the Variants team at the Data Sciences Platform at the\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nBroad for a productive and interesting collaboration as well as feedback on this document. We\nthank Lee Lichtenstein for data on the AoU VDS. This work is supported by 1U01MH115727-01,\n1UM1HG008900-01, U24HG011450-03, 1U01MH115727-01, 5U54DK105566-03,\n2R01MH107649-04, 5U01HG009088-03, U24HG011450, the Chan Zuckerberg Initiative, as well\nas The Stanley Center for Psychiatric Research, which together with the Neale Lab has provided\nan incredibly supportive and stimulating home. This work is also supported by the Novo\nNordisk Foundation (NNF21SA0072102) and R37MH107649. We thank the Jeremy M. and Joyce\nE. Wertheimer Foundation, whose strategic advice and generous philanthropy have been\nessential for growing the impact of Hail. The content is solely the responsibility of the authors\nand does not necessarily represent the o ﬃcial views of any funding source.\nReferences\nAll of Us Research Program. Genomic Research Data Quality Report\nAll of Us Curated Data Repository (CDR) release C2022Q4R9.\nhttps://support.researchallofus.org/hc/en-us/article_attachments/17973653017236/_QC_Rep\nort_v7_release.pdf\nFelsenfeld A. Centers for Common Disease Genomics. 2018-09-26. Retrieved 2023-11-20.\nhttps://www.genome.gov/Funded-Programs-Projects/NHGRI-Genome-Sequencing-Program\n/Centers-for-Common-Disease-Genomics.\nhttps://web.archive.org/web/20230710204657/https://www.genome.gov/Funded-Programs-\nProjects/NHGRI-Genome-Sequencing-Program/Centers-for-Common-Disease-Genomics.\nChao, Katherine, & gnomAD Production Team. “gnomAD v4.0”. 2023-11-01. Retrieved\n2023-11-20. https://gnomad.broadinstitute.org/news/2023-11-gnomad-v4-0/.\nhttps://web.archive.org/web/20231103034332/https://gnomad.broadinstitute.org/news/202\n3-11-gnomad-v4-0/.\nChen, S.*, Francioli, L. C.*, Goodrich, J. K., Collins, R. L., Wang, Q., Alföldi, J., Watts, N. A.,\nVittal, C., Gauthier, L. D., Poterba, T., Wilson, M. W., Tarasova, Y ., Phu, W., Yohannes, M. T.,\nKoenig, Z., Farjoun, Y ., Banks, E., Donnelly , S., Gabriel, S., Gupta, N., Ferriera, S., Tolonen, C.,\nNovod, S., Bergelson, L., Roazen, D., Ruano-Rubio, V ., Covarrubias, M., Llanwarne, C., Petrillo,\nN., Wade, G., Jeandet, T., Munshi, R., Tibbetts, K., gnomAD Project Consortium,\nO’Donnell-Luria, A., Solomonson, M., Seed, C., Martin, A. R., Talkowski, M. E., Rehm, H. L.,\nDaly , M. J., Tiao, G., Neale, B. M.†, MacArthur, D. G.† & Karczewski, K. J. A genomic mutational\nconstraint map using variation in 76,156 human genomes. Nature. 2024 Jan;625(7993):92-100.\nPMID: 38057664. doi: 10.1038/s41586-023-06045-0\nDanecek P , Bonﬁeld JK, Liddle J, Marshall J, Ohan V , Pollard MO, Whitwham A, Keane T,\nMcCarthy SA, Davies RM, Li H. Twelve years of SAMtools and BCFtools. Gigascience. 2021 Feb\n16;10(2):giab008. doi: 10.1093/gigascience/giab008. PMID: 33590861; PMCID: PMC7931819.\nhttps://doi.org/10.1093/gigascience/giab008\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nDegatano, Kylee. 2022. “Scaling variant discovery to a million genomes with the Genomic\nVariant Store”.\nhttps://terra.bio/scaling-variant-discovery-to-a-million-genomes-with-the-genomic-variant-sto\nre/.\nhttps://web.archive.org/web/20230324165159/https://terra.bio/scaling-variant-discovery-to-\na-million-genomes-with-the-genomic-variant-store/.\nEggertsson, Hannes. popvcf. https://github.com/DecodeGenetics/popvcf.\nhttps://web.archive.org/web/20220203165936/https://github.com/DecodeGenetics/popvcf\nFarjoun, Yosi. 2019. “Deﬁne Local Alleles in VCF to allow for sparser format #434”.\nhttps://github.com/samtools/hts-specs/pull/434.\nhttps://web.archive.org/web/20221230235434/https://github.com/samtools/hts-specs/pull/\n434.\nHowrigan, D., DeFelice, M., Grimsby , J., Blumenstiel, B., Holmes, L., Ferriera, S., ... & Buxbaum,\nJ. (2022). REPLACING GWAS ARRA YS: CAPTURING GENOMIC DIVERSITY WITH A NOVEL\nWHOLE-EXOME PLUS LOW-PASS WHOLE GENOME PRODUCT. European\nNeuropsychopharmacology , 63, e25-e26.\nIllumina. “Iterative gVCF Genotyper Analysis”. 2023. Retrieved 2023-11-20.\nhttps://support-docs.illumina.com/SW/dragen_v42/Content/SW/DRAGEN/gVCFGenotype\nr.htm.\nhttps://web.archive.org/web/20231120223510/https://support-docs.illumina.com/SW/drage\nn_v42/Content/SW/DRAGEN/gVCFGenotyper.htm\nKarczewski KJ, Francioli LC, Tiao G, Cummings BB, Alföldi J, Wang Q, Collins RL, Laricchia\nKM, Ganna A, Birnbaum DP , Gauthier LD, Brand H, Solomonson M, Watts NA, Rhodes D,\nSinger-Berk M, England EM, Seaby EG, Kosmicki JA, Walters RK, Tashman K, Farjoun Y , Banks\nE, Poterba T, Wang A, Seed C, Whi ﬃn N, Chong JX, Samocha KE, Pierce-Ho ﬀman E, Zappala Z,\nO'Donnell-Luria AH, Minikel EV , Weisburd B, Lek M, Ware JS, Vittal C, Armean IM, Bergelson\nL, Cibulskis K, Connolly KM, Covarrubias M, Donnelly S, Ferriera S, Gabriel S, Gentry J, Gupta\nN, Jeandet T, Kaplan D, Llanwarne C, Munshi R, Novod S, Petrillo N, Roazen D, Ruano-Rubio\nV , Saltzman A, Schleicher M, Soto J, Tibbetts K, Tolonen C, Wade G, Talkowski ME; Genome\nAggregation Database Consortium; Neale BM, Daly MJ, MacArthur DG. The mutational\nconstraint spectrum quantiﬁed from variation in 141,456 humans. Nature. 2020\nMay;581(7809):434-443. doi: 10.1038/s41586-020-2308-7. Epub 2020 May 27. Erratum in: Nature.\n2021 Feb;590(7846):E53. Erratum in: Nature. 2021 Sep;597(7874):E3-E4. PMID: 32461654; PMCID:\nPMC7334197. https://www.nature.com/articles/s41586-020-2308-7\nKoenig Z, Yohannes MT, Nkambule LL, Goodrich JK, Kim HA, Zhao X, Wilson MW, Tiao G,\nHao SP , Sahakian N, Chao KR; gnomAD Project Consortium; Rehm HL, Neale BM, Talkowski\nME, Daly MJ, Brand H, Karczewski KJ, Atkinson EG, Martin AR. A harmonized public resource\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nof deeply sequenced diverse human genomes. bioRxiv [Preprint]. 2023 Aug\n10:2023.01.23.525248. doi: 10.1101/2023.01.23.525248. PMID: 36747613; PMCID: PMC9900804.\nLeFaive J, Smith A V , Kang HM, Abecasis G. Sparse allele vectors and the savvy software suite.\nBioinformatics. 2021 Nov 18;37(22):4248-4250. doi: 10.1093/bioinformatics/btab378. PMID:\n33989384; PMCID: PMC9502232. https://doi.org/10.1093/bioinformatics/btab378\nLi H, Handsaker B, & Danecek P . “bcftools(1) Manual Page”. 2023-05-30. Retrieved 2023-11-20.\nhttps://samtools.github.io/bcftools/bcftools.html.\nhttps://web.archive.org/web/20231030232455/https://samtools.github.io/bcftools/bcftools.h\ntml.\nLin MF, Bai X, Salerno WJ, Reid JG. Sparse Project VCF: e ﬃcient encoding of population\ngenotype matrices. Bioinformatics. 2021 Apr 1;36(22-23):5537-5538. doi:\n10.1093/bioinformatics/btaa1004. PMID: 33300997; PMCID: PMC8016461.\nhttps://doi.org/10.1093/bioinformatics/btaa1004\nLek M, Karczewski KJ, Minikel EV , Samocha KE, Banks E, Fennell T, O'Donnell-Luria AH, Ware\nJS, Hill AJ, Cummings BB, Tukiainen T, Birnbaum DP , Kosmicki JA, Duncan LE, Estrada K, Zhao\nF, Zou J, Pierce-Hoﬀman E, Berghout J, Cooper DN, De ﬂaux N, DePristo M, Do R, Flannick J,\nFromer M, Gauthier L, Goldstein J, Gupta N, Howrigan D, Kiezun A, Kurki MI, Moonshine AL,\nNatarajan P , Orozco L, Peloso GM, Poplin R, Rivas MA, Ruano-Rubio V , Rose SA, Ruderfer DM,\nShakir K, Stenson PD, Stevens C, Thomas BP , Tiao G, Tusie-Luna MT, Weisburd B, Won HH, Yu\nD, Altshuler DM, Ardissino D, Boehnke M, Danesh J, Donnelly S, Elosua R, Florez JC, Gabriel\nSB, Getz G, Glatt SJ, Hultman CM, Kathiresan S, Laakso M, McCarroll S, McCarthy MI,\nMcGovern D, McPherson R, Neale BM, Palotie A, Purcell SM, Saleheen D, Scharf JM, Sklar P ,\nSullivan PF, Tuomilehto J, Tsuang MT, Watkins HC, Wilson JG, Daly MJ, MacArthur DG; Exome\nAggregation Consortium. Analysis of protein-coding genetic variation in 60,706 humans.\nNature. 2016 Aug 18;536(7616):285-91. doi: 10.1038/nature19057. PMID: 27535533; PMCID:\nPMC5018207. https://www.nature.com/articles/nature19057\nLu W, Gauthier LD, Poterba T, Giacopuzzi E, Goodrich JK, Stevens CR, King D, Daly MJ, Neale\nBM, Karczewski KJ. CHARR e ﬃciently estimates contamination from DNA sequencing data.\nAm J Hum Genet. 2023 Dec 7;110(12):2068-2076. doi: 10.1016/j.ajhg.2023.10.011. PMID:\n38000370; PMCID: PMC10327099.\nMichael F. Lin, Ohad Rodeh, John Penn, Xiaodong Bai, Je ﬀrey G. Reid, Olga Krasheninina,\nWilliam J. Salerno. GLnexus: joint variant calling for large cohort sequencing. bioRxiv 343970;\nhttps://www.biorxiv .org/content/10.1101/343970v1.\nOle Schulz-Trieglaﬀ, Andrew Lee, Zhuoyi Huang, and Cobus De Beer. 2023-04-18. “Genotyping\nvariants at population scale using DRAGEN gVCF Genotyper”. Retrieved 2023-11-20.\nhttps://www.illumina.com/science/genomics-research/articles/gVCF-Genotyper.html.\nhttps://web.archive.org/web/20230529012612/https://www.illumina.com/science/genomics\n-research/articles/gVCF-Genotyper.html.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint \n\nSun KY , Bai X, Chen S, Bao S, Kapoor M, Backman J, Joseph T, Maxwell E, Mitra G, Gorovits A,\nMansﬁeld A, Boutkov B, Gokhale S, Habegger L, Marcketta A, Locke A, Kessler MD, Sharma D,\nStaples J, Bovijn J, Gelfman S, Gioia AD, Rajagopal V , Lopez A, Varela JR, Alegre J, Berumen J,\nTapia-Conyer R, Kuri-Morales P , Torres J, Emberson J, Collins R; Regeneron Genetics Center;\nRGC-ME Cohort Partners; Cantor M, Thornton T, Kang HM, Overton J, Shuldiner AR, Cremona\nML, Nafde M, Baras A, Abecasis G, Marchini J, Reid JG, Salerno W, Balasubramanian S. A deep\ncatalog of protein-coding variation in 985,830 individuals. bioRxiv [Preprint]. 2023 May\n10:2023.05.09.539329. doi: 10.1101/2023.05.09.539329. PMID: 37214792; PMCID: PMC10197621.\nhttps://pubmed.ncbi.nlm.nih.gov/37214792/.\nYun T, Li H, Chang PC, Lin MF, Carroll A, McLean CY . Accurate, scalable cohort variant calls\nusing DeepVariant and GLnexus. Bioinformatics. 2021 Apr 5;36(24):5582-5589. doi:\n10.1093/bioinformatics/btaa1081. PMID: 33399819; PMCID: PMC8023681.\nhttps://academic.oup.com/bioinformatics/article/36/24/5582/6064144.\nKushal Datta, Karthik Gururaj, Mishali Naik, Paolo Narvaez, Ming Rutar, “GenomicsDB:\nStoring Genome Data as Sparse Columnar Arrays”. Intel White Paper. 2017\nhttps://www.intel.com/content/dam/www/public/us/en/documents/white-papers/genom\nics-storing-genome-data-paper.pdf\nPapadopoulos, Stavros & Datta, Kushal & Madden, Samuel & Mattson, Tim. (2016). The TileDB\narray data storage manager. Proceedings of the VLDB Endowment. 10. 349-360.\n10.14778/3025111.3025117.\nhttps://people.csail.mit.edu/stavrosp/papers/vldb2017/VLDB17_TileDB.pdf\nXunjieli. 2019-07-10. GitHub comment.\nhttps://web.archive.org/web/20201024170035/https://github.com/dnanexus-rnd/GLnexus/i\nssues/173. Original: https://github.com/dnanexus-rnd/GLnexus/issues/173#issue-466548036.\n.CC-BY 4.0 International licenseavailable under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made \nThe copyright holder for this preprintthis version posted January 10, 2024. ; https://doi.org/10.1101/2024.01.09.574205doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}