BUMHMM计算原理

参考文献：Robust statistical modeling inproves sensitivity of high-throughput RNA structure probing experiments

github: 链接

Drop-off Rate:

对于某一个位点，覆盖度（BD）是n，终止次数（RT）是k，则定义它的Drop-off Rate为：

$r=\frac{k}{n}$

Log-ratio of Drop-off Rates(LDR):

某一个位点，在样本 i 和 j 中，对应的Drop-off rate为$r_i$、$r_j$，则定义LDR为：

​ $log\frac{r_i}{r_j}=log(r_i)-log(r_j)$

数据

计算流程

1. 设定一个覆盖度的阈值，这个阈值在BUMHMM中默认是1（是不是太低了，单是可以保留更多的位点）

2. 计算每一个位点在四个样本中的Drop-off Rate，归一化到中位数一致；

3. 对于Control中的两个重复，计算这两个重复之间的LDR；对于Control-Case之间的四种组合，计算每种组合的LDR；如果两个样本之间的差异很小，则LDR会接近于0，否则会远大于0；

4. 用Control Reps之间求得的LDR作为NULL Distribution; 因此可以求得Control-Case LDR的Empirical P-value，每一个位点有四个P-value；P-Value越小表示这个位点更有可能被修饰；

5. 建立一个隐马科夫模型：

   

   • 转换概率：假设RNA中期望20个双链链碱基连续；5个单链碱基连续（这个概率文献没有过多的讨论）
   • emission prob：四个p-values

   当状态为无修饰状态时，四个pvalue服从0-1的均匀分布；当状态为修饰状态时，四个pvalue服从Beta(α, β)分布，其中的参数初始化为1和10，通过EM算法不断地修正这两个参数达到拟合。

   下面是一些公式：

   empirical P values的分布

   $P(v_t|h_t=0) \sim U(0,1)$

   $P(v_t|h_t=1) \sim Beta(\alpha, \beta), with\ \alpha=1; \beta=10$

   目标函数

   $<log\ p(v_{1:T},h_{1:T}|\alpha,\beta)>=<log\ p(h_1)>+<\sum_{t=1}^T{\sum_{n=1}^N{log\ p(v_t^n|h_t)}}>+<\sum_{t=1}^{T-1}{log\ p(h_{t+1}|h_t)}>$

   其中，t=1..T表示T个碱基；n=1...N表示N对treatment-control配对。其中，

   $<\sum_{t=1}^T{\sum_{n=1}^N{log\ p(v_t^n|h_t)}}>=$

   ​ $\sum_{t=1}^T{\sum_{n=1}^N{log\ p(v_t^n|h_t=0)p(h_t=0|v_{1:t}^n)}}+\sum_{t=1}^T{\sum_{n=1}^N}log\ p(v_t^n|h_t=1)p(h_t=1|v_{1:T}^n)$

   上式等号右边的第二项对应我们需要优化的一项，因为Beta分布的参数含在这一项中。注意$p(h_t=0|v_{1:t}^n)=1-p(h_t=1|v_{1:t}^n)$。我们分解右边第二项：

   $\sum_{t=1}^T{\sum_{n=1}^N}log\ p(v_t^n|h_t=1)p(h_t=1|v_{1:T}^n)=\sum_{t=1}^T{\sum_{n=1}^N{\gamma_t\ log\frac{(v_t^n)^{\alpha-1}(1-v_t^n)^{\alpha-1}}{B(\alpha, \beta)}}}=F$

   注意，其中，$\gamma_t=p(h_t=M|v_{1:T}^n)$。把上面的F函数对$\alpha$和$\beta$求导：

   $\frac{\delta F}{\delta \alpha}=\sum_{t=1}^T{\sum_{n=1}^N{\gamma_tlogv_t^n-\gamma_t(\varphi_0(\alpha)-\varphi_0(\alpha+\beta))}}=0$

   $\frac{\delta F}{\delta \beta}=\sum_{t=1}^T{\sum_{n=1}^N{\gamma_tlog(1-v_t^n)-\gamma_t(\varphi_0(\beta)-\varphi_0(\alpha+\beta))}}=0$

   可以发现当上面式子中的$\gamma_t$抵消了，这个优化过程与它无关。因此，我们可以用EM算法把beta分布的参数优化出来。

6. 根据HMM模型中的公式可以得到$p(h_t=0|v_{1:t}^n)$和$p(h_t=1|v_{1:t}^n)$。

   ​

   Bias

   这里我们主要讨论两类Bias，一类是Coverage Bias；另外一类是Sequence Bias

   我们假设Drop-off Count服从二项分布，参数为p,n。则Drop-off Count k的方差为：

   ​ $\sigma[k]=\sqrt{np(1-p)}$

   因此，LDR的方差为：

   ​ $\sigma[LDR]\propto \frac{\sigma[k]}{n}=\frac{\sqrt{p(1-p)}}{\sqrt{n}}$

   (这里的推导感觉有问题，不知道怎么推导出来的)

   因此，log-ratio和coverage之间的关系可以这样建模：$k\frac{1}{\sqrt n}+b$ (这里也有问题：每一个位点的p都不一样)

   这里的未知数是k和b，所以我们可以使用数据来拟合：

   ​ 把所有的位点按照Coverage从小到大排序，分成不同的bin，在每一个bin中，计算位点的Coverage的均值和：LDR的95%位点，把这两个值作为输入训练得到k和b。

   

   ​

   结果

   http://23.105.203.142/BUMHMM/index.html