optimalPath: optimalPath

In shimlab/riboHMM2: Analysis tool for Ribo-seq

Description

optimalPath performs Viterbi Algorithm on a single observed trajectures x to predict the optimal path according to the rule, Maximum A Posteriori (MAP).

Usage

optimalPath(x, RNA, E, v, m, trans)

Arguments

x	vector of single observed trajectures
RNA	0-1 vector. 1 if next 3-base is stop codon
E	scalar. Normalizing constant for the observed chain x.
v	vector of shape parameter in gamma distribution (from estimate())
m	vector of mean parameter in gamma distribution (from estimate())
trans	vector c(rho_u, rho, delta) (from estimate())

Value

List containing:

• the optimal path optPth

• the corresponding log-likelihood loglik

Examples

k <- 5  # number of transcripts from each group (see data ?uORF)
ii <- c()
for (i in 1:5){
  ii <- c(ii, ((i-1)*100+1):((i-1)*100+k))
}
print(ii)

# estimate parameters
df <- uORF[ii]
num <- length(df)
X <- RNA <- list(); E <- c()
for (i in 1:num){
  X[[i]] <-  df[[i]]$x
  RNA[[i]] <- df[[i]]$RNA
  E[i] <- df[[i]]$E
}
init <- init_generate(df, cutoff=10, r=4)
est <- estimate(X, RNA, E, init, tol=0.01)
print(est)

# using estimated parameters in \code{optimalPath}
res <- matrix(0, num ,3)
colnames(res) <- c("rate", "TRUE", "FASLE")
for (i in 1:num){
  n <- length(df[[i]]$x)
  optPth <- optimalPath(df[[i]]$x, df[[i]]$RNA, E[i], est$v_hat, est$m_hat, est$trans_hat)$optPth
  res[i,2] <- sum(optPth == df[[i]]$z)
  res[i,3] <- n - res[i,2]
  res[i,1] <- res[i,2]/n
}
print(res)

shimlab/riboHMM2 documentation built on May 19, 2019, 6:23 p.m.