tomwhite
commented
4 years ago I spent some time this week prototyping what this would look like.
The first piece is a function to read a (bgzip compressed) VCF file using cyvcf2 and convert it into sgkit xarray format. It does this by reading a chunk of N variants (N=10,000 for my experiment), creating an xarray dataset for those variants, then calling to_zarr on the dataset with append_dim="variants" to append to a Zarr store on disk. This is very similar to the way scikit-allel's vcf_to_zarr function is structured.
So far, everything is serial. To run in parallel, there needs to be a way to i) divide the VCF(s) into parts to be processed, and ii) merge the resulting Zarr stores.
For i), since VCF files are commonly tabix indexed, we could use genomic regions to partition the VCF(s) into roughly equal parts, each of which would be processed independently to create a Zarr store. We'd have to come up with some heuristics for doing a decent job on partitioning (e.g. like this).
Then, for ii), the problem is that each Zarr store will have a different number of variants, and even though the chunk sizes are the same, they need to be rechunked in the variant dimension.
For the sake of example, if the chunk size was 10, and the first Zarr store had (3, ) variants, and the second had (10, 6) variants, then they would be merged and rechunked as (3+7, 3+6) = (10, 9) variants.
This can be achieved in Xarray with
ds = xr.concat(datasets, dim="variants")
ds = ds.chunk({"variants": chunk_size})
ds.to_zarr(output, mode="w")However, there may be some caveats with running at scale, as discussed in https://github.com/dask/dask/issues/5105 (by @alimanfoo!), which may need work (rechunker may be an option too). I haven't looked at performance, but VCF parsing should be very fast with cyvcf2, since it is a thin interface to htslib (see the paper: https://academic.oup.com/bioinformatics/article/33/12/1867/2971439).
This two step-scheme should scale well (using Dask), since it is basically two embarrassingly parallel steps, mediated by a set of intermediate Zarr stores.
There are some build issues we'd probably want to address too:
- Installation requires installing htslib. I compiled it myself (on Mac), but using bioconda is another option, although we'd probably want to build wheels
- I successfully read a VCF file on GCP (using a
gs://orhttps://URL). I didn't try to authenticate, so that may need some work, or documentation at least, in the future
jeromekelleher
hammer
eric-czech
We can read data from VCF Zarr files with #40.
Ultimately we will need to read data from VCF files directly.
https://github.com/brentp/cyvcf2 looks to be a good option for this purpose.
Understanding our options here will also help us reach a conclusion on #65.