BUMHMM计算原理

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

github: 链接

Drop-off Rate:

对于某一个位点,覆盖度(BD)是n,终止次数(RT)是k,则定义它的Drop-off Rate为:

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

某一个位点,在样本 i 和 j 中,对应的Drop-off rate为,则定义LDR为:

数据

TypeSample1Sample2
CaseNAI 1NAI 2
ControlDMSO 1DMSO 2

计算流程

  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的分布

    目标函数

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

    上式等号右边的第二项对应我们需要优化的一项,因为Beta分布的参数含在这一项中。注意。我们分解右边第二项:

    注意,其中,。把上面的F函数对求导:

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

  6. 根据HMM模型中的公式可以得到

    Bias

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

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

    因此,LDR的方差为:

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

    因此,log-ratio和coverage之间的关系可以这样建模: (这里也有问题:每一个位点的p都不一样)

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

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

    结果

    http://23.105.203.142/BUMHMM/index.html