How to read in VCF files for use in tsinfer?

[hyanwong]

How do I import a VCF file into tsinfer?

[hyanwong]

Replying to my own question, there's code provided in the tutorial for reading VCF files, including reading in some of the metadata in the VCF (details further down that page). If you have a very large VCF file and it is taking a long time, then it is possible to read chunks of your VCF in parallel: see #277 (comment) for example code.

Note that in the longer term, we may be moving tsinfer to using sgkit as an import framework: when this happens this answer may become obsolete.

Also note that you if you are importing VCFs into tsinfer, you may wish to investigate how to infer ancestral states properly.

[npb596]

Hi Yan,

Have there been any updates to this? I ran a script with the exact code given in that tutorial link (except with file path changed to my VCF) on a VCF generated by SLiM (where I manually input a sequence length in header) and receive this error:

/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py:458: FutureWarning: The LMDBStore is deprecated and will be removed in a Zarr-Python version 3, see zarr-developers/zarr-python#1274 for more information.
return zarr.LMDBStore(self.path, subdir=False, map_size=map_size)
Traceback (most recent call last):
File "/mnt/shared/projects/uosa/Nick_Bailey/vcf_tskit.py", line 51, in
with tsinfer.SampleData(
^^^^^^^^^^^^^^^^^^^
File "/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py", line 952, in init
super().init(**kwargs)
File "/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py", line 384, in init
store = self._new_lmdb_store(max_file_size)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py", line 458, in _new_lmdb_store
return zarr.LMDBStore(self.path, subdir=False, map_size=map_size)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/zarr/storage.py", line 2331, in init
self.db = lmdb.open(path, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^
lmdb.MemoryError: /mnt/shared/projects/uosa/Nick_Bailey/BS_1_p1.samples: Cannot allocate memory

I can provide the exact script and VCF if needed but again both shouldn't have any substantive deviation from what's assumed in tutorial.

Cheers,
Nick

[benjeffery]

Add import tqdm to the top of your script.

[npb596]

That gave me this error:

TypeError: 'module' object is not callable. Did you mean: 'tqdm.tqdm(...)'?

and then I decided I don't really need/want a progressbar so I commented out the relevant lines. I now have this:

/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py:458: FutureWarning: The LMDBStore is deprecated and will be removed in a Zarr-Python version 3, see zarr-developers/zarr-python#1274 for more information.
return zarr.LMDBStore(self.path, subdir=False, map_size=map_size)
/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py:428: FutureWarning: The LMDBStore is deprecated and will be removed in a Zarr-Python version 3, see zarr-developers/zarr-python#1274 for more information.
store = zarr.LMDBStore(
Traceback (most recent call last):
File "/mnt/shared/projects/uosa/Nick_Bailey/vcf_tskit.py", line 54, in
add_diploid_sites(vcf, samples)
File "/mnt/shared/projects/uosa/Nick_Bailey/vcf_tskit.py", line 39, in add_diploid_sites
samples.add_site(pos, genotypes, alleles, ancestral_allele=ancestral_allele)
File "/mnt/apps/users/nbailey/conda/lib/python3.12/site-packages/tsinfer/formats.py", line 1832, in add_site
raise ValueError("Alleles must be distinct")
ValueError: Alleles must be distinct

[benjeffery]

Ah, sorry should have been from tqdm import tqdm.

That error means your alleles array for the given site has repeats in. Each allele needs to be unique.

[npb596]

Hmm, looking at the VCF I have a guess what's going on there. Some sites have the same ALT allele repeated twice like the A here (in order of CHROM, POS, ID, REF ALT):

1 892 . T A,A 1000 PASS

I guess it had to do with how SLiM generates the files. I could filter these out, does that seem like the right idea?

[npb596]

Yes that seemed to be the issue, along with other similarly odd sites. A clunky but viable solution right now is:

cat BS_1_p1.vcf | grep -v -E "A\s*\sA" | grep -v -E "C\s*\sC" | grep -v -E "G\s*\sG" | grep -v -E "T\s*\sT" | grep -v "[ACGT],[ACGT]" > test.vcf

where the first four grep conditions exclude sites where the ref and alt are the same allele and the final one excludes all multi-allelic sites, because it wasn't just the multi-allelic sites with the same two alt alleles that caused the issue. The output after that seems to be a valid samples file as per the tutorial. Thanks for all the quick assistance!