This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to merge SNP datasets aligned to different reference genome versions without re-calling variants?

I am working with SNP data from two different projects that used different versions of the reference genome.

Dataset 1:

  • 34,566,550 SNPs;
  • 1987 samples;
  • variants called using reference genome version X;
  • data available in VCF format;

Dataset 2:

  • 61,065,811 SNPs;
  • 658 samples;
  • variants called using an older version of the same reference genome (version Y);
  • data available in VCF format;

My goal is to combine the SNP information from both datasets so that all variants are represented in the same coordinate system, without performing variant calling again for the 658 samples (which is computationally unfeasible for me).

I am aware that one possible solution would be to re-call variants for the 658 samples using the same reference genome version as dataset 1, but I would like to avoid this step.

I have read that it might be possible to use liftover tools to convert variant coordinates between genome assemblies, but I am unsure whether liftover is appropriate for full VCF datasets or only for genomic coordinates.

Would it be appropriate to:

  1. liftover the VCF from the older genome assembly to the newer one
  2. then merge the VCF files (e.g. using bcftools merge or similar tools)?

Are there recommended best practices or potential pitfalls when combining SNP datasets generated using different genome assemblies?

Any suggestions or references would be greatly appreciated.

Ps: The samples from the two datasets are not the same sample, but they are of the same species

genome-assembly vcf snp variant-calling liftover

ok do you have SNPs aligned to hg19/grch37? why are you being secretive about the freeze?

Hello, Jeremy. I apologize if it seemed that I was omitting details. In fact, I was trying to be as general as possible so that everyone could understand my problem. However, as I mentioned in my reply to Giulio, after rereading my message I realized that I did not describe my problem appropriately.

I understand that my initial question may have suggested that I intended to harmonize two complete genotype callsets generated using different versions of the reference genome. However, my situation is somewhat different and more limited in scope.

In practice, one of the datasets (Dataset 2) has already been analyzed to identify a small subset of variants (on the order of dozens of variants) located in specific genomic regions of interest. These variants were identified using an older version of the reference genome (canFam3).

Independently, I am working with another genotype callset (Dataset 1), generated using a newer version of the reference genome (canFam4), and I would like to evaluate the same genomic regions in this newer dataset.

Therefore, my main objective is not to merge the two complete callsets, but rather to determine whether this small set of candidate variants (identified in canFam3 coordinates) can be mapped to canFam4 coordinates, so that genotype counts can also be evaluated in the second dataset.

Once again, I apologize if it seemed that I was omitting details, that was not my intention.

1 answer

Use BCFtools/liftover available on github. After you do that, you should keep for all downstream analyses a covariate tracking which samples had genotypes converted with BCFtools/liftover as there will be a lot of different types of biases that make the two genotype callsets not comparable

Hello Giulio, thank you very much for your suggestion regarding BCFtools/liftover; the paper was very helpful (I saw that you are the author, that’s amazing!). Thank you very much as well for your reply; I truly don’t have words to express my gratitude. However, after rereading my message, I realized that I did not explain my problem correctly.

I realize that my initial question may have suggested that I intended to harmonize two complete genotype callsets generated using different genome assemblies. However, my situation is slightly different and more limited in scope.

In practice, one of the datasets (Dataset 2) has already been analyzed to identify a small subset of variants (on the order of dozens of variants) located in specific genomic regions of interest. These variants were identified using an older version of the reference genome (canFam3).

Independently, I am working with another genotype callset (Dataset 1), generated using a newer version of the reference genome (canFam4), and I would like to evaluate the same genomic regions in this newer dataset.

Therefore, my main objective is not to merge the two complete callsets, but rather to determine whether this small set of candidate variants (identified in canFam3 coordinates) can be mapped to canFam4 coordinates, so that genotype counts can also be evaluated in the second dataset.

Given this more limited scope, I would appreciate your opinion on the following points:

1 - Do you think BCFtools/liftover would still be preferable over standard coordinate conversion tools (such as CrossMap or UCSC liftOver) when converting a relatively small set of variants?

2 - Would you recommend performing liftover only for this subset of candidate variants, or would there be methodological advantages to lifting over the full callset even if most variants are not directly relevant to the downstream analysis?

3 - In your experience, does comparing genotype counts for a small number of loci originally called on different assemblies introduce substantial bias, even if the variants are successfully lifted and validated?

4 - Are there particular types of bias or artifacts that you would recommend checking for when evaluating variants that were originally called on different genome assemblies?

5 - Would you recommend any specific QC steps after liftover to ensure that the converted variants are reliable (for example allele consistency, strand issues, or ambiguous mappings)?

I know that these are many questions and that it may seem like I am taking advantage of your goodwill, but the truth is that I have never performed liftOver of variants before, I have several doubts, and I am not sure where to clarify them. I hope it does not seem like I am abusing your help. Any guidance or best practices you could share would be greatly appreciated.

In any case, thank you very much for your assistance.

UCSC liftover is a tool to convert intervals rather than variants, so I would not suggest to use this at all. CrossMap will work mostly fine for SNPs but it cannot be trusted with indels. The main issue you can face with liftover tools is that some variants that you might have called if you had realigned your data might not show up at all after liftover conversion. Other than that the output of BCFtools/liftover is not going to be incorrect or in need of specific types of QC

Thank you very much for your response and for your help. I will try to use BCFtools/liftover. Thanks again!

Log in to answer this question.