R/mts_patternDetection.R

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

mts_patternDetection <- function(genepairs = NULL, mixModelClusters1, mixModelClusters2 = NULL, 
                                 SyntheticLethalityPrediction = TRUE, 
                                 p_adjustMethod = c("BY","BH"), qVal, 
                                 effectsize = TRUE,  
                                 effectsize_threshold = NULL, include_reverse_pairs = FALSE,
                                 directionality = c("depletion", "enrichment"),  verbose = FALSE, 
                                 cores) {
  
  # input check - 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")
  }
  
  # input check - ensure mixModelClusters1 is a list
  if (missing(mixModelClusters1) || !is.list(mixModelClusters1)) {
    stop("mixModelClusters1 must be a list of data frames generated from the mts_mixModelCluster() function")
  }
  
  # input check - if mixModelClusters2 is provided, ensure is a list
  if (!is.null(mixModelClusters2)) {
    if (!is.list(mixModelClusters2)) {
      stop("mixModelClusters2 must be a list of data frames generated from the mts_mixModelCluster() function")
    }
  }
  
  # input check - bidirectional analysis and remove atomic vectors
  if (!is.null(mixModelClusters2)) {
    bidirectionalAnalysis <- TRUE
    if (verbose) message("Performing Bidirectional Analysis")
    mixModelClusters1 <- prefilterMixModelClusters(mixModelClusters1)
    mixModelClusters2 <- prefilterMixModelClusters(mixModelClusters2)
    filtered_lists = GMM_CCL(mixModelClusters1, mixModelClusters2)
    mixModelClusters1 = filtered_lists[[1]]
    mixModelClusters2 = filtered_lists[[2]]
  } else {
    bidirectionalAnalysis <- FALSE
    mixModelClusters1 <- prefilterMixModelClusters(mixModelClusters1)
    filtered_lists = GMM_CCL(mixModelClusters1)
    mixModelClusters1 = filtered_lists[[1]]
    if (verbose) message("Performing One Directional Analysis")
  }
  
  # input check - generation of genepairs
  if (is.null(genepairs)) {
    genepairs <- generateGenePairs(mixModelClusters1, mixModelClusters2,
                                   bidirectionalAnalysis = bidirectionalAnalysis,
                                   include_reverse_pairs = include_reverse_pairs,
                                   verbose = verbose, cores=cores)
  } else {
    # input check - ensure user provided genepairs are a dataframe or matrix
    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)
  }
  
  # input check - choose the effect-size threshold 
  #___ effect size threshold handling ___ 
  if (is.null(effectsize_threshold)) {
    if (bidirectionalAnalysis) {
      effectsize_threshold <- 0.5029284
    } else {
      effectsize_threshold <- 0.9621389  
    }
    if (verbose) message("Using default effect size threshold = ", effectsize_threshold)
  } else if (identical(effectsize_threshold, FALSE)) {
    effectsize <- FALSE
    if (verbose) message(("Not filtering for effect size"))
  } else if (is.numeric(effectsize_threshold) && effectsize_threshold >= 0 && effectsize_threshold <= 1) {
    if (verbose) message("Using user-supplied effect size threshold =", effectsize_threshold)
  } else {
    stop("`effectsize_threshold` must be FALSE or a number between 0 and 1")
  }
  
  
  # input check -  decide whether to perform SL analysis
  if (!(identical(SyntheticLethalityPrediction, TRUE) || identical(SyntheticLethalityPrediction, FALSE))) {
    stop("Invalid input for SyntheticLethalityPrediction: please enter either TRUE or FALSE.")
  }
  allPatterns <- !SyntheticLethalityPrediction # if TRUE then set to FALSE for the cluster proportion calculation
  if (SyntheticLethalityPrediction) {
    if (verbose)  message("Synthetic Lethality Prediction")
  } else {
    if (verbose)  message("Predictions for all possible Gene Dependency Patterns")
  }
  
  
  # input check - selection of multiple test correction  
  if (missing(p_adjustMethod)) {
    p_adjustMethod <- "BY"
    if (verbose)  message("P-value correction method not entered, 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("Adjusting p-values for multiple comparisons using: ", p_adjustMethod)
  }
  
  # input check - q-value filter 
  if(missing(qVal))
  {
    qVal <- 0.05
    if (verbose) message("Q-value Threshold not entered, using 0.05")
  } else {
    if (verbose) message("Q Value Threshold: ", qVal)
  }
  
  # input check - determining the directionality of predictions 
  if (missing(directionality)) {
    # default to depletion
    directionality <- "less"
    if (verbose) message("Predicting depletion")
  } else {
    if (directionality == "depletion") {
      directionality <- "less"
      if (verbose) message("Predicting depletion")
    } else if (directionality == "enrichment") {
      directionality <- "greater"
      if (verbose) message("Predicting enrichment")
    } else {
      stop("Directionality specified is invalid: Please choose either 'depletion' or 'enrichment'")
    }
  }
  
  results <- data.frame(Gene1 = character(), Gene2 = character(), Actual_Count = integer(), Cluster_Combinations = integer(), Sample_Count = integer(), Expected_Count = double(),
                        p_value = double(), q_value = double(), stringsAsFactors = FALSE)
  
  if (bidirectionalAnalysis) {
    ### _____________________  bidirectional analysis (2-GMM)  _________________________________________
    # find unique genes and match to indices 
    genes1 <- unique(genepairs$Gene1)
    genes2 <- unique(genepairs$Gene2)
    idx1 <- match(genes1, names(mixModelClusters1))
    idx2 <- match(genes2, names(mixModelClusters2))
    
    #_____  calculate proportions for mixModelClusters1 _______
    prop_cluster1 <- pbmclapply(idx1, function(i) {
      GeneClusterProportions(mixModelClusters1, i, cores, allPatterns)
    }, mc.cores = cores)
    
    #_____  calculate proportions for mixModelClusters2 _______
    prop_cluster2 <- pbmclapply(idx2, function(i) {
      GeneClusterProportions(mixModelClusters2, i, cores, allPatterns)
    }, mc.cores = cores)
    
    #_____  format proportions to rows together and loop through genepairs  _______
    prop_cluster1 <- prop_cluster1[!sapply(prop_cluster1, is.null)]
    prop_cluster2 <- prop_cluster2[!sapply(prop_cluster2, is.null)]
    
    if (length(prop_cluster1) == 0) {
      if (verbose) message("No valid genes found to compute cluster proportions.")
      return(results[0,])  
    }
    
    if (length(prop_cluster2) == 0) {
      if (verbose) message("No valid genes found to compute cluster proportions.")
      return(results[0,])  
    }
    
    prop_cluster1 <- do.call("rbind", prop_cluster1)
    prop_cluster2 <- do.call("rbind", prop_cluster2)
    
    pd_res = pbmclapply(1:nrow(genepairs), function(i) {
      gene1 <- genepairs[i, 1]
      gene2 <- genepairs[i, 2]
      
      # _____ check if genes exist in the respective GMM objects 
      gene1_present <- gene1 %in% names(mixModelClusters1)
      gene2_present <- gene2 %in% names(mixModelClusters2)
      
      missing_msgs <- c()
      if (!gene1_present) missing_msgs <- c(missing_msgs, paste0("'", gene1, "' not found in mixModelClusters1"))
      if (!gene2_present) missing_msgs <- c(missing_msgs, paste0("'", gene2, "' not found in mixModelClusters2"))
      
      if (length(missing_msgs) > 0) {
        if (verbose) {
          message("Skipping pair '", gene1, "' - '", gene2, "': ", paste(missing_msgs, collapse = " and "))
        }
        return(NULL)
      }
      
      ### _____________________   binomial  _________________________________________
      prop1_gene <- prop_cluster1[prop_cluster1$gene == gene1, ]
      prop2_gene <- prop_cluster2[prop_cluster2$gene == gene2, ]
      
      if (nrow(prop1_gene) == 0 || nrow(prop2_gene) == 0) {
        return(NULL)
      }
      y1 = TwoGMMBinomial(gene1, gene2, mixModelClusters1, mixModelClusters2, prop1_gene, prop2_gene, NULL,directionality, effectsize, verbose)
      return(y1)
    }, mc.cores=cores)
  } else {
    ### _____________________  one-directional analysis (1-GMM) _________________________________________
    allgenes <- unique(c(genepairs$Gene1, genepairs$Gene2))
    idxs  <- match(allgenes, names(mixModelClusters1))
    # drop any genes not found 
    valid <- !is.na(idxs)
    idxs <- idxs[valid]
    allgenes <- allgenes[valid]
    
    #_____  calculate proportions for mixModelClusters1 _______
    prop_cluster1 <- pbmclapply(idxs, function(i) {
      GeneClusterProportions(mixModelClusters1, i, cores, allPatterns)
    }, mc.cores = cores)
    
    
    #_____  format proportions to rows together and loop through genepairs  _______
    prop_cluster1 <- prop_cluster1[!sapply(prop_cluster1, is.null)]
    
    if (length(prop_cluster1) == 0) {
      if (verbose) message("No valid genes found to compute cluster proportions.")
      return(results[0,])  
    }
    
    prop_cluster1 <- do.call("rbind", prop_cluster1)
    
    pd_res = pbmclapply(1:nrow(genepairs), function(i) {
      gene1 <- genepairs[i, 1]
      gene2 <- genepairs[i, 2]
      
      # _____ check if genes exist in the GMM objects
      gene1_present <- gene1 %in% names(mixModelClusters1)
      gene2_present <- gene2 %in% names(mixModelClusters1)
      missing_msgs <- c()
      if (!gene1_present) missing_msgs <- c(missing_msgs, paste0("'", gene1, "' not found"))
      if (!gene2_present) missing_msgs <- c(missing_msgs, paste0("'", gene2, "' not found"))
      if (length(missing_msgs) > 0) {
        if (verbose) {
          message("Skipping pair '", gene1, "' - '", gene2, "': ", 
                  paste(missing_msgs, collapse = " and "), " in mixModelClusters1.")
        }
        return(NULL)
      }
      ### _____________________   binomial  _________________________________________
      tmp_results_binomial <- OneGMMBinomial(gene1, gene2, mixModelClusters1, prop_cluster1,  NULL, directionality, effectsize)
      return(tmp_results_binomial)
    }, mc.cores = cores) 
  }
  
  pd_res = Filter(function(x) is.data.frame(x) && nrow(x) > 0, pd_res)
  if (length(pd_res) == 0) {
    if (verbose) message("No valid gene pairs were analysed. Stopping.")
    return(results[0,])
  }
  results = do.call(rbind, pd_res)
  
  # _____error check: if no gene pairs were analysed 
  if (nrow(results) == 0) {
    if (verbose) message("No valid gene pairs were analysed. Stopping.")
    return(results[0, ])  
  }
  
  # _____multiple test correction 
  p_adjusted = p.adjust(results$p_value, method = p_adjustMethod)
  results$q_value = p_adjusted
  results = results[results$q_value <= qVal,] 
  
  # _____ apply effect size filter 
  if (effectsize) {
    if (!identical(effectsize_threshold, FALSE)) {
      if (!"Effect_Size" %in% names(results)) {
        warning("No Effect_Size column present; skipping threshold filter")
      } else {
        results <- results[results[["Effect_Size"]] >= effectsize_threshold, ]
        if (verbose) message("Filtered results with Effect_Size >=", effectsize_threshold)
      }
    }
  }
  return(results)
}