LD pruning removes all variants

[blav]

I am using plink to remove variants in LD.  When I prune the study data, all the variants are removed.

PLINK v2.00aLM 64-bit Intel (2 Aug 2017)       www.cog-genomics.org/plink/2.0/

(C) 2005-2017 Shaun Purcell, Christopher Chang   GNU General Public License v3

Logging to qcdir/merge_germline.no_ac_gt_snps.log.

Options in effect:

  --bfile qcdir/merge_germline.no_ac_gt_snps

  --exclude range reference/high-LD-regions.txt

  --indep-pairwise 50 5 0.2

  --out qcdir/merge_germline.no_ac_gt_snps

Start time: Mon Aug 22 12:30:05 2022

122720 MB RAM detected; reserving 61360 MB for main workspace.

Using up to 38 threads (change this with --threads).

272 samples (0 females, 0 males, 272 ambiguous; 272 founders) loaded from

qcdir/merge_germline.no_ac_gt_snps.fam.

29765123 variants loaded from qcdir/merge_germline.no_ac_gt_snps.bim.

Note: No phenotype data present.

--exclude range: 981687 variants excluded.

272 samples (0 females, 0 males, 272 ambiguous; 272 founders) remaining after

main filters.

Calculating allele frequencies... done.

28783436 variants remaining after main filters.

--indep-pairwise (13 compute threads): 28783436/28783436 variants removed.

Variant lists written to qcdir/merge_germline.no_ac_gt_snps.prune.in and

qcdir/merge_germline.no_ac_gt_snps.prune.out .

End time: Mon Aug 22 12:30:26 2022

I know my samples and reference are both aligned to HG37.  To preprocess my samples (~300), I split multiallelic records from my VCFs, merged all patient VCFs into one, annotated the VCF, and converted the merged VCF to plink format.  The high-LD-regions.txt file is from plinkQC.  

Do you have any idea of what to try to solve this problem? 

Thanks for your help,

Brianne

[Christopher Chang]

1. Please update to a newer plink2 build.

2. If this still fails with a newer build, post the output of running --genotyping-rate on your data.

[blav]

Hi!

I redid with plink2/3.6 and it still fails.  The genotyping rate is  0.00662904, which I gather from some other posts is much too low.  This means I don't have enough mutations present in my samples to do this analysis correct?  Can you please provide some more information about this so I can look into why this may be (perhaps preprocessing error) and how to fix it.

Thanks,

Brianne

[Christopher Chang]

A common preprocessing mistake is generating a VCF which only contains genotypes that differ from the reference genome.  You must instead generate files that contain homozygous-reference genotypes as well, in order to distinguish them from truly missing genotype data.

[blav]

Hi!  Thanks for the advice that is definitely the problem here.  I have a merged VCF for 250 patients that only have mutation calls and I converted this to plink format.  Reading previous posts, it seems the following command should create a plink file with the reference FASTA alleles (plink2 --bfile data --fa ref.fa --make-bed --out data_ref).  Is this correct?  When I try this the genotyping rate is not improved.

Thanks,

Brianne

[Christopher Chang]

No, this problem generally occurs entirely upstream of plink.  If you use GATK for variant calling, see https://gatk.broadinstitute.org/hc/en-us/community/posts/360067785472-GenotypeGVCF-doesn-t-list-homozygous-positions- .  Other variant callers should have a similar flag to report homozygous-REF genotype calls.