癌症中显著突变通路的检测算法

### 癌症中显著突变通路的检测算法 在癌症研究中，识别显著突变的通路对于理解癌症的发生机制和开发有效的治疗方法至关重要。本文将介绍一种用于识别癌症中显著突变通路的方法，包括数学模型、影响图、组合模型、增强影响模型以及统计分析等方面，并展示相关的实验结果。 #### 1. 数学模型 我们使用图 $G = (V, E)$ 来建模交互网络，其中 $V$ 中的顶点代表单个蛋白质（及其相关基因），$E$ 中的边代表蛋白质 - 蛋白质或蛋白质 - DNA 相互作用。设 $T \subseteq V$ 是在样本集 $S$（患者）中进行突变检测或分析的基因子集。$T$ 的大小因研究而异，例如，一些近期的研究对数百个基因进行了重新测序，而另一些则检测了人类基因组中几乎所有已知的蛋白质编码基因。 我们假设每个基因 $g$ 在每个样本中被分配为突变或正常两种标签之一。设 $M_i$ 表示第 $i$ 个样本中 $T$ 里发生突变的基因子集，$S_j$ 表示基因 $g_j \in T$ 发生突变的样本集，$m = \sum_{i} |M_i|$ 是所有样本中观察到的改变基因的总出现次数。 我们将通路或子网络定义为 $G$ 的连通子图。需要注意的是，这个定义与生物学中该术语的常见用法相匹配，即通路在图中可以具有任意拓扑结构，而不限于顶点的线性链。通常，我们不知道在一个样本中是否必须有多个基因发生突变才能扰乱一个通路，因此我们假设如果通路中的任何基因发生突变，则该通路在样本中发生了突变。对于子集 $T' \subseteq T$，$S(T')$ 表示 $T'$ 中至少有一个基因发生突变的样本集。 #### 2. 影响图 我们的目标是识别相对于样本中突变基因集具有显著意义的子网络。子网络的显著性来自两个方面： - 子网络中的基因发生突变的样本数量。 - 子网络中基因在整个网络拓扑结构中的相互作用。 例如，考虑两种可能的突变节点场景： - **场景一**：两个突变节点是交互网络中线性链的一部分。 - **场景二**：两个突变节点通过一个高度数节点连接。 在场景一中，连接两个突变节点的路径只有一条，因此我们对这种局部突变聚类更为惊讶；而在场景二中，两个节点由一个存在于大量可能路径中的节点连接。 中心节点（hubs）是这种现象的极端情况，会导致许多“无趣”的子网络被认为具有显著性。由于许多高度突变的癌症基因（如 TP53）在交互网络中也具有高度数，因此在癌症突变数据分析中忽略这些基因是不可取的。 为了定义节点对之间严格的影响度量，我们在交互网络上使用扩散过程。考虑以下过程：以恒定速率将流体泵入源节点 $s$，流体通过边在图中扩散，每个节点以恒定的一阶速率 $\gamma$ 损失流体。设 $f_s^v(t)$ 表示时间 $t$ 时节点 $v$ 处的流体量，$f_s(t) = [f_s^1(t), \ldots, f_s^n(t)]^T$ 是所有节点处流体量的列向量。设 $L$ 是图的拉普拉斯矩阵，$L_{\gamma} = L + \gamma I$，则该连续时间过程的动态由向量方程 $\frac{df_s(t)}{dt} = -L_{\gamma} f_s(t) + b_s u(t)$ 控制，其中 $b_s$ 是第 $s$ 个位置为 1 其余为 0 的单位向量，$u(t)$ 是单位阶跃函数。当 $t \to \infty$ 时，系统达到稳态，图上流体密度的平衡分布为 $f_s = L_{\gamma}^{-1} b_s$。 我们将 $f_s^i$ 解释为基因 $g_s$ 对基因 $g_i$ 的影响。对所有测试基因计算扩散过程，我们可以得到每对基因 $g_j, g_k \in T$ 之间的影响 $i(g_j, g_k)$。一般来说，这种影响是非对称的，即 $i(g_j, g_k) \neq i(g_k, g_j)$。我们定义影响图 $G_I = (T, E_I)$，其节点集对应于测试基因集，边 $(g_j, g_k)$ 的权重为 $w(g_j, g_k) = \min[i(g_k, g_j), i(g_j, g_k)]$。如果交互网络中有 $n$ 个节点，则计算 $G_I$ 的成本主要由求一个 $n \times n$ 矩阵的逆的复杂度决定。 #### 3. 发现显著子网络：组合模型 给定基因之间的影响度量，发现显著子网络的第一个明显方法是在影响图 $G_I$ 中识别满足以下条件的节点集： - 通过具有高影响度量的边连接。 - 对应于在大量样本中发生突变的基因。 我们固定一个阈值 $\delta$，通过移除 $G_I$ 中权重 $w(g_i, g_j) < \delta$ 的所有边以及样本数据中没有突变的基因对应的所有节点，计算 $G_I$ 的简化影响图 $G_I(\delta)$。计算问题就简化为在 $G_I(\delta)$ 中识别连通子图，使得对应的基因集在大量患者中发生改变。 我们发现的连通子图的大小由阈值 $\delta$ 控制。我们选择足够小的 $\delta$，使得在零假设（即突变随机放置在对应于测试基因的节点中）下，我们的方法不太可能找到具有相似属性的连通子图。需要注意的是，$\delta$ 的值仅取决于零假设，而不取决于观察到的样本数据。 寻找 $k$ 个基因的连通子图，使其在最多样本中发生突变，这需要解决我们定义的连通最大覆盖问题： 给定一个定义在 $m$ 个顶点集 $V$ 上的图 $G$、一个元素集 $I$、一个子集族 $P = \{P_1, \ldots, P_m\}$（其中 $P_i \in 2^I$ 与 $v_i \in V$ 相关联）以及一个值 $k$，找到 $G$ 中具有 $k$ 个节点的连通子图 $C^* = \{v_{i_1}, \ldots, v_{i_k}\}$，使得 $|\cup_{j = 1}^k P_{i_j}|$ 最大。在我们的情况下，$G = G_I(\delta)$，$V$ 是 $T$ 中至少在一个样本中发生突变的基因子集，对于每个 $g_i \in V$，相关联的集合是 $S_i$。 连通最大覆盖问题与最大覆盖问题相关，但它是 NP 难问题，即使在简单的图（如星图）上也是如此。由于该问题是 NP 难的，我们转向近似解决方案。我们构造了一个多项式时间算法，当最优解 $C^*$ 的半径为 $r$ 时，该算法给出 $O(1/r)$ 近似。 我们的算法从每个节点 $v \in V$ 开始获得一个解 $C_v$（即一个连通子图），然后返回找到的最佳解。为了获得 $C_v$，算法执行以下步骤： 1. **探索阶段**：对于图 $G$ 中的每个节点 $u$，找到从 $v$ 到 $u$ 的最短路径 $p_v(u)$。设 $\ell_v(u)$ 是 $p_v(u)$ 中的节点集，$P_v(u)$ 是它们覆盖的 $I$ 中的元素。 2. **构建连通子图**：从 $v$ 开始构建连通子图 $C_v$，初始时 $C_v = \{v\}$，$P_{C_v}$ 是当前连通子图 $C_v$ 覆盖的元素集。 3. **迭代扩展**：当 $|C_v| < k$ 时，选择节点 $u \notin C_v$，使得 $\frac{|P_v(