同源基因表达水平的准确估计

### 同源基因表达水平的准确估计 #### 1. 引言 在RNA测序分析中，准确估计同源基因的表达水平是一个关键问题。传统的启发式方法在处理映射到多个基因组位置的多重读取（multireads）时存在局限性。本文介绍了一种基于期望最大化（EM）算法的方法，该方法能够更准确地估计基因表达水平，并联合估计表达水平和单核苷酸多态性（SNP）变体。 #### 2. 方法 ##### 2.1 表达水平估计的EM算法 由于似然函数是凹函数，上述EM算法保证收敛到最优解。尽管它不像某些方法那样具有多项式时间保证，但在实践中，它优于HAPLOFREQ方法，并为将最大似然估计问题扩展到测序样本与参考基因组不同的表达水平和变体联合估计情况提供了基本框架。 ##### 2.2 表达水平和SNP变体的联合估计 在实际中，测序的DNA可能与参考基因组略有不同，特别是在SNP位置。为了对SNP位置进行建模，引入变量$X_k = \{X_1^k, X_2^k\}$，其中$X_1^k, X_2^k \in \{A, C, T, G\}$ ，表示测序样本在该位置的基因型。$X_k$的值是未知的，需要进行推断。 可以假设$X_k$有一个先验分布，它对应于基因组中等位基因频率的分布。该分布可以从HapMap数据、ENCODE区域以及1000基因组计划的数据中进行经验估计。 如果读取$r_i$到区域$G_j$的合理比对跨越位置$X_1, \cdots, X_l$，假设测序错误在每个位置是独立的，则可以将$p_{ij}$写为： $p_{ij} = \prod_{k} \gamma(X_k, r_i^k, k)$ 其中， $\gamma(X_k, r_i^k, k) = \begin{cases} \epsilon(k), & \text{如果 } X_1^k \neq r_i^k, X_2^k \neq r_i^k \\ 1 - \epsilon(k), & \text{如果 } X_1^k = r_i^k, X_2^k = r_i^k \\ 0.5, & \text{否则} \end{cases}$ $\epsilon(k)$是读取中位置$k$的错误率函数。错误率对位置的依赖来自技术限制，因为错误率预计会随着读取长度的增加而增加。 基于此，表达水平和SNP变体的联合估计问题可以定义如下： 定义2（MLEI - SNP）：给定一组读取$r_1, \cdots, r_m$和一组区域$G_1, \cdots, G_n$，为每个区域$G_i$找到一个概率$P_i$，为每个位置$k$找到基因型$X_k = \{X_1^k, X_2^k\} \in \{A, C, T, G\}^2$，使得$\sum_{i} P_i = 1$，并且数据的似然$L = \prod_{j = 1}^{m} \sum_{i = 1}^{n} P_i p_{ij}$最大化，其中$p_{ij} = \prod_{k = 1}^{l} \gamma(X_k, r_i^k, k)$。 带有SNP变体的EM扩展：为了最大化数据的似然，需要同时找到$P = \{P_1, P_2, \cdots, P_n\}$（使得$\sum P_i = 1$）和每个基因组位置的基因型调用$X = \{x_1, \cdots, x_k\}$，使得数据的似然$L(P, X; R) = \prod_{j = 1}^{m} \sum_{i = 1}^{n} P_i p_{ij}$最大化。 EM算法可以如下调整： - **E步骤**： $Q(P, X|P^{(t)}, X^{(t)}) = E_{M|R,P^{(t)},X^{(t)}}[\log L(P, X; R, M)]$ $= E_{M|R,P^{(t)},X^{(t)}}[\sum_{i = 1}^{m} \log P_{M(i)} p_{iM(i)}]$ $= \sum_{i = 1}^{m} \sum_{j = 1}^{n} [(\log P_j p_{ij}) \times \frac{P^{(t)}_j p_{ij}^{X^{(t)}}}{\sum_{j = 1}^{n} P^{(t)}_j p_{ij}^{X^{(t)}}}]$ - **M步骤**： $(P^{(t + 1)}, X^{(t + 1)}) = \arg \max_{P,X} Q(P, X|P^{(t)}, X^{(t)})$ $= \arg \max_{P,X} [\sum_{i = 1}^{m} \sum_{j = 1}^{n} a_{ij} \log P_j p_{ij}]$ $= \arg \max_{P,X} [\sum_{i = 1}^{m} \sum_{j = 1}^{n} a_{ij} \log P_j + \sum_{i = 1}^{m} \sum_{j = 1}^{n} a_{ij} \log p_{ij}]$ 由于上述方程中的两项是独立的，可以分别对它们进行最大化。第一项在$P^{(t + 1)}_j = \frac{\sum_{i = 1}^{m} a_{ij}}{\sum_{i = 1}^{m} \sum_{j = 1}^{n} a_{ij}}$时最大化，其中$a_{ij} = \frac{P^{(t)}_j p_{ij}^{X^{