R/init_generate.R

In shimlab/riboHMM2: Analysis tool for Ribo-seq

#' init_generate
#'
#' \code{init_generate} generates a list of initial parameters
#' based on the ribo-seq.
#'
#' @export
#' @param df      a list of dataset of ribo-seq
#' @param cutoff  a scalar of cutoff point depending on histogram
#' @param r       a scalar. Ratio of chians with uORF and non-uORF.
#' @return        a list of initial parameters
#'
#' @examples
#' df <- uORF
#' init_generate(df, cutoff=10, r=4)

init_generate <- function(df, cutoff=10, r=4){
  X <- RNA <-list()
  E <- c()
  num <- length(df)
  for (i in 1:num){
    X[[i]] <-  df[[i]]$x
    RNA[[i]] <- df[[i]]$RNA
    E[i] <- df[[i]]$E
  }
  nx <- c()
  for (i in 1:num){
    nx <- c(nx, c(X[[i]]/E[i]))
  }
  m1 <- mean(nx[nx < 10])
  m6 <- mean(nx[nx > 10])
  init_m <- c(m1, rep(c(m6/3*4,m6,m6/3), 3), 
              m1, rep(c(m6/3*5,m6,m6/2), 3), m1)
  
  trans_reg <- nx[nx > 10]
  v <- m6^2/(var(trans_reg)-m6)
  init_v <- rep(v, 21)
  
  eL1 <- sum(nx < 10) /num/3
  rho <- 1/eL1/(r+1)
  rho_u <- r*rho
  delta <- rho + rho_u

  init <- list(v=init_v, m=init_m, trans=c(rho_u, rho, delta))
  return(init)
}