bcf conversion memory requirement

[alxsimon]

Hi,
I am encountering an issue at the first step of using sgkit.
I want to convert a bcf to zarr.
The bcf is 45Gb containing 18 samples for 268031806 sites.

First trying directly with vcf_to_zarr this fills my 512Go memory and is killed.

vcf_to_zarr(
	in_vcf_wgs,
	out_zarr,
	ploidy=4,
	mixed_ploidy=True, # code missing ploidy as -2
	fields=[
		'INFO/DP',
		'FORMAT/DP',
		'FORMAT/AD',
		'FORMAT/VAF',
	],
)

I imagine this is due to the rechunking step being memory intensive.

So I thought I understood that the way to go was to use partition_into_regions followed by vcf_to_zarrs to avoid having the full dataset into memory.

regions_wgs = partition_into_regions(in_vcf_wgs, target_part_size='500MB')
urls_wgs = vcf_to_zarrs(
	in_vcf_wgs,
	outdir,
	regions=regions_wgs,
	ploidy=4,
	mixed_ploidy=True, # code missing ploidy as -2
	fields=[
		'INFO/DP',
		'FORMAT/DP',
		'FORMAT/AD',
		'FORMAT/VAF',
	],
)

This is running right now but already taking more than 200Go and growing.
Could you explain to me what is going on under the hood?

What am I missing here? Is there a way to reduce the memory requirement of the conversion? What is the best way to obtain an sgkit formatted dataset from vcf for large data?

Thanks a lot!

[tomwhite]

Hi @alxsimon - thanks for raising this issue.

Could you explain to me what is going on under the hood?

The vcf_to_zarr function partitions the VCF (or BCF in this case) into a set of contiguous regions that cover the while file,
and then writes an intermediate Zarr file for each partition, in parallel (using Dask). The intermediate Zarr files are concatenated, rechunked, then written to the final output (again using Dask). We have had problems in the past with this rechunking step running out of memory, which is what seems to be happening in this case.

The underyling Dask issue has never been fixed (see dask/dask#6745), but we do have some (internal) code in sgkit that avoids this memory problem:https://github.com/pystatgen/sgkit/blob/718fe3b58b3da2d231bbfa6330a88bf209068c92/sgkit/io/utils.py#L111-L148

Unfortunately it's not hooked up to vcf_to_zarr - which is an oversight on my part. (It's used in vcfzarr_to_zarr which reads from scikit-allel's VCF Zarr format, and which is where I hit the Dask issue.)

I'm going to take a look at hooking up that code to vcf_to_zarr. But in the meantime, as a workaround to run the conversion sequentially, you could set target_part_size=None in vcf_to_zarr:

vcf_to_zarr(in_vcf_wgs, out_zarr, target_part_size=None, ...)

This will avoid writing the intermediate Zarr files, and hence shouldn't hit any memory issues due to rechunking, but it will not be able to take advantage of multiple CPUs on your machine.

Hope all of that makes sense! Please let us know how you get on or if you have any more questions.

BTW the reason that your second snippet of code didn't help is that vcf_to_zarr is equivalent to partition_into_regions with a target_part_size of auto (which is 100MB) followed by vcf_to_zarrs. So it's effectively the same, but with 500MB in your case.

[alxsimon]

Thanks a lot for the clear explanation!
I'll try the target_part_size option, but it would indeed be great to be able to take advantage of multiple CPUs.