R/generate_uORF.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' generate_uORF
#'
#' \code{generate_uORF} generates a single trajecture of observed chain x.
#'
#' @export
#' @param eL    a vector. Expected lengths of (5U, uE, 5U2, E, 3U).
#' @param eH    a vector. Expected heights of (5U, uI, uE, uT, 5U2, I, E, T, 3U).
#'              This is the mean parameter in gamma distribution.
#' @param v     a vector. Shape parameter in gamma distribution.
#' @param E     a scalar. Normalizing constant for the observed chain x.
#' @param r     a scalar. Ratio of chians with uORF and non-uORF.
#'
#' @return      A list:
#'                  \itemize{
#'                    \item z    :  a vector. Hidden chain z
#'                    \item x    :  a vector. Observed chain x
#'                    \item RNA  :  a 0-1 vector. 1 if next 3-base is stop codon
#'                    \item trans:  a vector. (rho, rho_u, delta)
#'                    \item m    :  a vector. mean parameter in gamma distribution for each state
#'                    \item v    :  a vector. Shape parameter in gamma distribution.
#'                    \item E    :  a scalar. Normalizing constant for the observed chain x
#'                    \item stp  :  a vector. (theta_u, theta) transition probability used
#'                    }

#' @examples
#' set.seed(2019)
#' eL <- c(24, 36, 24, 96, 20)
#' eH <- c(5, 50,35,20, 40,20,10, 50,35,20,
#'         5, 60,35,10, 55,25,10, 60,35,10, 5)
#' v <- rep(1000,21)
#' E <- 10
#' r <- 1
#' df <- generate_uORF(eL, eH, v, E, r)
#' length(df$x)
#' barplot(df$x, names.arg=df$z)

generate_uORF <- function(eL, eH, v, E, r){
  M <- 21
  # transition matrix
  rho <- 1/eL[1]/(r+1)
  rho_u <- r*rho
  delta <- 1/eL[3]
  theta_u <- 1/eL[2]*3
  theta <- 1/eL[4]*3
  ter <- 1/eL[5]
  trans <- c(rho_u, rho, delta)
  A <- transprob(trans, stp=theta_u, stp2=theta)

  m <- eH
  nb_prob <- v/(m*E+v)

  z <- 1  # start with state 5'U
  x <- rnbinom(1, size=v[1], prob=nb_prob[1])
  RNA <- c()
  flg <- 1
  t <- 2
  while (flg){
    z <- c(z, sample.int(M, size=1, prob=A[z[t-1],]))
    x <- c(x, rnbinom(1, size=v[z[t]], prob=nb_prob[z[t]]))
    if (z[t] == 8 || z[t] == 18){
      RNA <- c(RNA, 1)
    }else{
      RNA <- c(RNA, 0)
    }
    if (z[t] == 21 && runif(1) < ter){
      flg <- 0
    }
    t <- t + 1
  }
  RNA <- c(RNA, 0)
  stp <- c(theta_u, theta)
  return(list(z=z, x=x, RNA=RNA, trans=trans, m=m, v=v, E=E, stp=stp))
}