R/dRegulation.R

In scTenifoldKnk: In-Silico Knockout Experiments from Single-Cell Gene Regulatory Networks

#' @importFrom stats dist p.adjust pchisq
#' @importFrom MASS boxcox
#' @import pbapply
#' @import utils
#' @import methods

dRegulation <- function(manifoldOutput, gKO) {
  geneList <- rownames(manifoldOutput)
  geneList <- geneList[grepl("^X_", geneList)]
  geneList <- gsub("^X_", "", geneList)
  nGenes <- length(geneList)

  eGenes <- nrow(manifoldOutput) / 2

  eGeneList <- rownames(manifoldOutput)
  eGeneList <- eGeneList[grepl("^Y_", eGeneList)]
  eGeneList <- gsub("^Y_", "", eGeneList)

  if (nGenes != eGenes) {
    stop("Number of identified and expected genes are not the same")
  }
  if (!all(eGeneList == geneList)) {
    stop("Genes are not ordered as expected. X_ genes should be followed by Y_ genes in the same order")
  }

  dMetric <- sapply(seq_len(nGenes), function(G) {
    X <- manifoldOutput[G, ]
    Y <- manifoldOutput[(G + nGenes), ]
    I <- rbind(X, Y)
    O <- stats::dist(I)
    O <- as.numeric(O)
    return(O)
  })

  ### BOX-COX
  lambdaValues <- seq(-2, 2, length.out = 1000)
  lambdaValues <- lambdaValues[lambdaValues != 0]
  BC <- try(MASS::boxcox(dMetric[dMetric > 0] ~ 1, plot = FALSE, lambda = lambdaValues), silent = TRUE)
  if (inherits(BC, "try-error")) {
    nD <- dMetric
  } else {
    BC <- BC$x[which.max(BC$y)]
    if (BC < 0) {
      nD <- 1 / (dMetric^BC)
    } else {
      nD <- dMetric^BC
    }
  }
  Z <- scale(nD)
  dOut <- data.frame(
    gene = geneList,
    distance = dMetric,
    Z = Z
  )
  dOut <- dOut[order(dOut$distance, decreasing = TRUE), ]
  FC <- (dOut$distance^2) / mean((dOut$distance[-seq_len(length(gKO))]^2))
  pValues <- pchisq(q = FC, df = 1, lower.tail = FALSE)
  pAdjusted <- p.adjust(pValues, method = "fdr")
  dOut$FC <- FC
  dOut$p.value <- pValues
  dOut$p.adj <- pAdjusted
  dOut <- as.data.frame.array(dOut)
  return(dOut)
}