--- a +++ b/R/RandomWalkRestartMH_functions.R @@ -0,0 +1,381 @@ +# 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") +}