Tskit: Error with default tskit.write_vcf output in plink

[saurabhbelsare]

Unless alternatively specified, tskit's write_vcf function writes out default individual ids as tsk_0, tsk_1, … However, when I tried to input that vcf to plink, it would not accept “tsk_0” as it doesn’t accept “_0" at the end of an individual id. I fixed it simply by setting the individual ids for write_vcf as “tsk_0indv”, “tsk_1indv”… which plink is fine with.

This error stems from line 7846 of https://github.com/chrchang/plink-ng/blob/master/1.9/plink_data.c

Would it be possible to have a slightly modified default output format for tskit, or, if that wouldn't work, then having a note telling users to manually set individual ids if plink is part of their pipeline?

[benjeffery]

Thanks for the report @saurabhbelsare! Am I right in thinking using --id-delim=- in plink would fix this?
Also do you know if any other software uses this convention? We want things to work in the majority of cases, but we don't want to get into a game of VCF parser whack-a-mole!

[saurabhbelsare]

Hi @benjeffery,
Thanks for your reply. I tried using the --id-delim - option, but plink throws an error saying No '-' in sample ID, and the corresponding bed/bim/fam files are not generated. I think plink needs the sample ids to be in the form <familyID><delimiter><indvID>, and hence if it is tsk_0, it is going to look for whatever alternative delimiter we specify, and not find it in the name.

I'm not aware of any other software using this convention, in fact I only realized that plink had this issue today when I ran into this error. I agree though that with each parser having its own finicky requirements, it probably ends up being a whack-a-mole exercise. I think manually setting the sample ids to something else as I did is quite an easy work-around, so maybe just adding that as a note in the write_vcf documentation might be enough instead of having to change the existing default format? Thanks!

[petrelharp]

We might as well fix this - it would be easy - do you know what ID format would be more standard?

(but in the meantime: maybe --const-fid 0 would do it, which is default in plink 2.0?)

[jeromekelleher]

Thanks for the report @saurabhbelsare, this is really helpful information. Since we can specify the sample IDs manually, and people are probably depending on the current form, I don't think we want to change the default behaviour now.

A note in the documentation with a quick recipe for individual_names that'll work with plink is a good option I think.

[saurabhbelsare]

Thanks @petrelharp, @jeromekelleher. Just something as simple as this worked for me

indv_names = ['tsk_'+str(i)+'indv' for i in range(sample_size)]

with open('output.vcf', 'w') as vcf_file:
    tree_sequence.write_vcf(vcf_file, ploidy=1, individual_names=indv_names)

[jeromekelleher]

Good, looks like we all agree then. @saurabhbelsare, it would be super helpful if you could open a PR to add this to the documentation for write_vcf.

I think a sphinx .. note block would be good to alert plink users to this problem and give them a nice neat solution (we should probably focus on diploid though, for simplicity?)

See the developer docs on how to get a development environment working. To build the docs locally, go into the docs directory and type make.

[saurabhbelsare]

Hi @jeromekelleher, I've submitted a pull request including the comment in a note block in the documentation. Let me know if everything looks ok.

Thanks!

[jeromekelleher]

Thanks very much @saurabhbelsare, the note is perfect!

Closed in #541