tamlabscpipeline: A basic package for scRNAseq data analysis

Documented in checkdependencies

#check deps ---------------------------
.onAttach <- function(libname, pkgname) {

  packageStartupMessage('\nThank you for using this package. Checking dependencies:')

  pkgs <- c("Seurat", "sctransform", "ecp", "tidyverse", "cowplot", "DoubletFinder", "MAST", "fgsea")
  rv <- c("3.1.2", "0.2.1", "3.1.2", "1.3.0", "1.0.0", "2.0.2", "1.12.0", "1.12.0")

  suggests <- c("MAST")

  for(i in 1:length(pkgs)){

    pkg <- pkgs[i]
    msg <- paste0(pkg,': ',  rv[i])

    if(pkg %in% suggests){
      msg <- paste0(msg, ' (suggested);')
    }


    #test if pkg installed
    if(pkg %in% rownames(installed.packages()) ){



      #check if version installed is outdated
      if( packageVersion(pkg) < rv[i] ){
        msg <- paste0(msg, ' ', sprintf('‽'),
                      '\n  Warning: tamlabscpipeline is untested below version: ', rv[i])
      }

      #check if installed version is above ours
      if( packageVersion(pkg) > rv[i] ){
        msg <- paste0(msg, ' ', sprintf('‽'),
                      '\n  Warning: tamlabscpipeline is untested above version: ', rv[i],
                      "\n  Please contact Alex (alexander.ferrena@einsteinmed.org)")
      }

      #if installed version == latest we have:
      if( packageVersion(pkg) == rv[i] ){
        msg <- paste0(msg, ' installed ', sprintf('\u2714'))
      }


    } else{

      msg <- paste0(msg,' ',  sprintf('\u2715'), '\n\tWarning: pkg not installed or detected on your system ')
    }

    packageStartupMessage(msg)
  }

}

#' Check if Dependencies are installed
#'
#' A function to check if dependencies are installed. Installed dependncies must be loaded manually.
checkdependencies <- function() {

  message('\nThank you for using this package. Checking dependencies:')

  pkgs <- c("Seurat", "sctransform", "ecp", "tidyverse", "cowplot", "DoubletFinder", "MAST", "fgsea")
  rv <- c("3.1.2", "0.2.1", "3.1.2", "1.3.0", "1.0.0", "2.0.2", "1.12.0", "1.12.0")

  suggests <- c("MAST")

  for(i in 1:length(pkgs)){

    pkg <- pkgs[i]
    msg <- paste0(pkg,': ',  rv[i])

    if(pkg %in% suggests){
      msg <- paste0(msg, ' (suggested);')
    }

    #test if pkg installed
    if(pkg %in% rownames(installed.packages()) ){



      #check if version installed is outdated
      if( packageVersion(pkg) < rv[i] ){
        msg <- paste0(msg, ' ', sprintf('‽'),
                      '\n  Warning: tamlabscpipeline is untested below version: ', rv[i])
      }

      #check if installed version is above ours
      if( packageVersion(pkg) > rv[i] ){
        msg <- paste0(msg, ' ', sprintf('‽'),
                      '\n  Warning: tamlabscpipeline is untested above version: ', rv[i],
                      "\n  Please contact Alex (alexander.ferrena@einsteinmed.org)")
      }

      #if installed version == latest we have:
      if( packageVersion(pkg) == rv[i] ){
        msg <- paste0(msg, ' installed ', sprintf('\u2714'))
      }


    } else{

      msg <- paste0(msg,' ',  sprintf('\u2715'), '\n\tWarning: pkg not installed or detected on your system ')
    }

    message(msg)
  }

}

#grubbs test for outliers ------------------------------------
#' a test for outliers assuming normal distribution.
#' based on the grubbs.test() function from the outliers pkg by Lukasz Komsta
grubbs <- function (x, type = 10, opposite = FALSE, two.sided = FALSE)
{
  if (sum(c(10, 11, 20) == type) == 0)
    stop("Incorrect type")
  DNAME <- deparse(substitute(x))
  x <- sort(x[complete.cases(x)])
  n <- length(x)
  if (type == 11) {
    g <- (x[n] - x[1])/sd(x)
    u <- var(x[2:(n - 1)])/var(x) * (n - 3)/(n - 1)
    pval = 1 - pgrubbs(g, n, type = 11)
    method <- "Grubbs test for two opposite outliers"
    alt = paste(x[1], "and", x[n], "are outliers")
  }
  else if (type == 10) {
    if (xor(((x[n] - mean(x)) < (mean(x) - x[1])), opposite)) {
      alt = paste("lowest value", x[1], "is an outlier")
      o <- x[1]
      d <- x[2:n]
    }
    else {
      alt = paste("highest value", x[n], "is an outlier")
      o <- x[n]
      d <- x[1:(n - 1)]
    }
    g <- abs(o - mean(x))/sd(x)
    u <- var(d)/var(x) * (n - 2)/(n - 1)
    pval <- 1 - pgrubbs(g, n, type = 10)
    method <- "Grubbs test for one outlier"
  }
  else {
    if (xor(((x[n] - mean(x)) < (mean(x) - x[1])), opposite)) {
      alt = paste("lowest values", x[1], ",", x[2], "are outliers")
      u <- var(x[3:n])/var(x) * (n - 3)/(n - 1)
    }
    else {
      alt = paste("highest values", x[n - 1], ",", x[n],
                  "are outliers")
      u <- var(x[1:(n - 2)])/var(x) * (n - 3)/(n - 1)
    }
    g <- NULL
    pval <- pgrubbs(u, n, type = 20)
    method <- "Grubbs test for two outliers"
  }
  if (two.sided) {
    pval <- 2 * pval
    if (pval > 1)
      pval <- 2 - pval
  }
  RVAL <- list(statistic = c(G = g, U = u), alternative = alt,
               p.value = pval, method = method, data.name = DNAME)
  class(RVAL) <- "htest"
  return(RVAL)
}


qgrubbs <- function (p, n, type = 10, rev = FALSE)
{
  if (type == 10) {
    if (!rev) {
      return(((n - 1)/sqrt(n)) * sqrt(qt((1 - p)/n, n -
                                           2)^2/(n - 2 + qt((1 - p)/n, n - 2)^2)))
    }
    else {
      s <- (p^2 * n * (2 - n))/(p^2 * n - (n - 1)^2)
      t <- sqrt(s)
      if (is.nan(t)) {
        res <- 0
      }
      else {
        res <- n * (1 - pt(t, n - 2))
        res[res > 1] <- 1
      }
      return(1 - res)
    }
  }
  else if (type == 11) {
    if (!rev) {
      return(sqrt((2 * (n - 1) * qt((1 - p)/(n * (n - 1)),
                                    n - 2)^2)/(n - 2 + qt((1 - p)/(n * (n - 1)),
                                                          n - 2)^2)))
    }
    else {
      q <- p
      p <- vector()
      for (i in 1:length(q)) {
        if (q[i] > qgrubbs(0.9999, n, type = 11)) {
          pp <- 1
        }
        else if (q[i] < qgrubbs(2e-16, n, type = 11)) {
          pp <- 0
        }
        else {
          f <- function(x, q, n) {
            qgrubbs(x, n, type = 11) - q
          }
          pp <- uniroot(f, c(0.001, 0.9999), q = q[i],
                        n = n)$root
        }
        p <- c(p, pp)
      }
      return(p)
    }
  }
  else {
    if (n > 30)
      stop("n must be in range 3-30")
    pp <- c(0.01, 0.025, 0.05, 0.1, 0.15, 0.2, 0.4, 0.6,
            0.8, 0.9, 0.95, 0.975, 0.99)
    gtwo <- c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
              NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
              NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
              NA, NA, 1e-05, 2e-04, 8e-04, 0.0031, 0.007, 0.013,
              0.055, 0.138, 0.283, 0.399, 0.482, 0.54, 0.589, 0.0035,
              0.009, 0.0183, 0.0376, 0.058, 0.078, 0.169, 0.276,
              0.41, 0.502, 0.571, 0.63, 0.689, 0.0186, 0.0349,
              0.0565, 0.0921, 0.124, 0.153, 0.257, 0.361, 0.478,
              0.562, 0.626, 0.674, 0.724, 0.044, 0.0708, 0.102,
              0.1479, 0.186, 0.217, 0.325, 0.423, 0.53, 0.605,
              0.659, 0.701, 0.743, 0.075, 0.1101, 0.1478, 0.1994,
              0.238, 0.271, 0.376, 0.468, 0.568, 0.635, 0.684,
              0.722, 0.763, 0.1082, 0.1492, 0.1909, 0.2454, 0.285,
              0.319, 0.42, 0.507, 0.598, 0.658, 0.703, 0.738, 0.776,
              0.1415, 0.1865, 0.2305, 0.2863, 0.326, 0.358, 0.455,
              0.537, 0.622, 0.678, 0.721, 0.753, 0.787, 0.1736,
              0.2212, 0.2666, 0.3226, 0.363, 0.394, 0.487, 0.564,
              0.624, 0.695, 0.734, 0.765, 0.797, 0.2044, 0.2536,
              0.2996, 0.3552, 0.393, 0.424, 0.514, 0.585, 0.66,
              0.709, 0.745, 0.773, 0.802, 0.2333, 0.2836, 0.3295,
              0.3843, 0.421, 0.451, 0.537, 0.605, 0.676, 0.721,
              0.755, 0.782, 0.809, 0.2605, 0.3112, 0.3568, 0.4106,
              0.447, 0.477, 0.558, 0.622, 0.689, 0.733, 0.765,
              0.789, 0.817, 0.2859, 0.3367, 0.3818, 0.4345, 0.469,
              0.497, 0.576, 0.639, 0.701, 0.742, 0.773, 0.797,
              0.822, 0.3098, 0.3603, 0.4048, 0.4562, 0.491, 0.518,
              0.593, 0.653, 0.713, 0.751, 0.78, 0.803, 0.826, 0.3321,
              0.3822, 0.4259, 0.4761, 0.511, 0.536, 0.609, 0.666,
              0.723, 0.76, 0.787, 0.808, 0.831, 0.353, 0.4025,
              0.4455, 0.4944, 0.526, 0.552, 0.622, 0.676, 0.732,
              0.767, 0.793, 0.814, 0.835, 0.3725, 0.4214, 0.4636,
              0.5113, 0.543, 0.567, 0.635, 0.688, 0.74, 0.774,
              0.799, 0.818, 0.838, 0.3909, 0.4391, 0.4804, 0.5269,
              0.559, 0.582, 0.647, 0.697, 0.748, 0.781, 0.805,
              0.823, 0.843, 0.408, 0.457, 0.496, 0.542, 0.571,
              0.594, 0.657, 0.706, 0.755, 0.786, 0.81, 0.828, 0.847,
              0.425, 0.474, 0.512, 0.556, 0.584, 0.606, 0.668,
              0.715, 0.762, 0.792, 0.815, 0.834, 0.85, 0.442, 0.486,
              0.524, 0.568, 0.596, 0.618, 0.677, 0.723, 0.769,
              0.797, 0.819, 0.836, 0.853, 0.453, 0.5, 0.538, 0.581,
              0.608, 0.628, 0.686, 0.73, 0.774, 0.802, 0.823, 0.84,
              0.857, 0.466, 0.511, 0.547, 0.589, 0.616, 0.637,
              0.693, 0.736, 0.779, 0.807, 0.827, 0.843, 0.86, 0.482,
              0.525, 0.561, 0.601, 0.627, 0.647, 0.701, 0.743,
              0.784, 0.811, 0.83, 0.845, 0.861, 0.492, 0.536, 0.572,
              0.611, 0.636, 0.655, 0.709, 0.749, 0.789, 0.815,
              0.834, 0.849, 0.864, 0.505, 0.548, 0.583, 0.621,
              0.646, 0.664, 0.716, 0.755, 0.794, 0.819, 0.837,
              0.851, 0.866, 0.516, 0.558, 0.592, 0.629, 0.654,
              0.672, 0.722, 0.76, 0.798, 0.822, 0.84, 0.854, 0.869,
              0.528, 0.568, 0.602, 0.638, 0.661, 0.679, 0.728,
              0.765, 0.802, 0.826, 0.842, 0.856, 0.87)
    dim(gtwo) <- c(13, 30)
    if (!rev)
      res <- qtable(p, pp, gtwo[, n])
    else res <- qtable(p, gtwo[, n], pp)
    res[res < 0] <- 0
    res[res > 1] <- 1
    return(res)
  }
}

pgrubbs <- function (q, n, type = 10)
{
  qgrubbs(q, n, type, rev = TRUE)
}