R/comp-manipulate.R

In influential: Identification and Classification of the Most Influential Nodes

#=============================================================================
#
#    Computational manipulation of cells
#
#=============================================================================

#' Computational manipulation of cells
#'
#' This function works based on the SIRIR (SIR-based Influence Ranking) model and could be applied on the
#' output of the ExIR model or any other independent association network. For feature (gene/protein/etc.)
#' knockout the SIRIR model is used to remove the feature from the network and assess its impact on the
#' flow of information (signaling) within the network. On the other hand, in case of up-regulation a node similar
#' to the desired node is added to the network with exactly the same connections (edges) as of the original node.
#' Next, the SIRIR model is used to evaluate the difference in the flow of information/signaling after adding (up-regulating)
#' the desired feature/node compared with the original network. In case you are applying this function on the output of
#' ExIR model, you may note that as the gene/protein knockout would impact on the integrity of the under-investigation network
#' as well as the networks of other overlapping biological processes/pathways, it is recommended to select those features that
#' simultaneously have the highest (most significant) ExIR-based rank and lowest knockout rank. In contrast, as the up-regulation
#' would not affect the integrity of the network, you may select the features with highest (most significant) ExIR-based and
#' up-regulation-based ranks.
#' A shiny app has also been developed for Running the ExIR model, visualization of its results as well as computational
#' simulation of knockout and/or up-regulation of its top candidate outputs, which is accessible using
#' the `influential::runShinyApp("ExIR")` command.
#' You can also access the shiny app online at https://influential.erc.monash.edu/.
#' @param exir_output The output of the ExIR model (optional).
#' @param graph A graph (network) of the igraph class (not required if the exir_output is inputted).
#' @param ko_vertices A vector of desired vertices/features to knockout. Default is set to V(graph) meaning to assess
#' the knockout of all vertices/features.
#' @param upregulate_vertices A vector of desired vertices/features to up-regulate. Default is set to V(graph) meaning to assess
#' the up-regulation of all vertices/features.
#' @param beta Non-negative scalar corresponding to the SIRIR model. The rate of infection of an individual that is susceptible
#' and has a single infected neighbor. The infection rate of a susceptible individual with n
#' infected neighbors is n times beta. Formally this is the rate parameter of an exponential
#' distribution.
#' @param gamma Positive scalar corresponding to the SIRIR model. The rate of recovery of an infected individual.
#' Formally, this is the rate parameter of an exponential distribution.
#' @param no.sim Integer scalar corresponding to the SIRIR model. The number of simulation runs to perform SIR model on for the
#' original network as well perturbed networks generated by leave-one-out technique.
#' You may choose a different no.sim based on the available memory on your system.
#' @param node_verbose Logical; whether the process of Parallel Socket Cluster creation should be printed (default is FALSE).
#' @param loop_verbose Logical; whether the accomplishment of the evaluation of network nodes in each loop should be printed (default is TRUE).
#' @param ncores Integer; the number of cores to be used for parallel processing. If ncores == "default" (default), the number of 
#' cores to be used will be the max(number of available cores) - 1. We recommend leaving ncores argument as is (ncores = "default").
#' @param seed A single value, interpreted as an integer to be used for random number generation.
#' @return Depending on the input data, a list including one to three data frames of knockout/up-regulation rankings.
#' @keywords comp_manipulate
#' @family integrative ranking functions
#' @seealso \code{\link[influential]{exir}}, \code{\link[influential]{sirir}},
#' and \code{\link[igraph]{sir}} for a complete description on SIR model
#' @export comp_manipulate
#' @examples
#' \dontrun{
#' set.seed(1234)
#' My_graph <- igraph::sample_gnp(n=50, p=0.05)
#' GraphVertices <- V(My_graph)
#' Computational_manipulation <- comp_manipulate(graph = My_graph, beta = 0.5,
#'                                               gamma = 1, no.sim = 10, seed = 1234)
#'                                               }
#' @importFrom igraph vcount as_ids sir
#' @importFrom foreach %dopar% foreach

comp_manipulate <- function(exir_output = NULL,
                            graph = NULL,
                            ko_vertices = igraph::V(graph),
                            upregulate_vertices = igraph::V(graph),
                            beta = 0.5,
                            gamma = 1,
                            no.sim = 100,
                            node_verbose = FALSE,
                            loop_verbose = TRUE,
                            ncores = "default", 
                            seed = 1234) {
  
  ##**************************##
  # Take care of input graph
  if(!is.null(exir_output) && inherits(exir_output, "ExIR_Result")) {
    graph <- exir_output$Graph
  }
  
  ##**************************##
  # Over-expression function
  
  overexpr <- function(graph, vertices = upregulate_vertices, beta = beta, gamma = gamma,
                       no.sim = no.sim, loop_verbose = loop_verbose, 
                       ncores = ncores,
                       node_verbose = node_verbose, seed = seed) {
    
    suppressWarnings({
      
      # Make clusters for parallel processing
      cl <- parallel::makeCluster(ifelse(ncores == "default", parallel::detectCores() - 1, ncores), 
                                  outfile=ifelse(node_verbose, "", "NULL"))
      doParallel::registerDoParallel(cl)
      
      temp.loocr.table <- data.frame(difference.value = rep(NA, length(vertices)), rank = rep(NA, length(vertices)))
      if (inherits(vertices, "character")) {
        rownames(temp.loocr.table) <- vertices
      } else if (inherits(vertices, "igraph.vs")) {
        rownames(temp.loocr.table) <- igraph::as_ids(vertices)
      }
      
      set.seed(seed)
      all.included.spread <- igraph::sir(graph = graph, beta = beta,
                                         gamma = gamma, no.sim = no.sim)
      all.mean_spread <- foreach(i = 1:length(all.included.spread), .combine = "c", .verbose = loop_verbose) %dopar% {
        max(all.included.spread[[i]]$NR)
      }
      all.mean_spread <- mean(all.mean_spread)
      
      results <- foreach(s = 1:length(vertices), .combine = "c", .verbose = loop_verbose) %dopar% {
        all.edges <- as.data.frame(igraph::as_edgelist(graph))
        all.desired.edges.index <- unlist(apply(X = all.edges, MARGIN = 2,
                                                FUN = function(i) which(i %in% rownames(temp.loocr.table)[s])))
        
        all.edges <- all.edges[c(all.desired.edges.index),]
        
        all.edges.from.indices <- unlist(lapply(X = all.edges[,1],
                                                FUN = function(i) which(igraph::as_ids(V(graph)) %in% i)))
        
        all.edges.to.indices <- unlist(lapply(X = all.edges[,2],
                                              FUN = function(i) which(igraph::as_ids(V(graph)) %in% i)))
        
        all.edges.desired.indices <- data.frame(from = all.edges.from.indices, to = all.edges.to.indices)
        
        temp.graph <- igraph::add_vertices(graph = graph, nv = 1,
                                           name = paste0(rownames(temp.loocr.table)[s], "Fold1"))
        
        desired.vertex.index <- match(rownames(temp.loocr.table)[s], igraph::as_ids(V(temp.graph)))
        
        desired.vertexFold1.index <- match(paste0(rownames(temp.loocr.table)[s], "Fold1"),
                                           igraph::as_ids(V(temp.graph))
        )
        
        all.edges.desired.indices[which(all.edges.desired.indices[,1] == desired.vertex.index),1] <- desired.vertexFold1.index
        all.edges.desired.indices[which(all.edges.desired.indices[,2] == desired.vertex.index),2] <- desired.vertexFold1.index
        all.edges.desired.indices <- apply(X = all.edges.desired.indices, MARGIN = 1,
                                           FUN = function(i) paste(i[1], i[2], sep = ","))
        
        all.edges.desired.indices <- paste(all.edges.desired.indices, collapse = ",")
        
        all.edges.desired.indices <- as.numeric(unlist(strsplit(x = all.edges.desired.indices, split = ",")))
        
        temp.graph <- igraph::add_edges(graph = temp.graph, edges = all.edges.desired.indices)
        
        set.seed(seed)
        loocr.spread <- igraph::sir(graph = temp.graph, beta = beta,
                                    gamma = gamma, no.sim = no.sim)
        loocr.mean_spread <- foreach(h = 1:length(loocr.spread), .combine = "c", .verbose = loop_verbose) %dopar% {
          max(loocr.spread[[h]]$NR)
        }
        loocr.mean_spread <- mean(loocr.mean_spread)
        
        loocr.mean_spread - all.mean_spread
      }
      
      temp.loocr.table$difference.value <- results
      
      # Stop the parallel backend
      parallel::stopCluster(cl)
      
      temp.loocr.table$rank <- rank(-temp.loocr.table$difference.value, ties.method = "min")
      
    })
    
    return(temp.loocr.table)
  }
  
  ##**************************##
  # Knockout results
  if(!is.null(ko_vertices)) {
    base::suppressWarnings(
      ko_results <- influential::sirir(
        graph = graph,
        vertices = ko_vertices,
        beta = beta,
        gamma = gamma,
        no.sim = no.sim,
        loop_verbose = loop_verbose, 
        node_verbose = node_verbose,
        ncores = ncores,
        seed = seed
      )
    )
    
    ko_results <- cbind(Feature_name = rownames(ko_results),
                        ko_results,
                        Manipulation_type = "Knockout")
    
    rownames(ko_results) <- NULL
    colnames(ko_results)[3] <- "Rank"
    
    # correct the orders
    ko_results <- ko_results[order(ko_results[,3]),]
    
  } else {ko_results <- NULL}
  
  ##**************************##
  
  # Over-expression results
  if(!is.null(upregulate_vertices)) {
    base::suppressWarnings(
      overexpr_results <- overexpr(
        graph = graph,
        vertices = upregulate_vertices,
        beta = beta,
        gamma = gamma,
        no.sim = no.sim,
        loop_verbose = loop_verbose, 
        node_verbose = node_verbose,
        ncores = ncores,
        seed = seed
      )
    )
    
    overexpr_results <- cbind(Feature_name = rownames(overexpr_results),
                              overexpr_results,
                              Manipulation_type = "Up-regulation")
    
    rownames(overexpr_results) <- NULL
    colnames(overexpr_results)[3] <- "Rank"
    
    # correct the orders
    overexpr_results <- overexpr_results[order(overexpr_results[,3]),]
    
  } else {overexpr_results <- NULL}
  
  ##**************************##
  
  # Combined results
  if(!is.null(ko_results) & !is.null(overexpr_results)) {
    combined_results <- rbind(ko_results,
                              overexpr_results)
    
    # Correct the combined ranks
    combined_results[,3] <- rank(-combined_results[,2],
                                 ties.method = "min")
    
    # correct the orders
    combined_results <- combined_results[order(combined_results[,3]),]
    
  } else {
    combined_results <- NULL
  }
  
  ##**************************##
  
  # Correct results data frames
  
  ko_results <- ko_results[,-2]
  overexpr_results <- overexpr_results[,-2]
  combined_results <- combined_results[,-2]
  
  ##**************************##
  
  # Final results
  final.results <- list(Knockout = ko_results,
                        Up_regulation = overexpr_results,
                        Combined = combined_results
  )
  if(sum(sapply(final.results, is.null)) == 3) {
    cat("You should input the name of at least a vertix (feature/gene/etc)
        in the 'ko_vertices' or 'upregulate_vertices' argument.")
  } else {
    non_null_index <- which(sapply(final.results, function(i) (!is.null(i))))
    final.results <- final.results[non_null_index]
    names(final.results) <- names(non_null_index)
    return(final.results)
  }
}