R/mts_omics.R

In MultiSEp: Predict Synthetic Lethality and Other Gene Dependency Relationships from Multiomics Data

mts_omics <- function(dataMatrix = NULL, dataMatrix2 = NULL,
                      mixModelClusters1 = NULL, mixModelClusters2 = NULL,
                      categorical1 = FALSE, categorical2 = FALSE,
                      genepairs = NULL,
                      SyntheticLethalityPrediction = TRUE,
                      p_adjustMethod = c("BY", "BH"), qVal,
                      directionality = c("depletion", "enrichment"),
                      include_reverse_pairs = TRUE,
                      effectsize = TRUE, effectsize_threshold = NULL,
                      verbose = FALSE,
                      cores) {

  # Validate a pre-computed cluster list: non-empty list of data frames
  # with at least 3 columns (column 1 = sample id, column 3 = cluster int).
  validateClusterList <- function(x, name) {
    if (!is.list(x) || length(x) == 0) {
      stop(paste0(name, " must be a non-empty list of data frames produced by ",
                  "mts_mixModelCluster() or mts_formatMatrix()."))
    }
    ok <- vapply(x, function(e) is.data.frame(e) && ncol(e) >= 3, logical(1))
    if (!all(ok)) {
      stop(paste0(name, " entries must be data frames with at least 3 columns."))
    }
    invisible(TRUE)
  }

  # Returns a (possibly NULL) cluster list given the
  # raw data matrix, pre-computed clusters, and categorical flag.
  resolveSide <- function(dm, clusters, categorical, side_label,
                          required = FALSE) {
    if (!is.null(dm) && !is.null(clusters)) {
      stop(paste0("Provide only one of dataMatrix", side_label,
                  " or mixModelClusters", side_label, " (not both)."))
    }
    if (!is.null(clusters)) {
      validateClusterList(clusters,
                          paste0("mixModelClusters", side_label))
      return(clusters)
    }
    if (!is.null(dm)) {
      if (!is.data.frame(dm)) {
        stop(paste0("dataMatrix", side_label, " must be a data frame"))
      }
      if (isTRUE(categorical)) {
        if (verbose) message("Side ", side_label,
                             ": categorical input - using mts_formatMatrix()")
        return(mts_formatMatrix(matrix = dm, cores = cores))
      } else {
        if (verbose) message("Side ", side_label,
                             ": continuous input - using mts_mixModelCluster()")
        return(mts_mixModelCluster(dataMatrix = dm, cores = cores))
      }
    }
    if (required) {
      stop(paste0("No input provided for side ", side_label,
                  ": supply dataMatrix", side_label,
                  " or mixModelClusters", side_label, "."))
    }
    return(NULL)
  }

  ### Input checks
  # genepairs
  if (!is.null(genepairs)) {
    if (!(is.data.frame(genepairs) || is.matrix(genepairs))) {
      stop("The genepairs must be either a dataframe or a matrix; genepairs can be auto-generated by omitting the argument or setting genepairs to NULL.")
    }
    genepairs <- as.data.frame(genepairs)
  }

  # cores
  if (missing(cores)) {
    cores <- 1
    message("1 core selected")
  } else if (!is.numeric(cores) || cores < 1) {
    stop("cores should be >= 1")
  } else {
    message(cores, " cores selected")
  }

  # qVal
  if (missing(qVal)) {
    qVal <- 0.05
    if (verbose) message("Q-value Threshold not entered, using 0.05")
  } else {
    if (verbose) message(paste("Q-value Change Threshold: ", qVal, sep = ""))
  }

  # p_adjustMethod
  if (missing(p_adjustMethod)) {
    p_adjustMethod <- "BY"
    if (verbose) message("P-value correction using BY")
  } else {
    if (!(p_adjustMethod %in% c("BY", "BH", "fdr"))) {
      stop("Invalid p-value correction method. Please choose either 'BY', 'BH', or 'fdr'")
    }
    if (verbose) message("P-value correction using ", p_adjustMethod)
  }

  # SyntheticLethalityPrediction
  if (!(identical(SyntheticLethalityPrediction, TRUE) ||
        identical(SyntheticLethalityPrediction, FALSE))) {
    stop("Invalid input for SyntheticLethalityPrediction: please enter either TRUE or FALSE.")
  }
  if (verbose) {
    if (SyntheticLethalityPrediction) {
      message("Synthetic Lethality Prediction")
    } else {
      message("Predictions for all possible Gene Dependency Patterns")
    }
  }

  # effectsize
  if (!(identical(effectsize, TRUE) || identical(effectsize, FALSE))) {
    stop("Invalid input for effectsize: please enter either TRUE or FALSE.")
  }

  # categorical flags
  if (!(identical(categorical1, TRUE) || identical(categorical1, FALSE))) {
    stop("Invalid input for categorical1: please enter either TRUE or FALSE.")
  }
  if (!(identical(categorical2, TRUE) || identical(categorical2, FALSE))) {
    stop("Invalid input for categorical2: please enter either TRUE or FALSE.")
  }

  # directionality
  if (missing(directionality)) {
    directionality <- "depletion"
  } else if (!(directionality %in% c("enrichment", "depletion"))) {
    stop("Invalid value of directionality. Please specify either 'enrichment' or 'depletion'.")
  }

  ### Build clusters if needed
  # Determine if bidirectional analysis occurs (if data type 2 is input)
  bidirectional <- !is.null(dataMatrix2) || !is.null(mixModelClusters2)

  # Align samples when both input data types are data matrices
  if (bidirectional &&
      !is.null(dataMatrix) && !is.null(dataMatrix2)) {
    shared_cols <- intersect(colnames(dataMatrix), colnames(dataMatrix2))
    dataMatrix  <- dataMatrix[,  shared_cols, drop = FALSE]
    dataMatrix2 <- dataMatrix2[, shared_cols, drop = FALSE]
  }

  if (verbose) message("Resolving cluster inputs")

  mixModelClusters1 <- resolveSide(dataMatrix, mixModelClusters1,
                                   categorical1, "1", required = TRUE)
  mixModelClusters1 <- prefilterMixModelClusters(mixModelClusters1)

  if (bidirectional) {
    mixModelClusters2 <- resolveSide(dataMatrix2, mixModelClusters2,
                                     categorical2, "2", required = FALSE)
    mixModelClusters2 <- prefilterMixModelClusters(mixModelClusters2)
  } else {
    mixModelClusters2 <- NULL
    include_reverse_pairs <- FALSE
  }

  ### Predict gene dependency relationships
  findGDR <- mts_patternDetection(mixModelClusters1 = mixModelClusters1,
                                  mixModelClusters2 = mixModelClusters2,
                                  genepairs = genepairs,
                                  SyntheticLethalityPrediction = SyntheticLethalityPrediction,
                                  p_adjustMethod = p_adjustMethod,
                                  qVal = qVal,
                                  directionality = directionality,
                                  include_reverse_pairs = include_reverse_pairs,
                                  effectsize = effectsize,
                                  effectsize_threshold = effectsize_threshold,
                                  verbose = verbose,
                                  cores = cores)
  return(findGDR)
}