原理

其中，Xg,1—Xg,5代表 gene g 在时间点1—5的平均表达量（生物学重复的平均表达量），作为HMM的隐含层；而SgΔ，相对于上一个时间点基因的变化量作为HMM的观测层

1. 当 Xg,t-1 < Xg,t 时，定义SgΔ 为 UP；
2. 当 Xg,t-1 > Xg,t 时，定义SgΔ 为 DOWN；
3. 当 Xg,t-1 = Xg,t 时，定义SgΔ 为 EE；

那么在分析中我们可以观察到的状态是SgΔ ，而隐含层序列 Xg,t 和HMM的参数是未知的，所以该问题是已知观测序列，学习模型的参数，利用Baum-Welch算法求解

HMM的参数为：λ = ( Π，A，B)，其中Π是初始状态概率，A为转移概率矩阵，B为发射概率矩阵

假设观测序列为 O = (o1, o2, o3.....)；隐含序列为 I = (i1, i2, i3.....)；那么我们的目的是寻找一组合适的 λ 使得 P( O | λ ) 最大

如下图所示：

其中 O 为观测序列，即 "up，down，EE"，对于每个 gene，观测序列是已知的；而隐含层 I （基因的表达量）是未知的；HMM的参数 λ 是未知的
算法的目的是确定合适的 λ 使得 P( O | λ ) 最大；而P( O | λ )的统计学意义是在λ条件下，观测到 O 的概率

Baum-Welch算法的推导可以参见李航的《统计学习方法第二版》P203-207，该算法本质上是EM算法的一个变种，目的是学习HMM的参数，即确定一个参数使得似然函数值 P( O | λ ) 最大

代码分析

library(EBSeqHMM)

data(GeneExampleData)
str(GeneExampleData)

CondVector=rep(paste("t",1:5,sep=""),each=3)
print(CondVector)
Conditions=factor(CondVector, levels=c("t1","t2","t3","t4","t5"))
str(Conditions)


Sizes=MedianNorm(GeneExampleData)

EBSeqHMMGeneOut=EBSeqHMMTest(Data=GeneExampleData, sizeFactors=Sizes, Conditions=Conditions,
                             UpdateRd=5)

示例数据这里选用了5个时间点，那么意味着隐含序列为 I 有5个；而观测序列 O 有4个（观测序列为SgΔ）

# 核心代码：
## 其中Pi0为初始的状态概率
## LogB为利用Baum-Welch算法估计的发射概率矩阵B(对数化)

AlllogPost <- sapply(1:nrow(AllPosi),function(i){
    logPi0.use <- log(Pi0[AllPosi[i,1]])
    
    # 这里的每一个 LogB[[j]] 代表的是4个时间点中每一个时间点 t 对应的发射矩阵
    # AllPosi[i,j] 中的 i 代表81中组合的其中一中组合; j 代表4个时间点对应的其中一种状态(up，down和EE)
    # 取出每个时间点 t 对应其中一种状态的概率
    logB.use <- sapply(1:NumPoints,function(j)LogB[[j]][AllPosi[i,j],])

    # gene by time
    ## 取出每个时间点 t 对应其中一种状态的概率后，进行加和作为81种组合其中一种组合的概率
    logB.sumtime <- rowSums(logB.use)
    if(homo==TRUE)tran.use <- sapply(1:(NumPoints-1),function(j)Tran[as.numeric(AllPosi[i,j]),as.numeric(AllPosi[i,j+1])])
    if(homo==FALSE)tran.use <- sapply(1:(NumPoints-1),function(j)Tran[as.numeric(AllPosi[i,j]),as.numeric(AllPosi[i,j+1]),j])
    logtran.sum=sum(log(tran.use))
    # 加和作为最终的概率
    ## logPi0.use为起始的概率的对数值(默认为 1/3)
    ## tran.use为转移概率矩阵
    ## logtran.sum为转移概率矩阵的对数和
    ## logB.sumtime为发射概率矩阵的对数和
    logPi0.use+logB.sumtime+logtran.sum 
})

1. logPi0.use+logB.sumtime+logtran.sum 对数加和相当于乘积，表征的是对于 gene g 81中路径组合的概率

2. logPi0.use为起始的概率的对数值(默认为 1/3)，其中 Pi0 为初始的状态概率，对应三种状态

   
   

3. tran.use为转移概率矩阵，因为有5个时间点，则有4个SgΔ，所以有3次转移，每个基因为三维向量

   
   
   

4. logtran.sum为转移概率矩阵的对数和，因为有5个时间点，则有4个SgΔ，所以有3次转移（下图黑色横箭头）；

   
   
5. LogB 为发射概率矩阵，发射概率矩阵，因为有4个SgΔ（4个时间点），所以有4次发射，每个基因为四维的向量（图中展示的是取了log后的值），st1，st2，st3代表的是up，down和EE；而该发射概率矩阵是基于每个sample的基因表达情况，根据Baum-Welch算法估计出来的；其中行为状态，列为基因；它表征的意义是在已知隐含层Xg的条件下，观测到三种状态（up，down，EE）其中一个的概率
6. logB.use 表征的意义是对于81种组合的每一种组合，在已知在隐含层（Xg,1 Xg,2 Xg,3 Xg,4）其中一个的条件下，观测值为 （up，down，EE）其中一个的概率，每一列代表其中一种路径组合的三种观测状态（up，down，EE） 比方说第一条路经为"up，up，down，EE"。其中 LogB 为发射概率矩阵，那么上图第一列代表在第一个隐含层中，观测到 "up" 的那一行发射概率矩阵的数值 （行为状态，列为基因）；上图第二列代表在第二个隐含层中，观测到 "up" 的那一行发射概率矩阵的数值 ；上图第三列代表在第三个隐含层中，观测到 "down" 的那一行发射概率矩阵的数值；上图第四列代表在第四个隐含层中，观测到 "EE" 的那一行发射概率矩阵的数值

7. logB.sumtime为发射概率矩阵的对数和，因为有4个SgΔ，所以有4次发射（下图黑色斜箭头），对数求和即原矩阵相乘，同时发生事件的概率（路径概率，积事件）

   
   

其中：

发射概率矩阵参数的学习迭代公式如上图，代表HMM模型当前参数的估计值；P(O,i_t = j | ) 代表当前参数的条件下，观测值 O 属于隐含层 j 的概率，而这概率的初始值可以分别统计三种状态"up，down，EE"在各个观测层（该例子中有四个观测层；SgΔ1，SgΔ2，SgΔ3，SgΔ4）出现的频率（即每个观测层出现up，down，EE三种状态的基因频率；比方说一共100个基因，up状态一共有30个，那么up发射概率的初始值为0.3）

对于结果
LogB为Xg,2，Xg,3，Xg,4，Xg,5这几个状态的发射矩阵（一共有4个发射矩阵）；而 st1，st2，st3代表的是up，down和EE：

其中 1，2，3，4代表四个时间点；I1，I2...代表基因

AllPosi 为四个时间点状态观测值的可能性集合，3的4次方一共有81中组合；每一列分别代表 t1，t2，t3，t4 对应的状态观测值；每一行代表每一种组合：

最终的结果AlllogPost ：

每一列代表81种组合，每一行代表基因，数值代表其中一种组合的概率，而最佳路径，就是选取其中概率最大的那种组合即可

小节模型逻辑思想

1. 作者首先定义出三种状态"up，down，EE"，利用这三种状态
2. 然后用Baum-Welch算法估算当HMM的参数 λ 为多少时，观测到某一状态的概率最大，并且求出这个概率
3. 将每个gene的4个观测状态的概率相乘，那么概率最大的作为最佳路径

值得注意的是LogB 为发射概率矩阵，发射概率矩阵，因为有4个SgΔ，所以有4次发射，每个基因为四维的向量（图中展示的是取了log后的值），而该发射概率矩阵是基于每个sample的基因表达情况，根据Baum-Welch算法估计出来的；而 logB.use 表征的意义是已知在隐含层（Xg,1 Xg,2 Xg,3 Xg,4）其中一个的条件下，观测值为 （up，down，EE）其中一个的概率