How to make a small zarr file from a large VCF?

[hyanwong]

For testing, it can be useful to subset the zarr file that is created from a huge VCF. Is there a way to restrict the samples and/or sites in the VCF to a smaller workable subset, or is the approved way to subset the zarr file somehow later? If the latter, is there any example code to show how to do this?

I don't want to do this on-the-fly during analysis using masks, because the zarr file itself it too large to place in permanent storage (larger files in personal storage areas get deleted).

[jeromekelleher]

You can use the -max-variant-chunks option to encode to limit the number of variants in the output Zarr. You do have to explode the full VCF though.

I don't know of any code to do the post-hoc subsetting. This would be non-trivial if you want to exclude all sites that are fixed in the output and so on.

In general the recommended pattern would be to adapt the downstream software to work with a subset of the Zarr rather than always expect to work with the entire dataset.

[hyanwong]

Thanks. In similar cases, I suspect it could be helpful to pass around a smaller copy of the .vcz file, subset to a handful of samples. Passing around the whole .vcz file is prohibitive, but it seems wasteful to have to do this by (slow) subsetting of the VCF file first.

Mind you, this is a somewhat limited use-case, so perhaps it's not worth providing an API.

[tomwhite]

I've been playing with subsetting vcz files using bcftools-style regions, include/exclude expressions, and samples in this vcftools branch, following on from the work in sgkit-dev/vcztools#201.

The idea is to use xarray to do things like:

ds = xr.open_zarr(vcz, concat_characters=False)
root = zarr.open(vcz, mode="r")
ds = filter_regions(ds, root, regions="20:1230236-")
ds = filter_expressions(ds, include="FMT/DP>3")
ds = filter_samples(ds, samples="NA00002,NA00003")
# save ds...

Since it uses xarray it probably belongs in sgkit, with a dependency on vcztools?

[tomwhite]

Since it uses xarray it probably belongs in sgkit, with a dependency on vcztools?

See sgkit-dev/sgkit#1330