Bcftools homozygous calls have less read counts than heterozygotes

Question

[oselm]

Dear all,

I am trying to call SNPS across 150 individuals of a non-model species genotyped using a WGS resequencing approach.

In short: I aligned reads from each sample against the reference using BWA and subsequently used bcftools mpileup to calculate the counts and bcftools call to call genotypes. I performed this on each sample separatedly and allowed to call the consensus (i.e. equal to reference) genotypes. I then used bcftools merge to create a unique vcf file containing all the samples and filtered for missing rates.

I now want to perform a quality filter to remove genotypes with low read counts. The problem is that I noticed that heterozygous genotypes have usually more read counts than homozygotes. For this reason, filtering for read count (DP) produces a dataset where it is rare to observe a SNP with three genotypes, which doesn't make much sense...

Is it normal that heterozygotes genotypes have more DP, comapred to homozygotes? If not, what could be the cause? If yes, how can I deal with this during filtering of the vcf?

thank you in advance

OS

[finswimmer]

Hello,

can you please post some example lines of your vcf, where the differences you mentioned can be seen?

fin swimmer

[oselm]

It's not easy to generalize as there are 150 individuals x 10 M SNPs. Anyway, here I show three examples of SNPs predominantly homozygous-consensus, heterozygous and homozygous-alternative.

Homozygous-consensus:

CHR1    2725    .   C   T   225.009 PASS    MQSB=0.924584;MQ0F=0;AF1=0;AC1=0;MQ=60;FQ=-53.9861;SGB=-0.379885;RPB=1;MQB=1;BQB=1;PV4=1,1,1,1;VDB=0.06;DP=1030;DP4=338,417,37,121  GT:PL:DP:AD 0/0:0,.,.:9:9,. 0/0:23,.,.:6:5,.    0/0:24,.,.:4:3,.    0/0:0,.,.:6:5,. 0/1:48,0,167:9:7,2  0/0:0,.,.:5:5,. 0/1:33,0,28:2:1,1   0/0:0,.,.:6:6,. 0/0:0,.,.:0:0,. 1/1:41,3,0:1:0,1    0/0:0,.,.:5:5,. 0/0:0,.,.:9:9,. 0/0:20,.,.:7:6,.    1/1:69,6,0:2:0,2    0/0:0,.,.:7:7,. 0/0:0,.,.:7:7,. 0/0:0,.,.:1:1,. 0/0:18,.,.:1:0,.    0/0:0,.,.:11:11,.   0/0:29,.,.:2:1,.    0/0:0,.,.:10:10,.   0/0:0,.,.:7:7,. 0/0:0,.,.:2:2,. ...

Heterozygous:

CHR1    1282532 .   C   T   225.009 PASS    VDB=0.0609607;SGB=-0.693147;RPB=0.922626;MQB=0.0382567;MQSB=0.715533;BQB=0.385803;MQ0F=0.0684932;AF1=1;AC1=2;MQ=37;FQ=-37.9865;PV4=0.0111319,1,0.0117738,1;DP=10663;DP4=690,1958,407,6624   GT:PL:DP:AD 1/1:164,11,0:70:7,63    0/1:162,0,216:67:22,45  0/1:149,0,255:60:34,26  0/1:193,0,255:69:29,40  0/1:158,0,75:39:10,29   0/1:255,0,255:90:30,60  0/1:221,0,154:67:20,47  0/1:184,0,236:71:25,46  0/1:240,0,145:50:16,34  0/1:183,0,255:77:32,45  0/1:231,0,83:89:19,70   0/1:163,0,164:103:27,76 0/1:188,0,94:49:9,40    0/1:185,0,226:40:19,21  0/1:186,0,176:52:18,34  0/1:173,0,137:90:18,72  0/1:154,0,172:24:11,13  0/1:255,0,255:80:26,54  0/1:195,0,255:53:23,30  0/1:174,0,98:45:13,32   1/1:196,54,0:31:2,29    0/1:195,0,244:56:22,34  0/1:209,0,247:62:27,35  0/1:242,0,181:97:27,70  0/1:249,0,246:88:31,57  0/1:190,0,59:36:8,28    0/1:210,0,187:94:26,68  0/1:136,0,199:27:12,15  1/1:255,115,0:85:6,79       ...

Homozygous-alternative:

CHR2     2153   .   T   C,A 173.998 PASS    VDB=0.193383;SGB=-0.556411;MQ0F=0;AF1=1;AC1=2;MQ=26;FQ=-38.9864;RPB=1;MQB=1;BQB=1;PV4=1,5.2242e-05,1,0.422272;DP=944;DP4=5,0,902,0  GT:PL:DP:AD 1/1:77,12,0,.,.,.:4:0,4,.   1/1:150,24,0,.,.,.:8:0,8,.  1/1:147,20,0,.,.,.:10:1,9,. 1/1:193,39,0,.,.,.:13:0,13,.    1/1:142,21,0,.,.,.:7:0,7,.  ./.:.:.:.   1/1:115,21,0,.,.,.:7:0,7,.  1/1:40,6,0,.,.,.:2:0,2,.    1/1:133,18,0,.,.,.:6:0,6,.  1/1:160,36,0,.,.,.:12:0,12,.    1/1:142,18,0,.,.,.:6:0,6,.  1/1:134,27,0,.,.,.:9:0,9,.  1/1:69,9,0,.,.,.:3:0,3,.    1/1:59,6,0,.,.,.:2:0,2,.    1/1:57,6,0,.,.,.:2:0,2,.    1/1:96,12,0,.,.,.:4:0,4,.   1/1:41,6,0,.,.,.:2:0,2,.    1/1:81,5,0,.,.,.:4:1,3,.    1/1:65,12,0,.,.,.:4:0,4,.   1/1:38,6,0,.,.,.:2:0,2,.    1/1:87,22,10,72,0,69:5:0,4,1    1/1:189,45,0,.,.,.:15:0,15,.    1/1:122,18,0,.,.,.:6:0,6,.  ...