问题的描述
RSEM 是典型的基于转录本定量的方法,它的比对需要下载参考转录本的fa序列,与基因组比对不同,转录本比对往往是一个基因对应多个转录本序列,因此相同基因的不同转录本之间有很大的overlap
如上图所示,当reads比对到这些isoform的overlap区,如何确定每个isoform上的reads数呢?换句话来说有的reads可以mapping到isoform 1上,而有一些可以mapping到isoform 2上,所以需要解决的问题是确定isoform 1和isoform 2上reads的相对比例
类比双硬币模型来理解RSEM
假设 reads mapping到 isoform 1的概率是 θ1 ,而 reads mapping到 isoform 2的概率是 θ2(针对于 isoform 1 和 isoform 2 的overlap区)
假设我们一共有5组观测值,每组10个值,并且我们所观测到的值只是是否mapping到isoform上,而并不清楚是mapping到 isoform 1 还是 isoform 2 上,并且mapping到 isoform 上了的用 Y 表示,没有mapping上的用 N 表示:
step 1
随机设定初始值,θ1=0.6,θ2=0.5
step 2
对于第一组数据,设某reads分配给 isoform 1 的概率为0.45,记为p1:
那么,该reads分配给 isoform 2 的概率为0.55,记为p2
这里的"分配"作为隐含变量,在实际的意义中可以表示start position 等一些潜在变量,分配可以这样理解,reads会分配给isoform 1(或isoform 2),分配完成后会有两种可能性,一种是mapping到该isoform上,另一种是没有mapping到该isoform上
同理,依次计算第2到第5组:
step 3
由step 2计算的第一组,reads分配给isoform 1(isoform 2)的概率为p1=0.45(p2=0.55)
那么对于第一组的 isoform 1 :
mapping:5 × 0.45 = 2.25(mapping 到 isoform 的有5条reads); unmapping: 5 × 0.45 = 2.25(没有mapping 到 isoform 的也有5条reads)
对于第一组的 isoform 2 :
mapping:5 × 0.45 = 2.25(mapping 到 isoform 的有5条reads); unmapping: 5 × 0.45 = 2.25(没有mapping 到 isoform 的也有5条reads)
并且依次计算第2-5组的:
e.g:
以第二组为例:
对于第二组的 isoform 1 :
mapping:9 × 0.45 = 7.2(mapping 到 isoform 的有9条reads); unmapping: 5 × 0.45 = 2.25(没有mapping 到 isoform 的有1条reads)
对于第二组的 isoform 2 :
mapping:9 × 0.2 = 1.8(mapping 到 isoform 的有9条reads); unmapping: 1 × 0.2 = 0.2(没有mapping 到 isoform 的有1条reads)
step 4
重新定义θ1 和 θ2
θ1 = 19.64 / (19.64+8.56) =0.696;θ2 = 11.26 / (11.26+7.94) =0.586
将重新定义的θ1 和 θ2返回step 1—step 4,直至算法收敛,即θ1 和 θ2收敛
网友评论