Convert dosage to hardcall genotype, inconsistent result

[Zixuan Yu]

Dear PLINK2 team,

I am trying to convert dosage (from TOPMED imputation server and stored in vcf.gz files) to hardcall genotype using plink2 and have observed inconsistent result between dosage and hardcall.

Here are the code I used to test:

plink2 --vcf myvcf.dose.vcf.gz dosage=DS \
--hard-call-threshold 0.499 \
--double-id \
--make-bed \
--out test

plink1.9 --bfile test \ 

--allow-no-sex \ 

--keep-allele-order \ 

--recode A \ 

--out test_callbyplink2  

Then I used one SNP to see if the dosage information is consistent with hardcall, and here is a screenshot for the inconsistent example. I also did a boxplot and find this inconsistent is common for all 0,1,2 hardball genotype. Could you help me see if there is something wrong with my usage of converting dosage to hardcall? 

thanks,

Zixuan

[Christopher Chang]

Please post a small VCF (could have just a few samples and variants) and a pair of .log files so that someone else can reproduce what you are talking about.

[Zixuan Yu]

Hi  Christopher,

Thank you so much for your reply. I am working on creating a testing data and will post it once I have it.

In the meantime, could you please kindly look at the my code to examine if this is the right way to create genotype call from dosage? To be more specific, I want to discard the [GT] field in the vcf files and create the genotype using dosage information. Thanks!

  plink2 --vcf myvcf.dose.vcf.gz dosage=DS \
--hard-call-threshold 0.499 \
--double-id \
--make-bed \
--out test

Best,

Zixuan

[Zixuan Yu]

Here are the sample data, I use dosage info are saved in test1.plink.dosage.gz,  test1.plink.fam,  test1.plink.map, and I created hardcall genotype using PLINK2. I examine several SNPs and find their hardcalled genotype are not consistent with dosage info. For example, SNP 6:147750 are all called as '2' while its dosage value ranged from 1.249-1.375. 

I am looking forward to any help or suggestion on this issue. 

thanks,

Zixuan

[Chris Chang]

Please post the .log files from the commands you executed on the test input.

--
You received this message because you are subscribed to the Google Groups "plink2-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plink2-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/5d517eeb-15ac-44d6-b90c-5ec14012485cn%40googlegroups.com.

[Zixuan Yu]

Here are the relevant code and log. 

[Chris Chang]

Look more closely at your VCF.  All the GT values for 6:147750 are "0|0", despite dosages as high as 0.751 as you mention.  Because you did not specify "dosage=DS", plink2 just looked at the GT values.

To view this discussion on the web visit https://groups.google.com/d/msgid/plink2-users/00c910af-8566-4ff1-a2ad-6c14a336c9ban%40googlegroups.com.

[Zixuan Yu]

Hi  Christopher
Thank you so much. I was testing with the real data since yesterday and After adding the flag,  it is working! May I ask you if you have preferred way to conduct GWAS, i.e. GWAS based on hardcall or GWAS based on dosage information.

Best,

Zixuan

[Christopher Chang]

Dosages should be used in GWAS when they're available; see e.g. Zheng J, Li Y, Abecasis GR, Scheet P (2011) A Comparison of Approaches to Account for Uncertainty in Analysis of Imputed Genotypes (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3143715/ ).