R: HMM

HMM.R

## Demo for HMM
## Reproduce "What is HMM", Sean R Eddy, Nature Biotechnology 2004

S <- c(0,1,rep(0,3))
E <- c(0, 0.9, 0.1, 0, 0)
SS5 <- c(rep(0,3), 1, 0)
I <- c(rep(0,3), 0.9, 0.1)
N <- c(rep(0,4), 1)
myTransMtx <- rbind(S,E, SS5, I, N)
colnames(myTransMtx) <- rownames(myTransMtx)
myTransMtx
S <- rep(0,4)
E <- rep(1/4, 4)
SS5 <- c(0.05, 0, 0.95, 0)
I <- c(0.4, 0.1, 0.1, 0.4)
N <- c(rep(0, 4))
myEmisMtx <- rbind(S,E, SS5, I, N)
colnames(myEmisMtx) <- c("A", "C", "G", "T")
myEmisMtx

mySeq <- "CTTCATGTGAAAGCAGACGTAAGTCA"
myPath <- "EEEEEEEEEEEEEEEEEE5IIIIIII"

computeProb <- function(obsSeq, imgPath, transMtx, emisMtx, log=TRUE){
    obsSeq <- unlist(strsplit(toupper(obsSeq), ""))
    imgPath <- unlist(strsplit(imgPath, ""))
    imgPath[imgPath=="5"] <- "SS5"
    ## Previous is Start State
    prob <-  transMtx["S", imgPath[1]] * emisMtx[imgPath[1], obsSeq[1]]
    for (i in 2:length(obsSeq)){
        state_i <- imgPath[i]
        nucleotide_i <- obsSeq[i]
        prob_i <- transMtx[imgPath[i-1],state_i] * emisMtx[imgPath[i],nucleotide_i]
        prob <- prob * prob_i
    }
    ## End State
    prob <- prob * transMtx[imgPath[length(imgPath)],"N"]
    if (log){
        return (log2(prob)/log2(exp(1))) ## Why use e ?! Sean Eddy!!!
    } else {
        return (prob)
    }
}

giveMeMore <- function(obsSeq, transMtx, emisMtx ){
    sumProb <- 0
    allPath <- NULL
    allProb <- NULL
    print ("All State Path with non-zero Probability: ")
    for (i in 1:24) {
        path_i <- c(rep("E", i), "5", rep("I", 25-i))
        path_i <- paste(path_i, collapse="")
        allPath[i] <- path_i
        
        p <- computeProb(obsSeq, path_i,transMtx, emisMtx,log=FALSE)
        allProb[i] <- p 
        if (p != 0){
            print (path_i)
            print (p)
            sumProb <- sumProb + p 
        }
    }
    maxProb <- max(allProb, na.rm=TRUE)
    max <- which(allProb == maxProb)
    cat ("Best State Path and original Sequence", "\n")
    cat (allPath[max], "\n")
    cat (obsSeq,"\n")
    cat (paste("Maximum Probability: ", maxProb, sep=""),"\n")
    cat (paste("Maximum Prob (log):  ", log2(maxProb)/log2(exp(1)), sep=""),"\n")
    cat (paste("Posterior Decoding:  ", maxProb/sumProb, sep=""),"\n")
    return (allPath[max])
}