R/generate_ribo.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' generate_ribo
#'
#' \code{generate_ribo} 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 p1    a scalar. Probability of uORF is translated.
#' @param p2    a scalar. Probability of main region is translated.
#'
#' @return      A list:
#'                  \itemize{
#'                    \item z_base:  a vector. Hidden chain z
#'                    \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 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 p1, p2
#'                    }
#'
#'
#' @examples
#' set.seed(2019)
#' eL <- c(20, 36, 20, 96, 10)
#' eH <- c(5, 45,30,15, 40,20,10, 45,30,15,
#'         5, 90,75,30, 60,50,20, 90,75,30, 5)
#' v <- rep(10,21)
#' p1=0.8; p2=0.8
#' E <- 1
#'
#' ribo <- generate_ribo(eL, eH, v, E, p1,p2)

generate_ribo <- function(eL, eH, v, E, p1, p2){
  M <- 21
  A <- transprob_ribo(eL)
  ter <- 1/eL[5]

  # decide which type of trajectory
  type <- runif(1)
  if (type < p1*(1-p2)){
    eH[12:20] <- rep(eH[1], 9)
    v[12:20] <- rep(v[1], 9)
  } else if (type < p1+p2-2*p1*p2) {
    eH[2:10] <- rep(eH[1], 9)
    v[2:10] <- rep(v[1], 9)
  } else if (type < 1-p1*p2) {
    eH <- rep(eH[1], 21)
    v <- rep(v[1], 21)
  }

  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)

  z_base <- z
  if (type < p1*(1-p2)){
    z[z %in% 12:21] <- 11
  } else if (type < p1+p2-2*p1*p2) {
    z[z %in% 2:11] <- 1
  } else if (type < 1-p1*p2) {
    z <- rep(1, length(z))
  }

  return(list(z_base=z_base, z=z, x=x, RNA=RNA, m=m, v=v, E=E, p1=p1, p2=p2))
}