文献阅读 | Distribution of Parental Genome Blocksin Recombinant Inbred Lines

Martin O C. Distribution of Parental Genome Blocks in Recombinant Inbred Lines[J]. Genetics, 2011, 189(2): 645-654.

整体上，作者通过理论计算，得出了RIL群体中，单倍型块的长度分布、数量分布等统计学特征，并经过与实际实验的验证证明了理论结果的准确性。在计算上，作者通过马尔可夫过程来简化计算过程，并给出了相关代码文件。作者证明了块的长度分布类似指数分布，证明了用指数分布近似RIL情况是合理的，并证明了先前一些研究采用指数分布来近似是合理的。

1.概述部分

作者首先总结了前人的研究，介绍了单倍型在当时基于SNP的高密度基因组图谱逐渐出现的时代背景下，成为了群体基因组学研究的中心，在多种研究中被应用的“现”状。当时，已经有了许多对于单倍型块的统计特征研究。这些研究可被分为两类

• 在第一类中，生物学问题是研究一个或多个父母的基因组是如何在连续的世代中分裂成小块的，以及后代如何共享IBD块。
  • 在本类中，普遍的框架允许样本间随机交配，每个个体后代数量随机，并存在种群增长的可能。
  • 正由于这种随机性，分支数学理论可以起到关键作用。
• 与第一类相反，第二类假设所有的系谱都是完整的，因此单倍型块的特征只与受控杂交相关。
  • 正因这一限制，可以使用马尔科夫链跟踪IBD区块的统计数据，甚至可以跟踪所有创始父母的基因组是如何在后代间共享的

第一类的理论通常费舍尔的“交界理论”为基础。第二类起源于Donelly (1983)，这类研究将IBD问题映射到系谱相关图上随机行走的问题。在重组自交系的背景下，采用的是马尔科夫框架

随后，介绍了重组自交系（RIL）可以通过自体受精（植物）（SSD）或兄弟姐妹交配（动物）（SIB）来获得，已成为动植物研究（遗传图谱，QTL检测和关联研究）的首选工具。并介绍了block长度的均值，0.5 M in SSD and 0.25 M in SIB if the chromosome genetic length is large。

2. 材料与方法：junction定义与寻找

随后，在材料与方法部分，作者首先并说明了假设。明确了研究中采用的是二倍体的RIL材料，并说明由于每个染色体对都是独立的，可以着眼于某一对染色体而非全部。定义了F1是第零代。在从一代到下一代的过程中，认为SSD产生一个后代，而SIB产生一对。认为locus一旦确定了位置，它就会一直固定下来。 在SSD中，这仅要求两个同源染色体在该基因座具有相同的等位基因。 在SIB中，它要求所有四个染色体都具有相同的等位基因。

接下来，作者定义了“junction”。首先以摩根为单位，测量了染色体中的连续位置，并将染色体最左端作为原点，即此时x=0。跟随Fisher (1949, 1954, 1959), Bennett (1953), and Donelly (1983)的工作，交叉（block的转换边界）被称为“junction”，并在染色体上用精确的点进行标识。同时，作者假设交叉未收到干扰，那么可以认为区间[x,x+dx]中“junction”的产生与否同其他区间无关。此处，认为x代表遗传位置，dx无限小，junctions随着染色体上密度1出现。图1说明了几代之后的SSD群体中，利用junctions的方法。图中，使用二进制标签0和1标记block的来源（PA或Pa）。对于每一世代，都有0和1的标签列表来使我们能确定在任意位置x的IBD含量，如图1右侧所示。（Note：junctions的编号是自左向右进行的。）图S1给出了后续步骤：首先是对交界处及其编号进行布局，然后在每个交界处引入二进制标签，最后一个重构单倍型。SIB交配的情况类似，并且在每一代中，染色体都由基础父母的染色体的镶嵌体组成。

图S1给出了后续步骤：首先是对交界处及其编号进行布局，然后在每个交界处引入二进制标签，最后一个重构单倍型。SIB交配的情况类似，并且在每一代中，染色体都由基础父母的染色体的镶嵌体组成。

到目前为止，可以为任何谱系系统制定联结的引入方式。 以前的许多研究都做到了这一点，并将IBD问题映射到谱系相关图上的连续时间随机游走问题。在本文中，作者考虑了RIL，然后Markovian框架的谱系相关图就是一个超立方体。 任何给定基因座x处的二进制标记序列指定了超立方体的唯一顶点，随着x在染色体上的移动，顶点的改变决定了超立方体的结构。作者还指出，这一数学框架的某些部分数学框架在本质上与在存在重组的情况下研究聚结（coalescent）所使用的相似。 那里的中心对象是所谓的“祖先重组图”，而研究的问题是该图是如何沿着连续染色体的位置变化的。由于这一问题很难，研究者得出的确切结果相对较少。

在SSD和SIB RIL应用Markovian框架，需要考虑在高维超立方体上所有可能的连续时间步移。 作者分两个步骤进行。 首先，通过计算机枚举该超立方体上所有可能的离散时间随机游动。 然后，通过分析技术解决原始步行的连续等待时间。 最后，使用Mathematica（Wolfram1991）对这些分析表达式进行数值处理。 文件S2中提供了这些不同任务的C和Mathematica代码。

3. 材料与方法：连续时间马尔可夫过程

创建第1到g代时，在SSD配对方案中生产2g个配子，在SIB配对方案中生产4g个配子。用Nc表示该数目，因为它也是RIL构建中产生的（新）染色体的数目（请记住，我们仅追踪一对染色体）。相较于一代代地追踪配子，将所有Nc个配子同时考虑更为有用。这可以通过将所有世代的染色体对相互堆叠，然后从左到右扫描染色体栈，以查看junctions按x递增顺序的出现情况的形式进行可视化。

最重要的是，这些Nc个配子上的junctions在各方面都是独立的：在[x，x + dx]中的一个配子上具有junction不会影响在其他位置出现junctions的可能性（无论是在同一个配子上还是在另一个配子上）。正因这种独立性，人们可以将整个配子中junctions的产生看作是“连续时间”马尔可夫过程（Feller 1950），其中x扮演时间的角色。对于间隔[x，x + dx]，出现junction的概率为Nc * dx，然后真的出现了了此事件（event），则将junction随机分配给Nc个配子之一（每个配子的概率为1 / Nc）。 然后在间隔[x + dx，x + 2dx]中重复该操作，依此类推。

因此，我们有一个马尔可夫过程，其中，事件间隔是独立的，并且分布可表示为下式（方程式1）。同时，junction在染色体中的分布是均等的。

\[\rho(\Delta x) = N_{c}\cdot exp(-N_{c}\Delta x) \]

4.材料与方法：超立方体上的离散和连续时间随机游走

作者随机且统一地在\(x=0\)处初始化二进制标签，因为隔离是无偏的。然后，使用连续时间游走的马尔可夫过程将标签从\(x=0\)拓展至\(x\)。对于任意\(x\)，我们称\(M\)为\(N_{c}\)维超立方体\(H=\{0,1\}^{N_{c}}>\)到第\(g\)代的基因型的映射。这一映射过程可以视为对\(N_{c}\)维超立方体顶点的着色。颜色与第\(g\)代中的单基因座基因型一样多：SSD为4种，SIB为16种。然后，第\(g\)代的块模式（block pattern）可以通过马尔可夫过程在超立方体上游走时，顶点颜色的演替（succession of vertex colors）来“读取”。请注意，在\(x\)处出现junction，相当于在那个时候马尔科夫游走跳越到了\(H\)上的一个随机的相邻顶点。而每个顶点上的停顿时间是指数分布的（请参见方程式1）。

对于块的统计数据，作者希望找到给定颜色模式（given pattern of successive colors）下，在\(H\)上随机游走得到该模式的概率。为此，作者将与所需模式兼容的所有可能步行路线概率相加得到。关键点是，交界处值的连续变量会影响块的长度，但不会影响连续块的模式。这使得作者可以将块统计问题分解为两部分。

• 首先对于来自在H上访问的顶点序列的离散集（junction的“拓扑”结构）； 作者在超立方体上使用离散时间随机游动的主方程来追踪这些序列。
• 第二个因素与junction-junction间隔的连续性质有关，这涉及对从方程式1得出的已知概率分布求和。

5. 材料与方法：获取块长度分布

考虑最简单的情况：染色体上的第一个区块的长度的获取。

• 如果该块是杂合的，则其在SSD中的长度分布是第一个junction的距离，并由等式1给出。实际上，第\(g\)代的位点\(x = 0\)必须是杂合的。但是在SSD中，在超立方体上随机行走的第一跳处，杂合块就会结束，并开始一个被固定下来的块。
• 如果第一个块是纯合的（固定的），则这种情况更具指导意义。为了计算第一个块的长度分布，作者首先考虑从起始顶点（\(x=0\)）开始的不同游走的所有可能性。游走将在第\(g\)代保持纯合的结构经过几跳，然后在某一跳改变这种情况。如果\(k\)是结束该块的第一跳，我们可以将具有相同\(k\)的所有离散时间步长收集在一起。因此，我们将\(P^{(1)}(k)\)定义为在第\(g\)代、状态保持与\(x=0\)相同时，执行\(k-1\)次跳跃，然后在第\(k\)次跳跃时终止该块的概率。\(P^{(1)}(k)\)是\(H\)上所有恰好在\(k-1\)跳期间保持在\(x=0\)处的状态的离散时间游走的概率之和。由于此类遍历（遍历所有游走的可能路径）的数量随k呈指数增长，因此最好通过递归而不是枚举来确定此数量。使用关联的主方程式时，正是这样做的。该方程式的每次迭代都会更新超立方体上的向量，并生成后续的\(P^{(1)}(j)\)。要获得\(P^{(1)}(k)\)必须对主方程进行k次迭代。文件S1指定了该主方程，迭代向量的初始化以及迭代向量与\(P^{(1)}(j)\)之间的关系； 还提供了用于实现这些迭代的C程序（请参见文件S2）。

给定\(P^{(1)}(k)\)概率，可以重新引入在超立方体的每个顶点上花费的连续时间，以获得第一个块的长度分布。事实上，对于SSD和SIB，对于所有导致这种情况的随机游走，\(x\)都会从0游走到\(x_{k}\)，而\(x_{k}\)的分布遵循一个重新缩放过的伽玛分布

\[\rho_{k}(x_{k}) = \frac{N_{c}^{k} N_{k}^{k-1}exp(-N_{c}x_{k})}{(k-1)!} \]

这样的\(k\)个独立的指数分布之和。

第一个块的长度（假设它是固定的）\(l_{1}\)的分布，即第一个块的长度（假设它是固定的）由下式（等式3）给出

\[\mu^{(1)}(l_{1}) = \sum_{k=1}^{\infty} P^{(1)}(k) \rho_{k}(l_{1}) \]

该结果适用于无限长染色体。 对于长度为\(L\)的有限染色体，如果\(l>L\)，将会“跳出”有限染色体。因此，如果\(l_{1}\)（按等式3分布）的值大于\(L\)，在有限染色体上该块实际长度只有\(L\)。因此，为了将公式适用于有限染色体，作者在\(l_{1}<L\)时保持分布不变，而在\(l_{1}>L\)时，将\(l_{1}\)的值设为\(L\)。

6. 结果

6.1 块平均长度

块的平均长度\(\left \langle l \right \rangle\)与块末端密度\(\eta\)相关：对于非常大的、长度为\(L\)的染色体，当\(\left \langle l \right \rangle \approx \frac{L}{n}\)块的数量\(n\)满足\(\frac{n}{L} \approx \eta\)。对于SSD RIL，已知密度\(\eta\)为2，而在SIB中为4。通过考虑较小的[x1, x2]，并要求间隔是经过重组的，可以得到这一结果：在SSD中，其发生的概率为\(R=\frac{2r}{1+2r}\)其中\(r\)是[x1, x2]中的每减数分裂的重组率。令\(x_{2}-x_{1}\)为无穷小，我们就能得到\(r \approx x_{2}-x_{1}\)(Haldane 1919)，因此\(R \approx2( x_{2}-x_{1})\)。考虑到重组，在此间隔中存在一个块末端，该末端的密度为2。

令\(\eta=2\)，可以直接得到\(\left \langle l \right \rangle = \frac{1}{2}\)。这一结果在较大的代数\(g\)和大染色体上都有效。相同的推理，得出对于SIB，根据\(R=\frac{4r}{1+6r}\)，有\(\left \langle l \right \rangle = \frac{1}{4}\)

6.2 块长度分布

6.3 第一块block比后面的长

6.4 连续的块的长度间略有相关

6.5 与真实结果的比较

作者将他们理论计算的结果，同2006年发表的一份拟南芥的SSD RIL数据进行比较，并得出了理论结果与实验结果非常吻合的结果。

6.6 块的数量分布