`write_vcf` returning position 0

[andrewkern]

hey team-

working with @mufernando on some stdpopsim coding we stumbled upon a bug with write_vcf. using the default position_transform value, write_vcf can return position 0 which isn't allowed under the vcf spec i believe. e.g., this tree sequence

print(ts.tables.sites)
╔═════╤════════╤═══════════════╤════════╗
║id   │position│ancestral_state│metadata║
╠═════╪════════╪═══════════════╪════════╣
║0    │       0│              C│        ║
║1    │     534│              C│        ║
║2    │     991│              T│        ║
║3    │    1188│              G│        ║
║4    │    1281│              C│        ║
║5    │    1301│              A│        ║

will return the following vcf

##fileformat=VCFv4.2
##source=tskit 0.5.5
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=1,length=5000000>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	tsk_0	tsk_1	tsk_2	tsk_3	tsk_4	tsk_5	tsk_6	tsk_7	tsk_8	tsk_9	tsk_10	tsk_11	tsk_12	tsk_13	tsk_14	tsk_15	tsk_16	tsk_17	tsk_18	tsk_19
1	0	0	C	G	.	PASS	.	GT	0|0	0|0	0|0	0|0	0|0	0|0	0|0	1|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0
1	534	1	C	G	.	PASS	.	GT	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|0	0|1	0|0	0|0	0|0	0|0	0|0
1
.
.
.

i'm guessing we should not have been using the default position_transform here but instead legacy? either way i think this is a bug

[petrelharp]

legacy will still give you 0s: https://github.com/tskit-dev/tskit/blob/main/python/tskit/vcf.py
... and it looks like a workaround is to pass position_transform = lambda x: 1 + np.round(x) to position_transform? But, this is very annoying - if we made this the default, then the 'position' in the VCF file would be off-by-one from the position in the tree sequence.

Alternatively we could special-case the situation where there's a mutation at position 0, ie, do position_transform = lambda x: np.fmax(1, np.round(x)) (and throw a warning if there is such a mutation)?

[andrewkern]

i think i would prefer the latter option.

[jeromekelleher]

I don't think there's a right answer here tbh. The VCF spec says it wants 1-based coords, but I'm not sure anything bad will happen if there is a zero position in there. Do we know if, say, bcftools will accept a vcf containing a single site with 0 coord?

It would definitely create more problems if we +1'd all the coords by default. For example, in sc2ts we store 1-based coordinates, so these would all be wrong if we started incrementing all coordinates by one in vcf output.

[andrewkern]

actually we found the bug because a downstream program threw an error on the VCF

[andrewkern]

we are creating a VCF in the stdpopsim/analysis2 pipeline here that was messing stuff up

https://github.com/popsim-consortium/analysis2/blob/eb16df28c44a1c3a04a1a5801172e806280f5c2f/workflows/sweep_simulate.snake#L399-L413

For the time being we can create that VCF with @petrelharp's suggested lambda function (i.e. position_transform = lambda x: np.fmax(1, np.round(x))) and then keep that project moving.

[mufernando]

I think a lot of downstream code that handles VCFs will assume it's 1-based. Should I at least add a warning to the docs with this potential workaround?

[petrelharp]

That's a good suggestion - provide the solution in the docs.

[bhaller]

Note that SLiM adds 1 to positions (which are 0-based in SLiM) when saving VCF, and subtracts 1 from positions when reading from VCF (the inverse operation). This means that a round-trip from SLiM to VCF and back preserves the correct information. It also means that, since tskit saves VCF with 0 for the position sometimes, SLiM errors out when reading some tskit-saved VCF files. Came up in this slim-discuss thread: https://groups.google.com/g/slim-discuss/c/cWwRxX2zNb0/m/CotUV5JyAAAJ. I don't know what the right course of action is for tskit at this point. :->

[bhaller]

The more I ponder this, though, the less it seems like a doc bug. It's saving out non-spec-compliant VCF files, and that is causing problems with at least two software packages downstream of it. Seems like maybe a fix is in order, probably doing the same thing that SLiM does (which would also be nice for interoperability, of course).

[mufernando]

I'm going to echo @bhaller's opinion.

I agree it is confusing to have the positions in the VCF be off by one when compared to positions in the tree sequence, but this is something everyone dealing with bioinformatics is aware of. I think it is a lot more confusing to have a VCF that is 0-based, specially if we want tree sequences to be adopted more widely.

I don't fully understand the ramifications of changing the default behavior, but now we have more evidence that the current behavior is messing with other people's downstream codes.

[petrelharp]

Correct me if I'm wrong, but I don't think that the issue is the usual 0-or-1-based indexing. If it was, then "gene X extends from position 1000 to position 1200" would mean two different things, that differ at each end by one bp. I think the issue is only that tskit allows a nucleotide at position 0, which is a much smaller deal.

[bhaller]

Correct me if I'm wrong, but I don't think that the issue is the usual 0-or-1-based indexing. If it was, then "gene X extends from position 1000 to position 1200" would mean two different things, that differ at each end by one bp. I think the issue is only that tskit allows a nucleotide at position 0, which is a much smaller deal.

I don't understand the distinction you're drawing. Can you elaborate?

[jeromekelleher]

I think the issue is only that tskit allows a nucleotide at position 0, which is a much smaller deal.

I agree with Peter here. Msprime generates sites that can have a 0 position, so you can say it uses 0-based coordinates. Sc2ts uses 1-based coordinates from real data, so it is 1-based. I'm pretty sure the ARGs we infer using tskit used 1-based coordinates. All use tskit.

Adding 1 to all coordinates is incorrect behaviour, as it changes the actual coordinates for everything. So, if we did this, then a VCF output from a sc2ts inferred ARG would have wrong coordinates.

I think this is more correctly seen as an issue with simulators generating 0 or 1 based coordinates, and tskit shouldn't have an opinion about that. I accept that VCF doesn't allow for a coordinate of 0 in the spec, so we should probably emit a warning if there is a site with this coord.

[jeromekelleher]

(As a follow up, we could easily add an option to msprime to disallow sites with position=0, which I'd be happy to discuss over there)