R/RandomWalkRestartMH_functions.R

In abodein/netOmics: Multi-Omics (time-course) network-based integration and interpretation

# Credit to Valdeolivas et al.
# All the functions below come from the package RandomWalkRestartMH. 
# Due to build errors in Bioconducteur, and to avoid depreciation of the netOmics package, all functions imported by the package have been repatriated here.

# need ti keep:
# - RandomWalkRestartMH::create.multiplex
# - RandomWalkRestartMH::compute.adjacency.matrix
# - RandomWalkRestartMH::normalize.multiplex.adjacency
# - RandomWalkRestartMH::Random.Walk.Restart.Multiplex

#' @importFrom igraph set_vertex_attr
create.multiplex <- function(LayersList,...){
  
  if (!class(LayersList) == "list"){
    stop("The input object should be a list of graphs.")
  }
  
  
  Number_of_Layers <- length(LayersList)
  SeqLayers <- seq(Number_of_Layers)
  Layers_Name <- names(LayersList)
  
  # if (!all(sapply(SeqLayers, function(x) is.igraph(LayersList[[x]])))){
  #   stop("Not igraph objects")
  # }
  
  Layer_List <- lapply(SeqLayers, function (x) {
    if (is.null(V(LayersList[[x]])$name)){
      LayersList[[x]] <- 
        igraph::set_vertex_attr(LayersList[[x]],"name", 
                        value=seq(1,vcount(LayersList[[x]]),by=1))
    } else {
      LayersList[[x]]
    }
  })
  
  ## We simplify the layers 
  Layer_List <- 
    lapply(SeqLayers, function(x) simplify.layers(Layer_List[[x]]))
  
  ## We set the names of the layers. 
  
  if (is.null(Layers_Name)){
    names(Layer_List) <- paste0("Layer_", SeqLayers)
  } else {
    names(Layer_List) <- Layers_Name
  }
  
  ## We get a pool of nodes (Nodes in any of the layers.)
  Pool_of_Nodes <- 
    sort(unique(unlist(lapply(SeqLayers, 
                              function(x) V(Layer_List[[x]])$name))))
  
  Number_of_Nodes <- length(Pool_of_Nodes)
  
  Layer_List <-
    lapply(Layer_List, add.missing.nodes,Number_of_Layers,Pool_of_Nodes)
  
  # We set the attributes of the layer
  counter <- 0 
  Layer_List <- lapply(Layer_List, function(x) { 
    counter <<- counter + 1; 
    igraph::set_edge_attr(x,"type",igraph::E(x), value = names(Layer_List)[counter])
  })
  
  
  MultiplexObject <- c(Layer_List,list(Pool_of_Nodes=Pool_of_Nodes,
                                       Number_of_Nodes_Multiplex=Number_of_Nodes, 
                                       Number_of_Layers=Number_of_Layers))
  
  class(MultiplexObject) <- "Multiplex"
  
  return(MultiplexObject)
}



# internal 
#' @importFrom igraph as.undirected is_weighted E simplify
simplify.layers <- function(Input_Layer){
  
  ## Undirected Graphs
  Layer <- igraph::as.undirected(Input_Layer, mode = c("collapse"),
                         edge.attr.comb = igraph::igraph_opt("edge.attr.comb"))
  
  ## Unweighted or Weigthed Graphs
  if (igraph::is_weighted(Layer)){
    b <- 1
    weigths_layer <- igraph::E(Layer)$weight
    if (min(weigths_layer) != max(weigths_layer)){
      a <- min(weigths_layer)/max(weigths_layer)
      range01 <- (b-a)*(weigths_layer-min(weigths_layer))/
        (max(weigths_layer)-min(weigths_layer)) + a
      igraph::E(Layer)$weight <- range01
    } else {
      igraph::E(Layer)$weight <- rep(1, length(weigths_layer))
    }
  } else {
    igraph::E(Layer)$weight <- rep(1, ecount(Layer))
  }
  
  ## Simple Graphs
  Layer <- 
    igraph::simplify(Layer,remove.multiple = TRUE,remove.loops = TRUE, 
                     edge.attr.comb=mean)
  
  return(Layer)
}

#' @importFrom igraph add_vertices
add.missing.nodes <- function (Layers,Nr_Layers,NodeNames) {
  
  igraph::add_vertices(Layers,
               length(NodeNames[which(!NodeNames %in% igraph::V(Layers)$name)]),
               name=NodeNames[which(!NodeNames %in%  igraph::V(Layers)$name)])
}


#' @importFrom Matrix Diagonal bdiag
#' @importFrom igraph as_adjacency_matrix is_weighted

compute.adjacency.matrix <- function(x,delta = 0.5)
{
  if (!isMultiplex(x) & !isMultiplexHet(x)) {
    stop("Not a Multiplex or Multiplex Heterogeneous object")
  }
  if (delta > 1 || delta <= 0) {
    stop("Delta should be between 0 and 1")
  }
  
  N <- x$Number_of_Nodes_Multiplex
  L <- x$Number_of_Layers
  
  ## We impose delta=0 in the monoplex case.
  if (L==1){
    delta = 0
  }
  
  Layers_Names <- names(x)[seq(L)]
  
  ## IDEM_MATRIX.
  Idem_Matrix <- Matrix::Diagonal(N, x = 1)
  
  counter <- 0 
  Layers_List <- lapply(x[Layers_Names],function(x){
    
    counter <<- counter + 1;    
    if (igraph::is_weighted(x)){ 
      Adjacency_Layer <-  igraph::as_adjacency_matrix(x,sparse = TRUE, 
                                              attr = "weight")
    } else {
      Adjacency_Layer <-  igraph::as_adjacency_matrix(x,sparse = TRUE)
    }
    
    Adjacency_Layer <- Adjacency_Layer[order(rownames(Adjacency_Layer)),
                                       order(colnames(Adjacency_Layer))]
    colnames(Adjacency_Layer) <- 
      paste0(colnames(Adjacency_Layer),"_",counter)
    rownames(Adjacency_Layer) <- 
      paste0(rownames(Adjacency_Layer),"_",counter)
    Adjacency_Layer
  })
  
  MyColNames <- unlist(lapply(Layers_List, function (x) unlist(colnames(x))))
  MyRowNames <- unlist(lapply(Layers_List, function (x) unlist(rownames(x))))
  names(MyColNames) <- c()
  names(MyRowNames) <- c()
  SupraAdjacencyMatrix <- (1-delta)*(Matrix::bdiag(unlist(Layers_List)))
  colnames(SupraAdjacencyMatrix) <-MyColNames
  rownames(SupraAdjacencyMatrix) <-MyRowNames
  
  offdiag <- (delta/(L-1))*Idem_Matrix
  
  i <- seq_len(L)
  Position_ini_row <- 1 + (i-1)*N
  Position_end_row <- N + (i-1)*N
  j <- seq_len(L)
  Position_ini_col <- 1 + (j-1)*N
  Position_end_col <- N + (j-1)*N
  
  for (i in seq_len(L)){
    for (j in seq_len(L)){
      if (j != i){
        SupraAdjacencyMatrix[(Position_ini_row[i]:Position_end_row[i]),
                             (Position_ini_col[j]:Position_end_col[j])] <- offdiag
      }    
    }
  }
  
  SupraAdjacencyMatrix <- as(SupraAdjacencyMatrix, "dgCMatrix")
  return(SupraAdjacencyMatrix)
}
#' @importFrom Matrix t colSums
normalize.multiplex.adjacency <- function(x)
{
  if (!is(x,"dgCMatrix")){
    stop("Not a dgCMatrix object of Matrix package")
  }
  
  Adj_Matrix_Norm <- Matrix::t(Matrix::t(x)/(Matrix::colSums(x, na.rm = FALSE, dims = 1,
                                             sparseResult = FALSE)))
  
  return(Adj_Matrix_Norm)
}

Random.Walk.Restart.Multiplex <- function(x, MultiplexObject, Seeds, 
                                                  r=0.7,tau,MeanType="Geometric", DispResults="TopScores",...){

  
  L <- MultiplexObject$Number_of_Layers
  N <- MultiplexObject$Number_of_Nodes
  
  Seeds <- as.character(Seeds)
  if (length(Seeds) < 1 | length(Seeds) >= N){
    stop("The length of the vector containing the seed nodes is not 
         correct")
  } else {
    if (!all(Seeds %in% MultiplexObject$Pool_of_Nodes)){
      stop("Some of the seeds are not nodes of the network")
      
    }
  }
  
  if (r >= 1 || r <= 0) {
    stop("Restart partameter should be between 0 and 1")
  }
  
  if(missing(tau)){
    tau <- rep(1,L)/L
  } else {
    tau <- as.numeric(tau)
    if (sum(tau)/L != 1) {
      stop("The sum of the components of tau divided by the number of 
           layers should be 1")
    }
  }
  
  if(!(MeanType %in% c("Geometric","Arithmetic","Sum"))){
    stop("The type mean should be Geometric, Arithmetic or Sum")
  }
  
  if(!(DispResults %in% c("TopScores","Alphabetic"))){
    stop("The way to display RWRM results should be TopScores or
         Alphabetic")
  }
  
  ## We define the threshold and the number maximum of iterations for
  ## the random walker.
  Threeshold <- 1e-10
  NetworkSize <- ncol(x)
  
  ## We initialize the variables to control the flux in the RW algo.
  residue <- 1
  iter <- 1
  
  ## We compute the scores for the different seeds.
  Seeds_Score <- get.seed.scoresMultiplex(Seeds,L,tau)
  
  ## We define the prox_vector(The vector we will move after the first RWR
  ## iteration. We start from The seed. We have to take in account
  ## that the walker with restart in some of the Seed nodes, depending on
  ## the score we gave in that file).
  prox_vector <- matrix(0,nrow = NetworkSize,ncol=1)
  
  prox_vector[which(colnames(x) %in% Seeds_Score[,1])] <- (Seeds_Score[,2])
  
  prox_vector  <- prox_vector/sum(prox_vector)
  restart_vector <-  prox_vector
  
  while(residue >= Threeshold){
    
    old_prox_vector <- prox_vector
    prox_vector <- (1-r)*(x %*% prox_vector) + r*restart_vector
    residue <- sqrt(sum((prox_vector-old_prox_vector)^2))
    iter <- iter + 1;
  }
  
  NodeNames <- character(length = N)
  Score = numeric(length = N)
  
  rank_global <- data.frame(NodeNames = NodeNames, Score = Score)
  rank_global$NodeNames <- gsub("_1", "", row.names(prox_vector)[seq_len(N)])
  
  if (MeanType=="Geometric"){
    rank_global$Score <- geometric.mean(as.vector(prox_vector[,1]),L,N)    
  } else {
    if (MeanType=="Arithmetic") {
      rank_global$Score <- regular.mean(as.vector(prox_vector[,1]),L,N)    
    } else {
      rank_global$Score <- sumValues(as.vector(prox_vector[,1]),L,N)    
    }
  }
  
  if (DispResults=="TopScores"){
    ## We sort the nodes according to their score.
    Global_results <- 
      rank_global[with(rank_global, order(-Score, NodeNames)), ]
    
    ### We remove the seed nodes from the Ranking and we write the results.
    Global_results <- 
      Global_results[which(!Global_results$NodeNames %in% Seeds),]
  } else {
    Global_results <- rank_global    
  }
  
  rownames(Global_results) <- c()
  
  RWRM_ranking <- list(RWRM_Results = Global_results,Seed_Nodes = Seeds)
  
  class(RWRM_ranking) <- "RWRM_Results"
  return(RWRM_ranking)
}

#' @method print RWRM_Results
#' @export
print.RWRM_Results <- function(x,...)
{
  cat("Top 10 ranked Nodes:\n")
  print(head(x$RWRM_Results,10))
  cat("\nSeed Nodes used:\n")
  print(x$Seed_Nodes)
}

get.seed.scoresMultiplex <- function(Seeds,Number_Layers,tau) {
  
  Nr_Seeds <- length(Seeds)
  
  Seeds_Seeds_Scores <- rep(tau/Nr_Seeds,Nr_Seeds)
  Seed_Seeds_Layer_Labeled <- 
    paste0(rep(Seeds,Number_Layers),sep="_",rep(seq(Number_Layers), 
                                                length.out = Nr_Seeds*Number_Layers,each=Nr_Seeds))
  
  Seeds_Score <- data.frame(Seeds_ID = Seed_Seeds_Layer_Labeled,
                            Score = Seeds_Seeds_Scores, stringsAsFactors = FALSE)
  
  return(Seeds_Score)
}


geometric.mean <- function(Scores, L, N) {
  
  FinalScore <- numeric(length = N)
  
  for (i in seq_len(N)){
    FinalScore[i] <- prod(Scores[seq(from = i, to = N*L, by=N)])^(1/L)
  }
  
  return(FinalScore)
}

regular.mean <- function(Scores, L, N) {
  
  FinalScore <- numeric(length = N)
  
  for (i in seq_len(N)){
    FinalScore[i] <- mean(Scores[seq(from = i, to = N*L, by=N)])
  }
  
  return(FinalScore)
}

sumValues <- function(Scores, L, N) {
  
  FinalScore <- numeric(length = N)
  
  for (i in seq_len(N)){
    FinalScore[i] <- sum(Scores[seq(from = i, to = N*L, by=N)])
  }
  
  return(FinalScore)
}

isMultiplex <- function (x) 
{
  is(x, "Multiplex")
}

isMultiplexHet <- function (x) 
{
  is(x, "MultiplexHet")
}