Reading VCZ#
This page is a worked introduction to reading VCF Zarr (VCZ) datasets from
Python. Every code cell below runs against a small example dataset,
data/sample.vcz.zip (9 variants across 3 samples, plus a null sample at index
1 used to illustrate masking; see Null samples).
For exact signatures of everything used here, see the API reference.
Opening a dataset#
vcztools.open_zarr() is the recommended entry point. It returns a
zarr.Group and dispatches to one of four storage backends; see
Storage backends for the full description of each. A local directory
or .zip archive needs no explicit backend:
import vcztools
root = vcztools.open_zarr("data/sample.vcz.zip")
root
<Group zip://data/sample.vcz.zip>
Remote URLs require selecting a backend, and storage_options is forwarded to
it (the fsspec filesystem, obstore.store.from_url, or the chosen Icechunk
storage constructor):
root = vcztools.open_zarr(
"s3://bucket/sample.vcz",
backend_storage="fsspec",
storage_options={"anon": True},
)
An already-built zarr.Group or zarr.abc.store.Store is
accepted and passed through unchanged, which is useful when you’ve configured a
store yourself:
import zarr
store = zarr.storage.MemoryStore()
# ... populate ``store`` ...
root = vcztools.open_zarr(store)
Creating a reader#
vcztools.VczReader wraps an opened group and provides metadata and
variant iteration. Use it as a context manager so resources are released
deterministically:
with vcztools.VczReader(root) as reader:
print(reader.num_variants, "variants")
print(reader.num_samples, "samples")
9 variants
4 samples
Displaying a reader in a notebook renders a summary of the store and the reader’s current state — its dimensions, null-sample count, sample and variant selection, any configured filter, and the available fields:
reader = vcztools.VczReader(root)
reader
| Store | ZipStore('zip://data/sample.vcz.zip') |
|---|---|
| Source | bio2zarr-0.2.0 |
| Variants | 9 (chunk 4) |
| Samples | 4 (chunk 4) |
| Null samples | 1 |
| Sample selection | all non-null (default) |
| Variant selection | all (default) |
| Filter | none |
| Contigs | 3 |
| Filters | 3 |
| Readahead | 32 workers, 256.0 MiB |
| Fields | 27 stored, 8 virtual |
Fields (31)
| name | dtype | shape | dims |
|---|---|---|---|
| call_DP | int8 | (9, 4) | variants, samples |
| call_GQ | int8 | (9, 4) | variants, samples |
| call_HQ | int8 | (9, 4, 2) | variants, samples, FORMAT_HQ_dim |
| call_genotype | int8 | (9, 4, 2) | variants, samples, ploidy |
| call_genotype_mask | bool | (9, 4, 2) | variants, samples, ploidy |
| call_genotype_phased | bool | (9, 4) | variants, samples |
| contig_id | StringDType() | (3,) | contigs |
| filter_description | StringDType() | (3,) | filters |
| filter_id | StringDType() | (3,) | filters |
| region_index | int32 | (4, 6) | region_index_values, region_index_fields |
| sample_id | StringDType() | (4,) | samples |
| variant_AA | StringDType() | (9,) | variants |
| variant_AC | int8 | (9, 2) | variants, INFO_AC_dim |
| variant_AF | float32 | (9, 2) | variants, INFO_AF_dim |
| variant_AN | int8 | (9,) | variants |
| variant_DB | bool | (9,) | variants |
| variant_DP | int8 | (9,) | variants |
| variant_H2 | bool | (9,) | variants |
| variant_NS | int8 | (9,) | variants |
| variant_allele | StringDType() | (9, 4) | variants, alleles |
| variant_contig | int8 | (9,) | variants |
| variant_filter | bool | (9, 3) | variants, filters |
| variant_id | StringDType() | (9,) | variants |
| variant_id_mask | bool | (9,) | variants |
| variant_length | int8 | (9,) | variants |
| variant_position | int32 | (9,) | variants |
| variant_quality | float32 | (9,) | variants |
| variant_F_MISSING (virtual) | float64 | (9,) | variants |
| variant_N_ALT (virtual) | int64 | (9,) | variants |
| variant_N_MISSING (virtual) | int64 | (9,) | variants |
| variant_index (virtual) | int64 | (9,) | variants |
The reader starts with every sample and every variant selected. The remaining sections each open a fresh reader so their selections stay independent.
Inspecting metadata#
reader = vcztools.VczReader(root)
print("sample_ids:", reader.sample_ids)
print("contig_ids:", reader.contig_ids)
sample_ids: ['NA00001' 'NA00002' 'NA00003']
contig_ids: ['19' '20' 'X']
sample_ids lists only the real samples, so it is shorter than num_samples
(printed above), which counts every slot in the store including null samples;
see Null samples.
field_names lists the arrays stored in the dataset, while
virtual_field_names lists fields computed on demand (allele counts and
frequencies, missingness, and so on) that are never emitted unless you request
them by name:
print("stored fields:", sorted(reader.field_names))
print()
print("virtual fields:", sorted(reader.virtual_field_names))
stored fields: ['call_DP', 'call_GQ', 'call_HQ', 'call_genotype', 'call_genotype_mask', 'call_genotype_phased', 'contig_id', 'filter_description', 'filter_id', 'region_index', 'sample_id', 'variant_AA', 'variant_AC', 'variant_AF', 'variant_AN', 'variant_DB', 'variant_DP', 'variant_H2', 'variant_NS', 'variant_allele', 'variant_contig', 'variant_filter', 'variant_id', 'variant_id_mask', 'variant_length', 'variant_position', 'variant_quality']
virtual fields: ['variant_AC', 'variant_AF', 'variant_AN', 'variant_F_MISSING', 'variant_NS', 'variant_N_ALT', 'variant_N_MISSING', 'variant_index']
Field metadata#
get_field_info() returns a vcztools.FieldInfo
snapshot describing a single field’s dtype, shape, dimensions and attributes:
info = reader.get_field_info("call_DP")
info
| dtype | int8 |
|---|---|
| shape | (9, 4) |
| dims | variants, samples |
| description | Read Depth |
Iterating variants#
variants() yields one dict per variant row. Pass
fields to restrict what is read; here we also request the virtual
variant_AF (allele frequency) field. We take the first three rows:
import itertools
rows = itertools.islice(
reader.variants(fields=["variant_position", "variant_allele", "variant_AF"]),
3,
)
for row in rows:
print(row["variant_position"], row["variant_allele"], row["variant_AF"])
111 ['A' 'C' '' ''] [nan nan]
112 ['A' 'G' '' ''] [nan nan]
14370 ['G' 'A' '' ''] [0.5 nan]
variant_chunks() yields one dict per variant chunk,
with each field as a NumPy array whose first axis is the variants in that chunk.
This is the efficient path for bulk processing:
for chunk in reader.variant_chunks(fields=["variant_position", "call_genotype"]):
positions = chunk["variant_position"]
genotypes = chunk["call_genotype"]
print(f"chunk of {len(positions)} variants, "
f"call_genotype shape {genotypes.shape}")
chunk of 4 variants, call_genotype shape (4, 3, 2)
chunk of 4 variants, call_genotype shape (4, 3, 2)
chunk of 1 variants, call_genotype shape (1, 3, 2)
Missing and fill values#
The arrays you get back are raw: they still contain the sentinel values VCZ uses to encode two distinct ideas, and you must interpret them yourself.
Missing — a value that is absent (“no data”). Detect it with
vcztools.is_missing().Fill (end-of-vector) — padding that makes a ragged field rectangular. A VCF INFO or FORMAT field can hold a variable number of values per variant or sample — for example one
AFper ALT allele, so a biallelic site has one value and a triallelic site has two. VCZ stores these in a fixed-width array and pads the unused tail of each row with a fill sentinel. Detect it withvcztools.is_fill(), or drop the trailing fill from a single 1-D vector withvcztools.trim_fill().
The sentinels differ by dtype (-1/-2 for integers, distinct
not-a-number values for floats, "."/"" for strings), so printing a raw
array often can’t tell the two apart by eye — the helpers can. The VCF and
query output paths already trim fill and render missing as .; these helpers
let you do the same when working with the arrays directly.
variant_AF is ragged over ALT alleles. Trimming the fill and then marking the
missing entries recovers the true per-variant vector:
reader = vcztools.VczReader(root)
fields = ["variant_position", "variant_allele", "variant_AF"]
for row in reader.variants(fields=fields):
af = row["variant_AF"]
trimmed = vcztools.trim_fill(af)
clean = [
"." if missing else round(float(value), 3)
for value, missing in zip(trimmed, vcztools.is_missing(trimmed))
]
alleles = [a for a in row["variant_allele"] if a != ""]
print(f"POS {row['variant_position']:>7} {'/'.join(alleles):8} "
f"raw {af} -> AF {clean}")
POS 111 A/C raw [nan nan] -> AF ['.', '.']
POS 112 A/G raw [nan nan] -> AF ['.', '.']
POS 14370 G/A raw [0.5 nan] -> AF [0.5]
POS 17330 T/A raw [0.017 nan] -> AF [0.017]
POS 1110696 A/G/T raw [0.333 0.667] -> AF [0.333, 0.667]
POS 1230237 T raw [nan nan] -> AF ['.', '.']
POS 1234567 G/GA/GAC raw [nan nan] -> AF ['.', '.']
POS 1235237 T raw [nan nan] -> AF ['.', '.']
POS 10 AC/A/ATG/C raw [nan nan] -> AF ['.', '.']
Multiallelic sites (e.g. A/G/T) keep both allele frequencies, biallelic sites
have their single value followed by one trimmed-away fill, and sites where AF
is absent come back as missing (.).
The same two sentinels appear in genotypes. Missing marks a no-call (. in
VCF), while fill pads a call shorter than the array’s ploidy — here one haploid
call in a diploid field:
reader = vcztools.VczReader(root)
missing_cells = 0
fill_cells = 0
fill_example = None
for row in reader.variants(fields=["variant_position", "call_genotype"]):
genotype = row["call_genotype"]
missing_cells += int(vcztools.is_missing(genotype).sum())
fill = vcztools.is_fill(genotype)
fill_cells += int(fill.sum())
if fill.any() and fill_example is None:
fill_example = (row["variant_position"], genotype)
print("missing (no-call) cells:", missing_cells)
print("fill (ploidy-padding) cells:", fill_cells)
position, genotype = fill_example
print(f"POS {position} has a haploid call padded with fill (-2):")
print(genotype)
missing (no-call) cells: 4
fill (ploidy-padding) cells: 1
POS 10 has a haploid call padded with fill (-2):
[[ 0 -2]
[ 0 1]
[ 0 2]]
Selecting samples and variants#
Call set_samples() with sample IDs before iterating
to restrict the sample selection. The output follows the order you request:
reader = vcztools.VczReader(root)
reader.set_samples(["NA00003", "NA00001"])
print("selected samples:", reader.sample_ids)
selected samples: ['NA00003' 'NA00001']
Pass complement=True to select every sample except those named, or
ignore_missing_samples=True to warn-and-drop unknown names instead of raising:
reader = vcztools.VczReader(root)
reader.set_samples(["NA00002"], complement=True)
print("selected samples:", reader.sample_ids)
selected samples: ['NA00001' 'NA00003']
set_variants() restricts the variant selection with a
sorted array of global variant indexes:
import numpy as np
reader = vcztools.VczReader(root)
reader.set_variants(np.array([0, 1, 2]))
positions = [row["variant_position"]
for row in reader.variants(fields=["variant_position"])]
print("positions:", positions)
positions: [np.int32(111), np.int32(112), np.int32(14370)]
Null samples#
Some VCZ datasets contain null (masked) samples: slots whose sample_id is
the empty string "". They must never be read, filtered on, or output. The
example dataset has one at index 1. The reader hides them automatically —
sample_ids and every iterated call_* array exclude null samples, and
sample-dependent virtual fields (AC/AN/AF/NS) ignore them.
The raw picture is available alongside the filtered view:
reader = vcztools.VczReader(root)
print("num_samples (raw, counts nulls):", reader.num_samples)
print("raw_sample_ids:", reader.raw_sample_ids)
print("sample_ids (nulls hidden):", reader.sample_ids)
num_samples (raw, counts nulls): 4
raw_sample_ids: ['NA00001' '' 'NA00002' 'NA00003']
sample_ids (nulls hidden): ['NA00001' 'NA00002' 'NA00003']
Iterating drops the null sample’s columns, so call_* arrays have one column
per real sample:
chunk = next(reader.variant_chunks(fields=["call_genotype"]))
print("call_genotype sample columns:", chunk["call_genotype"].shape[1])
call_genotype sample columns: 3
A null sample’s sample_id is the empty string, which is reserved and never
names a real sample. Passing it to set_samples()
raises a ValueError, so null slots are unreachable by ID:
reader = vcztools.VczReader(root)
try:
reader.set_samples([""])
except ValueError as exc:
print(exc)
The empty string is not a valid sample name: it is the reserved ID of a null sample, which is never selectable.
For the lower-level case where you select by raw integer index, use
set_sample_indexes(). Indexes are positions in the raw
sample_id array, and selecting a null index raises a ValueError:
reader = vcztools.VczReader(root)
try:
reader.set_sample_indexes([1]) # index 1 is the null sample
except ValueError as exc:
print(exc)
sample index refers to a null sample (sample_id == ''): [1]
Filtering#
vcztools.BcftoolsFilter compiles a bcftools -i/-e expression into a
filter. Construct it from the reader (so bare VCF names resolve against the
dataset’s fields), then attach it with
set_variant_filter():
reader = vcztools.VczReader(root)
variant_filter = vcztools.BcftoolsFilter(reader, include="QUAL>10")
reader.set_variant_filter(variant_filter)
for row in reader.variants(fields=["variant_position", "variant_quality"]):
print(row["variant_position"], row["variant_quality"])
14370 29.0
1110696 67.0
1230237 47.0
1234567 50.0
The fixed-width writers and encoders in Format conversion cannot
iterate a still-configured filter directly; resolve it into a fixed selection
first with materialise_variant_filter().