Document different ways of ancestral state specification

[hyanwong]

@AprilWei001 was finding it harder to use the VariantData class than the original SampleData class. Some of this was based around the extra step of providing ancestral states as a FASTA file. However, it should be just as easy to provide either a values in the VCF file (e.g. the info_AA field) or simply use the VCF reference sequence as the ancestral allele (although not entirely recommended).

We should make sure we document this in https://tskit.dev/tsinfer/docs/latest/usage.html#variantdata-and-ancestral-alleles.

[AprilWei001]

Agree that this is doable by rewriting the VCF files (which is what we do now), but it partially reduced the whole benefit of using VCF.zarr. Because we need to have extra I/Os related to VCF before running tsinfer. It also seems that we can only convert VCF.gz to Zarr, not VCF directly.
In our experience, tsinfer is somewhat sensitive to the ancestral derived state.

[jeromekelleher]

Can you outline your workflow a little here please @AprilWei001, so we can understand where the bottlenecks are and how we can alleviate?

It's true that bio2zarr requires indexed, bgzipped VCF as this is necessary to support parallel reading of the VCF .

[AprilWei001]

Previously:
For each line in VCF file:
Add variant to sample data (if not multi-allelic), and choose the ancestral allele to be the REF (this is how simulated data is exported, and how our polarization works)

Currently:
bgzip the VCF file
use tabix to index it
vcf2zarr convert the file
Load the ZARR file in Python and construct an ancestral state array using the REF allele via ancestral_state = vcf_zarr["variant_allele"][:, 0]

[jeromekelleher]

This is more steps all right, but it should be substantially faster overall if you're using several cores. Going from uncompressed VCF isn't a common use-case (nearly everyone bgzips it), so I'm curious about where the simulated data is coming from. Is it in tskit format? If so, we have a tskit2zarr program in the pipelines: sgkit-dev/bio2zarr#232

[hyanwong]

It looks to me like the extra steps are less to do with the ancestral state specification, and more in the gzip and tabix steps. Indeed, I have previously found it a real pain that there isn't an easy way to create index files from vcfs without installing all the SAMtools stuff (including tabix).

Going from uncompressed VCF isn't a common use-case

I'm not sure that's entirely true. That might be true for real datasets, but there are lots of programs that simulate VCFs and simply output uncompressed files (but not tree sequences). I have previously wanted to do all of this without having to resort to installing extra software, and it's been problematic, especially e.g. on windows where the bioinformatic tools are less mature.

[jeromekelleher]

We could support uncompressed VCF in bio2zarr if there's a real need - it wouldn't be that hard, it's just that the VCF would have to be read sequentially. Please do open an issue there if you need it.

However, I'm not sure that bio2zarr will ever work on Windows (because of the dependencies like cyvcf2).

[hyanwong]

FYI I found it was less to do with the compression and more to do with the extra software required to generate the index file (tabix) that hindered me. It also made it harder to set up a self-contained conda env to run tests on a cluster.

Requiring an index file means that you can't just template everything in python, but have to start fiddling with running external commands, which slows development in my experience.

[hyanwong]

However, I'm not sure that bio2zarr will ever work on Windows (because of the dependencies like cyvcf2).

When teaching people to use tsinfer, I found that several students (esp from developing counties) were trying to follow instructions for VariantData use on their windows laptops, and having problems with installation etc. I didn't have this problem when going via the SampleData format.

[jeromekelleher]

That can't be helped I'm afraid, we can't make vital dependencies work ourselves. WSL works pretty well these days, I think that's a good workaround for most people.

[benjeffery]

Do we know what is blocking cyvcf2 on windows? I assume it is possible? I had a quick search but couldn't see what it was.

[jeromekelleher]

I don't know TBH, I don't think it's worth spending a lot of time on. WSL is good, people should just use it if they really want to do inference from real data.

[hyanwong]

Oh yes, I'm not really concerned about inference from real data. It's the pathway for doing inference from simulated data that still feels less-than-smooth. I think that's where @AprilWei001 is coming from too? I suspect that tsinfer will be run more times on simulated data than inferred data, TBH (mainly because even analysis of real data is likely to also include a lot of testing on simulations too).

[jeromekelleher]

Whether bio2zarr supports windows is a side issue here, let's take it somewhere else if people really want this.

[jeromekelleher]

FWIW bio2zarr will support uncompressed/unindexed VCF soon: sgkit-dev/bio2zarr#337

I assume that the requirement that each VCF files contains records for a single contig isn't a deal-breaker on these workflows?

[hyanwong]

@AprilWei001 : the most recent release of bio2zarr allows you to use unindexed VCFs, and I've added the description of using the REF as the ancestral state to the doves, so I think the number of steps is now basically same for SampleData and VariantData methods:

SampleData

For each line in VCF file:
Add variant to sample data (if not multi-allelic), and choose the ancestral allele to be the REF

VariantData

run vcf2zarr, then

vcfzarr = zarr.load(vcz_fn)
inferre_ts = tsinfer.infer(
    tsinfer.VariantData(vcfzarr, ancestral_allele= vcfzarr["variant_allele"][:, 0])

If that works for you, can I close this issue?