阅读 85

【群体遗传学】 π (pi)的计算

杂合度 heterozygosity

某个位点的第 个等位基因的样本频率为 ,那么该位点所有等位基因的频率和应该是1。先考虑二倍体的双等位基因,那就是 。衡量单个多态位点变异(variation)的一个方法是计算样本杂合度(heterozygosity),公式如下:

在公式中, 代表的是样本中序列的数量。

上面这个公式是针对一个位点的,如果是正对一条序列的话,那其实就就是将整条序列的杂合度加起来即可。

其中 表示的是分离位点的数量, 表示的是第 个分离位点的杂合度。在Wright-Fisher模型(无限位点的二倍体)下, ,因此有时这个统计量也叫 。我们需要注意的是在单态位点(monomorphic site)时杂合度是0。

先看这样一个例子:

假设现在有4个样本,15个位点,但是只有6个位点是分离位点,我们先计算每个分离位点的杂合度:

根据公式可知,对分离位点1(图中的第二列序列),有两个等为位点,分别是T和C,其中T有3个,C有1个,那么对T来说,它的频率就是0.75,对C来说它的频率就是0.25。根据公式可得:

我们以此计算就能得到其他5个分离位点的杂合度分别为:0.667,0.5,0.667,0.5,0.5。

那么就能计算 值了:

但是我们通常关注的是每个位点 的均值:

我们将 的计算进行推广就能得到下面这个公式:

其中 表示的是第 条序列和第 条序列之间不同核苷酸的数量,分母表示的是 个序列之间进行比较的唯一次数(非重复比较)。现在我们将这个公式应用到上面的序列中。

现在是有4条序列,所以 . 然后以此进行比较:

第一条VS第二条:3个不同的核苷酸

第一条VS第三条:4个不同的核苷酸

第一条VS第四条:3个不同的核苷酸

第二条VS第三条:5个不同的核苷酸

第二条VS第四条:0个不同的核苷酸

第三条VS第四条:5个不同的核苷酸

所以,

需要注意的是当数据量很大的时候,使用公式 计算更快

正如前面说到的,我们在计算序列之间的差异时通常是省略indel将其变成缺失值进行处理的。当使用公式 并且将indel变成缺失值时,针对不同位点 是不同的。使用公式 的话,通常会省略gap位置。

比如这个例子:

如果用第一个公式,那么 ,但是如果用第二个公式的话, 。原因是第一个公式将indel当作缺失值进行处理,而第二个公式将indel当作gap直接省略了这些位点(哪怕是在这些位点并不是分离位点)。不同的公式给出的结果也不一样,尤其是正对平均的每个位点时。因此,在处理基因组这种大数据时,通常使用 这个公式。

我们可以把 的期望方差表示成参数为 的函数。虽然在中性进化模型下,这个参数没啥用?。

如果没有重组发生的话:

从公式可以看出,和 相关的方差很大,即使样本很大时,方差也不接近于0。

通常叫 \pi S S S$进行校正:

对类似于Wright-Fisher模型处于平衡状态且有无限突变位点的群体, 也是 的估计量。

那么综上:

将这个公式应用到这个例子上:

可以看到这个公式得到的结果和前面公式计算得到的3.33很接近。

还是和前面说的一样,遇到indel不同的处理方式得到的结果不一样:

  1. 如果将indel当作缺失值进行处理,那
  2. 如果将indel当作gap进行处理,那

将这两种不同方法得到的结果相加:

同样,我们可以用参数为 的函数来表示 的期望方差(Wright-Fisher模型,没有重组发生):

如果是自由重组的话,就只是前半部分。

还可以从这个公式推断出:

我们通常会看到关于 的两种估计值: 和 ,测序错误等会造成不同的影响,因此通常需要两个值都看,还有更多的统计参数可以使用(如Tajima's D)。

作者:研究僧小蓝哥

原文链接:https://www.jianshu.com/p/53727b4f646b

文章分类
后端
文章标签
版权声明:本站是系统测试站点,无实际运营。本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 XXXXXXo@163.com 举报,一经查实,本站将立刻删除。
相关推荐