Yuliang Zhang

Nucleotide Diversity pi 值的计算

张宇亮 / 2017-12-09


$\pi$值在群体遗传学中的定义是用来衡量用来衡量一个群体的遗传信息的差异程度的,它的计算方法也很直观,这篇博文推导一下单倍体物种$\pi$值的计算方法。

$$\pi = \sum\limits_{ij}^q x_i x_j d_{ij}$$ 定义为:一个群体中任意两个等位基因之间的平均核苷酸差异数目,其中q为总的等位基因数目,Xi为第i个等位基因的频率,dij为i和j两类等位基因差异核苷酸的数目。根据定义计算可以分为两步:

推导

首先我们假设存在3类等位基因A、B、C,他们的长度相同,假设为10bp,但在某些位点的核苷酸组成不同,其基因数目分别为$n_A=3$条、$n_B=1$条、$n_C=2$条,总数$n_{total}=6$条,如下表所示:

Allel Sequence
A ‘ATGCATGCAA’
A ‘ATGCATGCAA’
A ‘ATGCATGCAA’
B ‘ACGCATGCTA’
C GTTCATGCAA’
C GTTCATGCAA’

首先计算第一步——核苷酸的差异数。可以采取的方法:从第一个碱基开始统计一直到第十位碱基为止,最后把1-10位碱基的统计结果加起来即得到总的差异数目。以第一位碱基的统计为例,此位置可以形成六类等位基因的配对,分别是A-B、A-C、B-C、A-A、B-B、C-C。对于A-B pair,我们定义一个新的变量$d_{ABP}$取值为0、1,用来表示A、B两个等位基因在第P(postion)个位点上碱基是否相同,0表示相同,1表示不相同,则在第一位碱基处所有可能的配对形成的差异核苷酸的数目为($\frac {n_A (n_A-1)}{2} d_{AA1}$$d_{AA1}$为零,略去不写):

$$N_{diff 1}=n_A n_B d_{AB1} +n_B n_C d_{BC1} +n_A n_C d_{AC1} $$

再将1到10为的差异核苷酸数目加起来,可以证明$\sum \limits_{P=1}^{10}d_{ABP} = d_{AB}$,得到:

$$N_{all} = n_A n_B d_{AB} +n_B n_C d_{BC} +n_A n_C d_{AC} $$

$N_{pair}^{all} = \frac{n (n-1)}{2} \ \ (n=n_A+n_B+n_C)$

因此: $$\pi =\frac{N_{all}}{N_{pair}^{all}}=\frac{2 n_A n_B}{n} d_{AB} +\frac{2 n_B n_C}{n} d_{BC} +\frac{2 n_A n_C}{n} d_{AC}$$

因为原始公式进行了重复计算,所以不需要乘以2,且将n-1约等于n。当然也可以由三个等位基因推广到n个等位基因,证明方法同上。

到此我们就证明了π值的求法,目前vcftools可以通过vcf文件计算π值,也可通过MEGA软件计算π值。