参考文献: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为:
数据
Type | Sample1 | Sample2 |
---|---|---|
Case | NAI 1 | NAI 2 |
Control | DMSO 1 | DMSO 2 |
计算流程
设定一个覆盖度的阈值,这个阈值在BUMHMM中默认是1(是不是太低了,单是可以保留更多的位点)
计算每一个位点在四个样本中的Drop-off Rate,归一化到中位数一致;
对于Control中的两个重复,计算这两个重复之间的LDR;对于Control-Case之间的四种组合,计算每种组合的LDR;如果两个样本之间的差异很小,则LDR会接近于0,否则会远大于0;
用Control Reps之间求得的LDR作为NULL Distribution; 因此可以求得Control-Case LDR的Empirical P-value,每一个位点有四个P-value;P-Value越小表示这个位点更有可能被修饰;
建立一个隐马科夫模型:
当状态为无修饰状态时,四个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分布的参数优化出来。
根据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。