How to write a Dataset

How to write a Dataset#

Depending on your needs, you’ll potentially need to normalize variants and prepare a mapping from sample names to BigWig files before you can write a Dataset.

Preparing input data#

Normalizing variants#

Before passing variants to GenVarLoader, they must be:

  • left-aligned

  • bi-allelic

  • atomic (no MNPs or compound MNP-indels)

  • for VCFs, indexed with bcftools index or tabix

In general, any VCF can be preprocessed to meet these requirements using bcftools norm, for example:

bcftools norm -f $reference \
    -a --atom-overlaps . \
    -m -any --multi-overlaps . \
    -O b -o $out $in

Alternatively, if your data is already in the PGEN format you can use PLINK 2.0:

plink2 --make-bpgen --pfile $input --out $intermediate
plink2 --make-pgen --normalize --ref-from-fa --fa reference.fa --bpfile $intermediate --out $normalized

See the PLINK 2.0 documentation for details.

BigWig table#

To process BigWig data, GenVarLoader needs a mapping from sample names to BigWig file paths. This can be provided either as a dictionary using BigWig() or a table using BigWig.from_table(). If using a table, it must at least contain the columns sample and path, for example:

sample

path

Aang

/data/bigwigs/aang.bw

Katara

/data/bigwigs/katara.bw

Sokka

/data/bigwigs/sokka.bw

Regions of interest (ROIs)#

Whether you’re working with variants, BigWigs, or both, you will need regions of interest specified in a BED-like table. This means your table must have the columns chrom, chromStart, and chromEnd as 0-based coordinates. You can also include a strand column consisting of +, -, and . (treated as +) which will let you choose to reverse and/or reverse complement data on negative strands. Any extra columns will be kept in the table and exist in Dataset.regions after writing. An example of a BED-like table:

split

chrom

chromStart

chromEnd

strand

train

chr1

0

5

-

val

chr2

3

8

+

test

chr3

20

21

.

Writing data#

Once your data is prepared, you can use gvl.write() to convert the data for regions of interest to a format that gvl.Dataset can open. You can include a single source of variants and/or any number of BigWig files. Some examples:

import genvarloader as gvl

gvl.write(
    path="1000_genomes_haplotypes.gvl",
    bed="tiling_windows.bed",
    variants="all_chroms.bcf",
    # OR variants='all_chroms.pgen',
)

This dataset would have haplotypes available for all samples in all_chroms.bcf.

gvl.write(
    path="1000_genomes_lncRNA.gvl",
    bed="lncRNA.bed",  # can be varying length regions
    variants="all_chroms.bcf",
    tracks=[
        gvl.BigWigs.from_table("pos", "pos_strands.tsv"),
        gvl.BigWigs.from_table("neg", "pos_strands.tsv"),
    ],
)

This dataset would have both haplotypes and two tracks (pos and neg) available for samples that exist in both all_chroms.bcf and the BigWig tables (i.e. gvl.write() performs an inner join on samples).

Variants from a genoray sparse store (.svar / .svar2)#

Besides BCF/VCF and PGEN, variants= also accepts a genoray sparse columnar variant store — either the original .svar format or the newer .svar2 format. Build one from a normalized VCF/BCF with genoray:

from genoray import VCF, SparseVar, SparseVar2

# .svar (SVAR1): a VCF reader + a memory budget
SparseVar.from_vcf("all_chroms.svar", VCF("normed.bcf"), max_mem="4g")

# .svar2 (SVAR2): the VCF/BCF path + a reference FASTA (or no_reference=True)
SparseVar2.from_vcf("all_chroms.svar2", "normed.bcf", reference="ref.fa")

Then pass the resulting store to gvl.write:

gvl.write(
    path="1000_genomes_haplotypes.gvl",
    bed="tiling_windows.bed",
    variants="all_chroms.svar2",  # or "all_chroms.svar", or a SparseVar/SparseVar2 instance
)

Both formats store a back-reference in the dataset’s metadata.json instead of duplicating per-variant arrays, so the source store must remain accessible when the dataset is later opened with gvl.Dataset.open() (override its location with svar=/svar2= if it has moved).

.svar2 additionally produces a write-time cache under <path>/genotypes/svar2_ranges/ and reads back through an all-Rust, read-bound path with no interval-search-tree build and no dense-union rebuild per read — see the FAQ for the read-path and on-disk-size tradeoffs, and the format reference for the on-disk layout, including the size formula. This cache is small at cohort scale: ~504 MB for the All of Us chr22 grid, since it stores 28 bytes only for the (region, sample, ploid) windows that hold a variant. max_mem bounds the genoray chunk stream gvl.write reads while producing the cache, but not the per-contig entry accumulator: the per-chunk (region int32, cell int32, 24-byte entry) blocks and the merged (cell, entry) output it is sorted into live at the same time, peaking at roughly 60 bytes per entry on the largest contig — about 1.1 GB at All of Us chr22 and ~3.4 GB at chr19. The permanent range cache on disk is governed by disk space, not max_mem (see the format reference). .svar2 currently has a Phase-1 scope: a handful of output combinations (annotated haplotypes, min_af/max_af, spliced variant-window/track outputs, etc.) aren’t wired yet and raise NotImplementedError — see the genvarloader skill or the format reference for the full list. Haplotype and variants output support splicing and var_filter="exonic".

Sample order#

Whether you pass samples= or leave it None, gvl.write stores the selection in lexicographic order, and Dataset.samples reports that same order.

Warning

Lexicographic order is not numeric order. A cohort with integer-like IDs of mixed digit counts sorts as "1000" < "999", whereas the phenotype table you want to join against is usually in numeric order. Align external tables to Dataset.samples by name, never by position.

Reusing a variant index across cohorts#

If several cohorts are sample subsets of one parent PGEN, do not pre-split the PGEN with plink2 --keep. Point gvl.write at the parent and pass samples=:

gvl.write("cohort_A.gvl", bed, "parent.pgen", samples=cohort_A_samples)
gvl.write("cohort_B.gvl", bed, "parent.pgen", samples=cohort_B_samples)

genoray caches the parent PGEN’s variant index on disk (.pvar.gvi, keyed by the .pvar’s mtime); samples= only subsets at the genotype-extraction level and never touches that index, so a second gvl.write against the same parent path reuses the cached index instead of rescanning the PVAR — and gvl hardlinks that cached index file into each dataset’s genotypes/variants.arrow. Splitting the PGEN first (plink2 --keep) produces a distinct .pvar per cohort, which forces the expensive per-cohort index rebuild.

Merging datasets#

gvl.concat() merges datasets that were written from one shared variant source (the same PGEN/VCF variant table, or the same .svar/.svar2 store — merging datasets built from different variant sources raises), along either axis:

gvl.concat("merged.gvl", ["chr1.gvl", "chr2.gvl"], axis="regions")
gvl.concat("merged.gvl", ["cohortA.gvl", "cohortB.gvl"], axis="samples")

axis="regions" requires identical samples in identical order across every input, and concatenates their regions. axis="samples" requires identical regions across every input and disjoint sample sets, and merges the samples into sorted order. Both axes merge tracks and annotation tracks alongside the genotypes, and require at least two input datasets.

gvl.concat streams data at the byte level rather than re-deriving anything from the source variants, so its cost is dominated by I/O rather than computation. As an order-of-magnitude expectation (not a measured figure — this hasn’t been benchmarked), it moves roughly the full size of the merged dataset at sequential-IO speed, so a merge of a very large dataset (e.g. on the order of a terabyte) could plausibly take hours rather than minutes. It’s worth doing against the alternative of re-extracting genotypes from scratch, not as a routine step.

gvl.concat’s planning is streaming: it derives the merge plan from the sorted merge order rather than building a map over every (region, sample, ploid) slot, so its memory stays proportional to regions plus samples rather than to their product. The merged offsets.npy it writes is still one int64 per merged (region, sample[, ploid]) slot — about 16 GB for per-sample tracks (no ploidy axis) or 32 GB for diploid genotypes, at a 3,734-region by 535,662-sample grid — and that is a property of the ragged on-disk format, not of the merge. Peak resident memory during the copy is roughly double that figure: the copy also holds each input’s own offsets array for the duration, and since the inputs partition the merged grid, their combined size is about the same as the merged array itself — so about 32 GB for per-sample tracks or 64 GB for diploid genotypes, at the same grid.

For contig-sharded .svar2 workflows there is a cheaper path at the variant-store layer: genoray.SparseVar2.concat(output, sources, mode="copy") merges disjoint-contig .svar2 stores that share identical samples/ploidy/fields, then a single gvl.write call over the merged store only needs to populate .svar2’s small per-(region, sample, ploid) cache ranges rather than touch the bulk variant data. This is a different operation from gvl.concat, which merges already written .gvl dataset directories — SparseVar2.concat merges the upstream .svar2 stores before gvl.write ever runs.