R/model_cpm_counts.R

In leandroroser/RNASeqFunctions: Several functions for rna-seq analysis

#' Linear models for cpm vs counts and counts vs cpm
#' @param x count matrix
#' @param n sample size for features, included in the model the model
#' @export

model_cpm_counts <- function(x,  n = 1000) {

  samp <- sample(seq_len(nrow(x)), n)

  x <- x + 1
  y <- normalize_counts_extra(x, "cpm")

  # if(!is.null(groups)) {
  #   x <- t(apply(x, 1, function(y) tapply(y, groups, mean)))
  #   y <- t(apply(y, 1, function(y) tapply(y, groups, mean)))
  # }

  x<- as.numeric(unlist(x[samp, ]))
  y <- as.numeric(unlist(y[samp,]))

  p_cpm <- lm(log(y) ~ log(x))
  p_counts <- lm(log(x) ~ log(y))


  predict_cpm <- function(counts) exp(p_cpm$coefficients[1] + (p_cpm$coefficients[2] * log(counts + 1)))
  predict_counts <- function(cpm) exp(p_counts$coefficients[1] + (p_counts$coefficients[2] * log(cpm))) - 1

  list(predict_cpm = predict_cpm, predict_counts = predict_counts,
       model_counts = p_counts, model_cpm = p_cpm)
}