32、基因组学中的变异检测与泛基因组分析

基因组学中的变异检测与泛基因组分析

1. 基因组学基础与变异检测概述

在基因组学研究中，参考基因组代表了一个群体的平均基因组，研究群体中常见的变异类型具有普遍意义。同时，特定个体（供体）的基因组也备受关注。例如，在癌症遗传学中，对携带遗传性疾病的家族进行测序，寻找该家族特有的变异（即种系突变）是常规操作。体细胞突变则是不影响生殖细胞的变异，有时会影响细胞代谢并扩散到附近细胞，可能与现有的种系突变共同导致癌症。

从测序数据中识别变异的过程称为变异检测，也常被称为变异检测/分析。根据所寻找变异的规模不同，采用的方法也有很大差异，大多数方法基于读取比对技术和后续的读取堆积分析。

2. 小变异检测

• 单核苷酸多态性（SNPs） ：读取堆积能直接提供供体基因组中SNP的信息。如果T中位置j被rj条读取覆盖，其中p%的读取表明T在位置j包含核苷酸a，其余表明包含核苷酸b，就可以判断这是由于杂合多态性还是测量误差。测量误差通常是独立事件，在同一位置观察到多个误差的概率会呈指数下降，因此通过足够高的测序覆盖度，可以较为准确地获得供体基因组的SNP图谱，同时分析时还可以考虑映射准确性，对可信的比对赋予更高权重。
• 短插入缺失（indels） ：虽然基本思路与SNP检测类似，但需要更多注意。首先，如果读取比对是基于k误差搜索问题进行的，只能观察到大小至多为k的indel。其次，由于读取是独立比对的，它们在indel的确切位置上可能存在分歧。为解决这个问题，可以对覆盖位置j的rj条读取进行多重比对（重新比对），将indel定位到一个共识位置。这种多重比对的规模通常是可管理的，甚至可以对每个子问题应用最优算法。

3. 大变异检测

• 读取堆积分析 ：如果供体基因组中存在比读取比对误差阈值更长的缺失，读取堆积中应该有与缺失长度相同的未覆盖区域；如果存在插入，读取堆积中应该有一对连续位置（j, j + 1），没有读取同时覆盖j和j + 1。此外，应该有读取的前缀与某些Tj′..j匹配，后缀与某些Tj+1..j′′匹配。对于插入区域，可以使用片段组装方法，考虑所有未映射到基因组任何位置的读取，以及长度为m（读取长度）的子串Tj−m+1..j和Tj+1..j+m+1，尝试创建一个以Tj−m+1..j开头、以Tj+1..j+m+1结尾的重叠群。
• 双端读取法 ：考虑双端读取对（L1, R1）, （L2, R2）, …, （Lr, Rr），假设每个Li覆盖参考基因组的同一位置j，已知输入双端读取之间的预期距离k及其方差。可以计算所有Li和Ri映射位置之间观察到的距离的平均值k′和方差，通过测试距离偏差是否具有统计学意义来确定预期的indel长度为|k - k′|。插入长度对于片段组装方法也很有用，它对要组装的重叠群施加了约束，这与间隙填充问题相同。对于比双端读取长度更长的插入，可以通过分析有一端未映射到任何位置的读取对来找到插入区域的边界，但插入大小无法预测。
• 分裂读取比对法 ：跨越供体序列中缺失区域的读取可以与参考基因组比对，使其前缀映射到缺失区域之前，后缀映射到缺失区域之后。当许多这样的分裂读取比对都指向相同的边界时，就可以准确检测到缺失区域。

4. 双端读取法与图优化问题

在二倍体生物中，indel可以是纯合或杂合的，杂合indel还可能重叠。因此，需要将每个双端读取分配给恰好支持一个变异假设，并选择最少的假设集来形成预测。

假设我们有一组通过分裂读取映射、读取堆积下降或双端比对找到的合理变异假设H，每个假设是一个对（j, g），其中j是indel的起始位置，g是其预期长度。如果双端读取r的左部分比对到j的左侧，右部分比对到j的右侧，且比对长度与预期indel大小的差异g′满足|g - g′| ≤ t（t为给定阈值），则称r的比对支持假设（j, g）。这样就得到一个加权二分图G = (R ∪ H, E)，其中r ∈ R是双端读取，h = (j, g) ∈ H是假设，(r, h) ∈ E表示r的比对支持h，边（r, h）的成本为c(r, h) = |g - g′|。

由此引出了以下问题：
- 问题14.1：最小成本最小覆盖分配 ：给定二分图G = (R ∪ H, E)和成本函数c : E → Q+，找到G的边的子集M ⊆ E，使得每个r ∈ R恰好与M中的一条边e相关联，并且H中与M中的边相关联的顶点数量最少。在满足这两个条件的所有集合M中，找到使总成本∑e∈M c(e)最小的集合。
- 定理14.1 ：问题14.1是NP难的。证明通过将其归约为NP难的顶点覆盖问题。给定无向图G = (V, E)，创建二分关联图G∗ = (E ∪ V, A)，将所有边的成本设为零。可以证明G有大小至多为k的顶点覆盖当且仅当问题14.1在G∗上有满足dM(V) ≤ k的解。
- 问题14.2：最小成本集合覆盖 ：如果放宽问题14.1中每个读取应恰好支持一个假设的要求，就得到了经典的最小成本集合覆盖问题。给定集合{r1, …, rm}的子集集合H = {H1, …, Hn}和成本函数c : {H1, …, Hn} → Q+，找到子集C ⊆ H，使得每个ri至少属于C中的一个Hj，并使∑Hj∈C c(Hj)最小。对于变异检测问题，可以将每个假设的成本设为支持它的读取的成本之和，该问题存在一个多项式时间算法，其产生的解的总成本至多是最优解的O(log|H|)倍。

以下是大变异检测方法的总结表格：
| 方法 | 原理 | 优点 | 缺点 |
| — | — | — | — |
| 读取堆积分析 | 通过观察读取堆积中的未覆盖区域和不连续位置检测变异 | 直观，可初步判断变异情况 | 对背景噪声敏感 |
| 双端读取法 | 计算双端读取之间的距离偏差检测变异 | 可估计indel长度 | 对于长插入的插入大小难以预测 |
| 分裂读取比对法 | 利用跨越缺失区域的读取检测缺失 | 能准确检测缺失边界 | 对读取质量要求较高 |

下面是大变异检测的mermaid流程图：

graph TD
    A[开始] --> B[读取堆积分析]
    B --> C{是否有未覆盖区域或不连续位置}
    C -- 是 --> D[判断缺失或插入]
    C -- 否 --> E[双端读取法]
    E --> F[计算距离偏差]
    F --> G{偏差是否显著}
    G -- 是 --> H[估计indel长度]
    G -- 否 --> I[分裂读取比对法]
    I --> J[检测缺失边界]
    D --> K[结束]
    H --> K
    J --> K

5. 泛基因组上的变异检测

在之前讨论了考虑群体中已发现变异的读取比对方法后，接下来探讨基于这些比对预测供体序列的两种方法，其核心思路是先预测供体序列中包含的群体已知变异，再预测供体特有的变异，供体特有的变异预测采用前面提到的方法，以预测的供体序列为参考。

5.1 基于个体基因组集合的比对

假设我们有d个参考序列的读取比对，并且已知它们的多重比对，可创建一个d × n矩阵C，其中ci,j表示在多重比对的第j列中，第i个基因组被读取覆盖的次数。

• 简单预测策略 ：一种简单的预测供体序列的方法是在每一列j中选择ci,j值最大的行i，从多重比对中提取相应的核苷酸并去除间隙符号，即可得到预测的供体基因组。但这种方法会任意重组底层序列，为了更贴近生物过程，可以通过给定参数限制重组次数，由此引出以下问题。
• 问题14.3：带重组约束的最重路径 ：给定一个d × n矩阵C（ci,j ≥ 0）和阈值K < n，找到C中最多遍历K行的路径P，即一个细胞序列P = (i1, 1), (i1, 2), …, (i1, j1), (i2, j1 + 1), (i2, j1 + 2), …, (i2, j2), …, (iK, jK−1 + 1), (iK, jK−1 + 2), …, (iK, n)，使路径上所有细胞的ci,j之和最大。

该问题可以通过动态规划解决，设v(i, j, k)表示子矩阵C1..i,1..j中所有重组k - 1次的路径的最优解值，其递推公式为：
[v(i, j, k) = \max{v(i, j - 1, k) + c_{i,j}, \max{v(i’, j - 1, k - 1) + c_{i,j} | 1 \leq i’ \leq d \text{ 且 } i’ \neq i}}]
初始条件为v(i, 1, 1) = ci,1（1 ≤ i ≤ d），整个矩阵C的最优解值为max{v(i, n, K) | 1 ≤ i ≤ d}。在合适的计算顺序下，矩阵v可以在O(dnK)时间内计算完成，通过回溯可以提取最优路径，采用空间高效的实现方法，回溯可以在O(dnK)时间和O(d√nK)空间内完成。

以下是该算法的复杂度分析表格：
| 操作 | 时间复杂度 | 空间复杂度 |
| — | — | — |
| 计算矩阵v | O(dnK) | O(dnK) |
| 回溯提取路径（简单实现） | O(dnK) | O(dnK) |
| 回溯提取路径（空间高效实现） | O(dnK) | O(d√nK) |

下面是带重组约束的最重路径问题求解的mermaid流程图：

graph TD
    A[开始] --> B[初始化矩阵v]
    B --> C[按递推公式计算v]
    C --> D[找到最大v(i, n, K)]
    D --> E[回溯提取最优路径]
    E --> F[结束]

5.2 基于群体标记有向无环图（DAG）的比对

回顾泛基因组索引的方法，参考基因组和已知变异可以用标记DAG G表示。在提取每个读取比对的路径时，可以为每个被读取覆盖的顶点v增加计数器c(v)。

由于此时无法获取每个个体基因组的信息，假设供体序列是G中从源点到汇点的最重路径，即顶点权重之和最大的路径。之前为有向无环图上的最短路径问题（问题4.1）推导的动态规划递推公式同样适用于此问题，算法的运行时间与G的弧数成线性关系。也可以通过将G中的每个顶点v拆分为vin和vout，创建一个新的有向无环图G∗，使所有v的入邻接点成为vin的入邻接点，所有v的出邻接点成为vout的出邻接点，并添加权重为 -c(v)的弧(vin, vout)，G∗中从源点到汇点的最短路径就是G中的最重路径。

以下是基于标记DAG的变异检测步骤列表：
1. 对读取进行比对，提取比对路径。
2. 为被读取覆盖的顶点v增加计数器c(v)。
3. 创建有向无环图G∗。
4. 在G∗中找到从源点到汇点的最短路径，即供体序列的最重路径。

下面是基于标记DAG的变异检测的mermaid流程图：

graph TD
    A[开始] --> B[读取比对并提取路径]
    B --> C[更新顶点计数器c(v)]
    C --> D[创建有向无环图G*]
    D --> E[在G*中找最短路径]
    E --> F[得到供体序列]
    F --> G[结束]