R/mle.getPermutations.r

In lashmore/MolEndoMatch: Molecular Endophenotype Patient Matchmaking

#' Generate the "adaptive walk" node permutations, starting from a given perturbed variable
#'
#' This function calculates the node permutation starting from a given perturbed variable in a subset of variables in the background knowledge graph.
#' @param sig.nodes - The names (as a character vector) of the variables that were significantly perturbed in the given patient.
#' @param patient - The row number (integer) in data.pvals associated with the patient being processed.
#' @param ig.pt - An igraph object of the background knowledge graph.
#' @keywords probability diffusion algorithm
#' @keywords network modularity
#' @keywords weighted Bonferroni correction
#' @export mle.getPermutations
#' @examples
#' mle.getPermutations(sig.nodes, patient, ig.pt)
mle.getPermutations <- function(sig.nodes, patient, ig.pt) {
  # Phase 1: Do adaptive permutations ahead of time for each possible startNode in subGraphS.
  # Get all node names, so we know what all possible startNodes are.
  permutationByStartNode <- list()
  for (n in 1:length(sig.nodes)) {
    #print(sprintf("Calculating permutation %d of %d for patient %d...", n, length(sig.nodes), patient))
    # Draw all nodes in graph
    current_node_set <- NULL
    stopIterating=FALSE;
    startNode <- sig.nodes[n]
    hits <- startNode
    numMisses <- 0
    currentGraph <- igraphObjectG
    while (stopIterating==FALSE) {
      current_node_set <- c(current_node_set, startNode)
      #print(sprintf("Draws=%d, Hits=%d/%d", length(current_node_set), length(hits), length(sig.nodes)))
      #STEP 4: Diffuse p1 to connected nodes from current draw, startNode node
      sumHits <- as.vector(matrix(0, nrow=length(V(ig.pt)$name), ncol=1))
      names(sumHits) <- V(ig.pt)$name
      for (hit in 1:length(hits)) {
        #For unseen nodes, clear probabilities and add probability (p0/#unseen nodes)
        for (t in 1:length(currentGraph)) {
          currentGraph[[t]] = 0; #set probabilities of all nodes to 0
        }
        #determine base p0 probability
        baseP = p0/(length(currentGraph)-length(current_node_set));
        for (t in 1:length(currentGraph)) {
          if (!(names(currentGraph[t]) %in% current_node_set)) {
            currentGraph[[t]] = baseP;  #set probabilities of unseen nodes to diffused p0 value, baseP
          } else {
            currentGraph[[t]] = 0;
          }
        }
        p0_event <- sum(unlist(currentGraph[!(names(currentGraph) %in% current_node_set)]))
        currentGraph <- graph.diffuseP1(p1, hits[hit], currentGraph, currentGraph[current_node_set], 1, verbose=FALSE)
        p1_event <- sum(unlist(currentGraph[!(names(currentGraph) %in% current_node_set)]))
        if (abs(p1_event-1)>thresholdDiff) {
          extra.prob.to.diffuse <- 1-p1_event
          currentGraph[names(current_node_set)] <- 0
          currentGraph[!(names(currentGraph) %in% names(current_node_set))] <- unlist(currentGraph[!(names(currentGraph) %in% names(current_node_set))]) + extra.prob.to.diffuse/sum(!(names(currentGraph) %in% names(current_node_set)))
        }
        sumHits <- sumHits + unlist(currentGraph)
      }
      sumHits <- sumHits/length(hits)

      #Set startNode to a node that is the max probability in the new currentGraph
      # When there are ties, choose amongst winners for highest degree.
      maxProb <- names(which.max(sumHits))
      # Break ties: TODO. Just pick top of alphabet in maxProb vector.
      startNode <- names(currentGraph[maxProb[1]])
      if (startNode %in% sig.nodes) {
        hits <- c(hits, startNode)
        numMisses <- 0
      } else {
        numMisses<-numMisses+1
      }

      if (all(sig.nodes %in% hits) || numMisses>thresholdDrawT || length(c(startNode,current_node_set))>=(length(V(ig.pt)$name))) {
        current_node_set <- c(current_node_set, startNode)
        stopIterating = TRUE
      }
    }
    permutationByStartNode[[n]] <- current_node_set
  }
  names(permutationByStartNode) <- sig.nodes
  return(permutationByStartNode)
}