I have ~20K single-sample VCFs, all covering the same sites, that I want to merge into a single ~20K-sample VCF.
In the past, I have used plink/plink2 to do this sort of merge, but I'm having trouble keeping REF/ALT alleles and ploidy straight[1].
I am looking for a tool that (unlike plink 1.x) does not mess at all with the values of the REF and ALT columns (even for indel variants). I'm not sure what's going on with chrX ploidy (0/0 genotypes getting converted to 0, 1/1 to 1, etc.), but I'd prefer a tool that (unlike plink2) automatically keeps the genotypes in my input VCFs.
Thank you in advance!
[1] The original version of this post included the plausibly-incorrect claim that plink/plink2 was designed for use-cases that are significantly different from the ones I am interested in. See @chrchang523's replies below for details on why my use-case might or might not be covered. Regardless, for this question, I do not want to elicit responses that aim at fixing my use of plink/plink2. My goal is to find a tool that is a better fit to my work.
use
bcftools merge.I am trying to merge ~20K single-sample VCF files. Each of these files has around 1.7M sites.
Can
bcftools mergehandle a merge of this many files quickly (i.e. in under 24h)?The way I do this sort of thing within a reasonable time (though, unfortunately, I use plink/plink2 for this) is to convert the ~20K VCF files to BED files, I then merge the ~20K BED files (which is much quicker than merging the ~20K VCF files directly), and finally I convert the merged BED to VCF.
(I use plink for the first two steps above, and plink2 for the final step. Unfortunately, the final VCF has problems, as I mentioned.)
sure, by parallelizing per chromosome, first merge per group of sqrt(20,000)=142 files and the merge the 142 files to create the final vcf per chromosome, and at the end , use
bcftools concatto concatenate each chromosome.I am intrigued by this idea, but I am unable to fill in the omitted details. Would you mind posting somewhat more explicit code?
see my answer below... How to merge 20K single-sample VCFs *without* using plink or plink2?
Please clarify that you were making the conscious choice to continue using plink 1.x as part of your old workflow, and that that was the cause of your REF/ALT issues. One of the main reasons for plink 2.0’s existence is addressing that problem. (Due to the simplicity of the .bed format, the currently-missing plink2 merge functionality can be replaced with a short shell script in this case.)
But as for 0/0 -> 0 and the like, yes, if your VCFs have fine-grained ploidy information that you want to keep, that is indeed outside of plink2’s scope. bcftools is the standard solution here.
If you knew how much I resisted and fought using two different versions of plink in the same pipeline, you may have guessed that I must not have been able to find a way around it.
Indeed, the parts of the pipeline that I have been doing with plink 1.x are ones for which I have found explicit references online that they cannot be done with plink2 at the moment. I don't have these references onhand at the moment, but I remember that it was only after finding these explicit references that I finally threw in the towel, and accepted the necessity of using both plink 1.x and plink2 in the same pipeline.
As mentioned in my previous comment, the merge step that you're currently using plink 1.x for, because it "cannot be done with plink2 at the moment", can be done with a short shell script in your case, because all your files have the same variants in the same order. The key is to export the sample-major .bed format, rather than the regular variant-major format. From there, a for-loop involving
catandtail -c +4does almost everything you need.Any other commands where plink 1.x must be used instead of plink2 should be usable in a way where your REF/ALT allele ordering is unaffected.
I will assume for now that you really do have fine-grained ploidy information in your VCF that can't be handled properly by plink2, rather than the far more common scenario of pseudoautosomal vs. non-pseudoautosomal regions of chrX which is handled properly by
--split-parduring VCF import. Again, if you do have fine-grained ploidy information, then plink2 is not the right tool for the job; you won't be able to use thecat/tail -c +4merge script that would otherwise be >100x as efficient asbcftools merge.But please edit the question to make it clear that plink2 does not flip REF/ALT alleles. It's fine if you don't want to continue using plink even if it does cover your use cases -- especially if it's due to poor documentation/support that makes it impractical for you to use the program properly! -- but a major purpose of this forum is to provide accurate information to others, so I will edit away the misleading parts of your question if you don't.
Thank you for your comments. I find them extremely helpful and informative.
I'm happy to correct the errors in my original post. Is this rewording of what I wrote in that post accurate?
The first part is fine.
As for the second part, it would be useful to clarify how your chrX use case differs from the norm, because plink2
--split-pardoes in fact handle the normal chrX use case properly and you haven't given any indication so far that you understand this; you've only said that it "mess[es] with the genotypes" without describing the pattern. (It seems likely that the actual problem is that plink2 does not print a warning/error when--vcfis used without--split-paron a chrX-containing file; this will be addressed in the next release.)I wish I understood plink 1.x and plink2 well enough to write what you'd like me to write, but unfortunately I don't. (In fact, I suspect that if I had this level of understanding, it is pretty unlikely that I would have had to post this question in the first place, because I would have been able to Google for its answer far more precisely and effectively.)
Instead, I've done the best I can do, namely to alert the readers of my post to these shortcomings, and point them to your replies in this thread.
Okay, I have edited your post to remove misleading claims about plink/plink2, while keeping the "goal is to find a tool that is a better fit to my work" bottom line.