background

最近我导兴致大起准备进行组间交流开展群体遗传学的读书交流会，为了贪图内容的简易我先选了遗传漂变的部分

但并不是说要水过，多少要有点b格。

在遗传漂变的部分里有一个重要的模型Wright-Fisher 模型，很多群体遗传学教科书都会提到这个模型，但是关于该模型中的一个细节

假设存在一个单倍体或双倍体群体，当种群中不存在选择作用，种群迁入迁出，突变，每个世代的种群成员不会和上一个世代的成员重叠，该世代的成员完全由上一个世代的成员的克隆随机抽取而来，并且种群大小是有限且固定大小的时候，又假设种群中只有两种基因型经过无限长的世代更迭后，种群中两个基因型必然有一个基因型被淘汰另一个基因型在种群被固定.

我模糊的记得这个命题的一种证明方式是通过双倍体群体杂合度随世代递减进行证明的，但是在这里我又想了一个通过马尔可夫链进行证明的一个方法，这样遗传漂变导致的基因固定也可以在单倍体群体上进行阐述。

Markov chain

我们假设我们手上有一个种群大小为N的单倍体细胞群体，经过一个世代后，种群的成员由随机选取上一个世代的所有成员的克隆进行更迭.

我们假设这个种群中只有两种基因型A，a，没有突变，没有种群的迁入迁出，没有选择作用

这样种群中不同基因型数量的组合一共有N种,这N种基因型数量组合构成了种群数量变化的事件空间，for example：

\[\begin{array}{l} A:N & a:0,\\ A:N-1 & a:1,\\ A:N-2 & a:2, \\ A:N-3 & a:3,\\ A:N-4 & a:4, \\ A:N-5 &a:5,\\ ... A:0 & a:N,\\ \end{array} \]

通过上图，对于只有两个基因型的单倍体种群我们其实只用观察其中一种基因型的数量，就可以很好的描述种群此时所处的状态了

我们约定观察基因型A，当A在种群中的数量为0的时候，记为基因型A丢失，当A在种群中的数量为N的时候，记为基因型A被固定我们在这里统称基因A丢失和基因A被固定为种群的稳定态(当基因型A丢失的时候种群的基因型也就被基因型a完全接管了)，当种群中基因型A和基因型a都不为0的时候我们称此时种群的数量状态为过渡态。

如果基因型A进入到稳定态，由于不存在突变的作用，在之后的世代里种群的基因型数量不会再变化，因此种群从稳定态转移到过渡态的概率设为0

又因为假设中提到我们是随机选取上一个世代的所有成员的克隆进行更迭，这是一个可放回式的抽样，假设种群的大小为N，其中基因型为A的单倍体个体数量为i,那么我们抽取到基因型A数量为j的概率就为

\[ P(j|K) = C_{N}^{j}(\frac{K}{N})^{j}(\frac{N-K}{N})^{N-j} \]

通过上式我们可以构建一个概率转移矩阵\(P\)

为了直观起见我们以N=6为例构建一个概率转移矩阵

\[\begin{array}{l} \left. C_{6}^{0} \frac{0}{6}^{0} \frac{6}{6}^{6} \right. & \left. C_{6}^{1} \frac{0}{6}^{1} \frac{6}{6}^{5} \right. & \left. C_{6}^{2} \frac{0}{6}^{2} \frac{6}{6}^{4} \right. & \left. C_{6}^{3} \frac{0}{6}^{3} \frac{6}{6}^{3} \right. & \left. C_{6}^{4} \frac{0}{6}^{4} \frac{6}{6}^{2} \right. & \left. C_{6}^{5} \frac{0}{6}^{5} \frac{6}{6}^{1} \right. & \left. C_{6}^{6} \frac{0}{6}^{6} \frac{6}{6}^{0} \right.\\ \left. C_{6}^{0} \frac{1}{6}^{0} \frac{5}{6}^{6} \right. & \left. C_{6}^{1} \frac{1}{6}^{1} \frac{5}{6}^{5} \right. & \left. C_{6}^{2} \frac{1}{6}^{2} \frac{5}{6}^{4} \right. & \left. C_{6}^{3} \frac{1}{6}^{3} \frac{5}{6}^{3} \right. & \left. C_{6}^{4} \frac{1}{6}^{4} \frac{5}{6}^{2} \right. & \left. C_{6}^{5} \frac{1}{6}^{5} \frac{5}{6}^{1} \right. & \left. C_{6}^{6} \frac{1}{6}^{6} \frac{5}{6}^{0} \right.\\ \left. C_{6}^{0} \frac{2}{6}^{0} \frac{4}{6}^{6} \right. & \left. C_{6}^{1} \frac{2}{6}^{1} \frac{4}{6}^{5} \right. & \left. C_{6}^{2} \frac{2}{6}^{2} \frac{4}{6}^{4} \right. & \left. C_{6}^{3} \frac{2}{6}^{3} \frac{4}{6}^{3} \right. & \left. C_{6}^{4} \frac{2}{6}^{4} \frac{4}{6}^{2} \right. & \left. C_{6}^{5} \frac{2}{6}^{5} \frac{4}{6}^{1} \right. & \left. C_{6}^{6} \frac{2}{6}^{6} \frac{4}{6}^{0} \right.\\ \left. C_{6}^{0} \frac{3}{6}^{0} \frac{3}{6}^{6} \right. & \left. C_{6}^{1} \frac{3}{6}^{1} \frac{3}{6}^{5} \right. & \left. C_{6}^{2} \frac{3}{6}^{2} \frac{3}{6}^{4} \right. & \left. C_{6}^{3} \frac{3}{6}^{3} \frac{3}{6}^{3} \right. & \left. C_{6}^{4} \frac{3}{6}^{4} \frac{3}{6}^{2} \right. & \left. C_{6}^{5} \frac{3}{6}^{5} \frac{3}{6}^{1} \right. & \left. C_{6}^{6} \frac{3}{6}^{6} \frac{3}{6}^{0} \right.\\ \left. C_{6}^{0} \frac{4}{6}^{0} \frac{2}{6}^{6} \right. & \left. C_{6}^{1} \frac{4}{6}^{1} \frac{2}{6}^{5} \right. & \left. C_{6}^{2} \frac{4}{6}^{2} \frac{2}{6}^{4} \right. & \left. C_{6}^{3} \frac{4}{6}^{3} \frac{2}{6}^{3} \right. & \left. C_{6}^{4} \frac{4}{6}^{4} \frac{2}{6}^{2} \right. & \left. C_{6}^{5} \frac{4}{6}^{5} \frac{2}{6}^{1} \right. & \left. C_{6}^{6} \frac{4}{6}^{6} \frac{2}{6}^{0} \right.\\ \left. C_{6}^{0} \frac{5}{6}^{0} \frac{1}{6}^{6} \right. & \left. C_{6}^{1} \frac{5}{6}^{1} \frac{1}{6}^{5} \right. & \left. C_{6}^{2} \frac{5}{6}^{2} \frac{1}{6}^{4} \right. & \left. C_{6}^{3} \frac{5}{6}^{3} \frac{1}{6}^{3} \right. & \left. C_{6}^{4} \frac{5}{6}^{4} \frac{1}{6}^{2} \right. & \left. C_{6}^{5} \frac{5}{6}^{5} \frac{1}{6}^{1} \right. & \left. C_{6}^{6} \frac{5}{6}^{6} \frac{1}{6}^{0} \right.\\ \left. C_{6}^{0} \frac{6}{6}^{0} \frac{0}{6}^{6} \right. & \left. C_{6}^{1} \frac{6}{6}^{1} \frac{0}{6}^{5} \right. & \left. C_{6}^{2} \frac{6}{6}^{2} \frac{0}{6}^{4} \right. & \left. C_{6}^{3} \frac{6}{6}^{3} \frac{0}{6}^{3} \right. & \left. C_{6}^{4} \frac{6}{6}^{4} \frac{0}{6}^{2} \right. & \left. C_{6}^{5} \frac{6}{6}^{5} \frac{0}{6}^{1} \right. & \left. C_{6}^{6} \frac{6}{6}^{6} \frac{0}{6}^{0} \right.\\ \end{array} \]

我们可以给出具体的数值

\[\begin{array}{l} 1.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000\\ 0.33490 & 0.40188 & 0.20094 & 0.05358 & 0.00804 & 0.00064 & 0.00002\\ 0.08779 & 0.26337 & 0.32922 & 0.21948 & 0.08230 & 0.01646 & 0.00137\\ 0.01562 & 0.09375 & 0.23438 & 0.31250 & 0.23438 & 0.09375 & 0.01562\\ 0.00137 & 0.01646 & 0.08230 & 0.21948 & 0.32922 & 0.26337 & 0.08779\\ 0.00002 & 0.00064 & 0.00804 & 0.05358 & 0.20094 & 0.40188 & 0.33490\\ 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 1.00000\\ \end{array} \]

当我给定一个当前的种群数量状态或者一个种群数量状态的概率向量\(\pi\)(行向量)

我们可以通过概率向量\(\pi\)右乘概率矩阵，得到经过一个世代的繁衍后种群所处的种群数量状态

我们求概率转移矩阵100次幂，也就是给定经过100个世代的概率转移矩阵

以上图的矩阵为例

\[\begin{array}{l} 1.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000\\ 0.83333 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.16667\\ 0.66667 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.33333\\ 0.50000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.50000\\ 0.33333 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.66667\\ 0.16667 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.83333\\ 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 0.00000 & 1.00000\\ \end{array} \]

可以发现经过100个世代过后，给定任意一个初始种群数量状态，种群基本没有可能处于过渡态

要么固定要么丢失，所以我们该如何严谨的通过这个概率转移矩阵去证明固定有限大小的种群，当时间趋近于无穷的时候种群的遗传组成必然处于稳定态呢？

我们首先对矩阵进行一下初等变换，先对基因型A固定的条件概率所在的列和矩阵的第二列进行交换，然后将当种群数量状态为N时概率所在的行与矩阵的第二行进行交换，这种变换并不存在数值上的变化仅仅是更方便我们进行表示，它仍然是一个概率转移矩阵，前两行表示当种群数量状态为稳定态时的转移概率，后N-2行表示种群数量状态为过渡态时的转移概率

此时矩阵的分块表示为

\[\left(\begin{array}{l} I & O\\ A' & A\\ \end{array}\right) \]

\(I\)是一个\(2\times2\)的单位矩阵

\(P\)是一个\((N-2)\times(N-2)\)的矩阵，根据概率转移矩阵的设定我们能知道这个\(A\)矩阵的元素是大于0小于1的，同时行上元素之和大于0小于1

当我们求矩阵的N次幂的时候

矩阵的分块表示为

\[\left(\begin{array}{l} I & O\\ A' & A\\ \end{array}\right)^{n} = \left(\begin{array}{l} I & O\\ X & A^{n}\\ \end{array}\right) \]

X太过冗长这里不做表示，我们想要证明当世代趋近于无穷的时候，种群基因型A要么固定要么丢失，只需要证明给定任意初始状态\(A^{n}\)元素全为0，也就是无穷世代后转移到任意过渡态的概率为0就可以了

我们其实不太好通过之前给出的概率转移矩阵的具体解析形式直接推导

我想到的一个方法就是我们从矩阵\(A\)的性质入手

首先，我们能知道这个\(A\)矩阵的元素是大于0小于1的，同时行上元素之和大于0小于1

于是我们来证明一条定理，也就是对于两个矩阵A,B，A和B的元素都是大于等于0小于1的，同时行上元素之和大于等于0小于1，它们的矩阵乘\(C = A\times B\)的行上元素之和大于等于0小于1，而且C得行元素之和要小于矩阵A的对应行元素之和

为直观起见这里给出A，B的矩阵表示

\[\left(\begin{array}{l} a_{11} & a_{12} & ... & a_{1n}\\ a_{21} & a_{22} & ... & a_{2n}\\ \vdots & \vdots & & \vdots\\ a_{n1} & a_{n2} & ... & a_{nn} \end{array}\right) \left(\begin{array}{l} b_{11} & b_{12} & ... & b_{1n}\\ b_{21} & b_{22} & ... & b_{2n}\\ \vdots & \vdots & & \vdots\\ b_{n1} & b_{n2} & ... & b_{nn} \end{array}\right) \]

对于矩阵C第j行的元素我们有

\[\begin{aligned}{} & c_{j1} = a_{j1}b_{11}+a_{j2}b_{21}+a_{j3}b_{31}+...a_{jn}b_{n1}\\ & c_{j2} = a_{j1}b_{12}+a_{j2}b_{22}+a_{j3}b_{32}+...a_{jn}b_{n2}\\ & c_{j3} = a_{j1}b_{13}+a_{j2}b_{23}+a_{j3}b_{33}+...a_{jn}b_{n3}\\ & ...\\ & c_{jn} = a_{j1}b_{1n}+a_{j2}b_{2n}+a_{j3}b_{3n}+...a_{jn}b_{nn}\\ \end{aligned} \]

它的第j行的元素之和为

\[\begin{aligned}{} a_{j1}(b_{11}+b_{12}+b_{13}+...+b_{1n})+\\ a_{j2}(b_{11}+b_{12}+b_{13}+...+b_{1n})+\\ a_{j3}(b_{11}+b_{12}+b_{13}+...+b_{1n})+\\ ...\\ a_{jn}(b_{11}+b_{12}+b_{13}+...+b_{1n})\\ \end{aligned} \]

我们之前提到矩阵B的行向量之和是大于0小于1的，矩阵C的行向量之和相当于矩阵A对应行向量元素，每个元素乘上一个小于1大于0的数再作和，所以矩阵C的行向量元素之和一定是小于A的对应行向量元素之和大于0的.

现在回到群体遗传学wright-fisher model的马尔可夫链的那个问题里，在此问题中的矩阵A的n次幂的任意行向量之和所构成的序列根据上述定理是一个递减序列，同时我们又知道这个和肯定大于0所以根据我们高数中所学的定理，单调有界序列收敛定理，当n趋近于无穷的时候这个行向量元素之和趋近于一个固定值的，另一方面我们可以进一步通过极限的定义\(\delta-\sigma\)语言结合刚才证明的定理，进一步说明这个值是0.

最后我们知道行元素它是大于0的，另一方面行元素肯定是小于行元素之和的，所以通过夹逼定理我们可以知道，当世代数n趋近于无穷的时候矩阵A的n次幂，里的元素都是趋近于0的，也就是当给定任意初始状态的时候经过无穷世代，种群几乎不可能为过渡态，证明完毕