Nothing
##' @name make_state
##' @aliases make_state_matrix
##' @rdname make_state
##'
##' @title Make State Matrix
##'
##' @description Functions to compute the matrix of states (1 for activating and -1 for inhibiting)
##' for link signed correlations, from a vector of edge states to a signed adjacency matrix for use
##' in \code{\link[graphsim]{generate_expression}}. This resolves edge states to determine the sign
##' of all correlations between nodes in a network. These are computed interally for sigma matrices
##' as required.
##' @param graph An \code{\link[igraph:aaa-igraph-package]{igraph}} object. May be directed or weighted as long as
##' a shortest path can be computed.
##' @param state numeric vector. Vector of length E(graph). Sign used to calculate state matrix,
##' may be an integer state or inferred directly from expected correlations for each edge. May be
##' applied a scalar across all edges or as a vector for each edge respectively. May also be
##' entered as text for "activating" or "inhibiting" or as integers for activating (0,1) or
##' inhibiting (-1,2). Compatible with inputs for \code{\link[graphsim]{plot_directed}}. Vector
##' input is supported either directly calling the function with a value for each edge in
##' \code{E(graph)} or as an edge "attribute" in the igraph object (using
##' \code{E(g)$state <- states}).
##' @keywords graph network igraph mvtnorm simulation
##' @importFrom igraph graph.edgelist get.edge.attribute get.edgelist
##' @import igraph
##'
##' @family graphsim functions
##' @family generate simulated expression functions
##' @seealso
##' See also \code{\link[graphsim]{generate_expression}} for computing the simulated data,
##' \code{\link[graphsim]{make_sigma}} for computing the Sigma (\eqn{\Sigma}) matrix,
##' and
##' \code{\link[graphsim]{make_distance}} for computing distance from a graph object.
##'
##' See also \code{\link[graphsim]{plot_directed}} for plotting graphs or
##' \code{\link[gplots]{heatmap.2}} for plotting matrices.
##'
##' See also \code{\link[graphsim]{make_laplacian}}, \code{\link[graphsim]{make_commonlink}},
##' or \code{\link[graphsim]{make_adjmatrix}} for computing input matrices.
##'
##' See also \code{\link[igraph:aaa-igraph-package]{igraph}} for handling graph objects.
##'
##' @author Tom Kelly \email{tom.kelly@@riken.jp}
##'
##' @examples
##'
##' # construct a synthetic graph module
##' library("igraph")
##' graph_test_edges <- rbind(c("A", "B"), c("B", "C"), c("B", "D"))
##' graph_test <- graph.edgelist(graph_test_edges, directed = TRUE)
##'
##' # compute state matrix for toy example
##' state_matrix <- make_state_matrix(graph_test)
##'
##' # construct a synthetic graph network
##' graph_structure_edges <- rbind(c("A", "C"), c("B", "C"), c("C", "D"), c("D", "E"),
##' c("D", "F"), c("F", "G"), c("F", "I"), c("H", "I"))
##' graph_structure <- graph.edgelist(graph_structure_edges, directed = TRUE)
##'
##' # compute state matrix for toy network
##' graph_structure_state_matrix <- make_state_matrix(graph_structure)
##' graph_structure_state_matrix
##'
##' # compute state matrix for toy network with inhibitions
##' edge_state <- c(1, 1, -1, 1, 1, 1, 1, -1)
##' # edge states are a variable
##' graph_structure_state_matrix <- make_state_matrix(graph_structure, state = edge_state)
##' graph_structure_state_matrix
##'
##' # compute state matrix for toy network with inhibitions
##' E(graph_structure)$state <- c(1, 1, -1, 1, 1, 1, 1, -1)
##' # edge states are a graph attribute
##' graph_structure_state_matrix <- make_state_matrix(graph_structure)
##' graph_structure_state_matrix
##'
##' library("igraph")
##' graph_test_edges <- rbind(c("A", "B"), c("B", "C"), c("B", "D"))
##' graph_test <- graph.edgelist(graph_test_edges, directed = TRUE)
##' state_matrix <- make_state_matrix(graph_test)
##'
##' # import graph from package for reactome pathway
##' # TGF-\eqn{\Beta} receptor signaling activates SMADs (R-HSA-2173789)
##' TGFBeta_Smad_graph <- identity(TGFBeta_Smad_graph)
##'
##' # compute sigma (\eqn{\Sigma}) matrix from geometric distance directly from TGF-\eqn{\Beta} pathway
##' TFGBeta_Smad_state <- E(TGFBeta_Smad_graph)$state
##' table(TFGBeta_Smad_state)
##' # states are edge attributes
##' state_matrix_TFGBeta_Smad <- make_state_matrix(TGFBeta_Smad_graph)
##' # visualise matrix
##' library("gplots")
##' heatmap.2(state_matrix_TFGBeta_Smad , scale = "none", trace = "none",
##' dendrogram = "none", Rowv = FALSE, Colv = FALSE,
##' col = colorpanel(50, "blue", "white", "red"))
##'
##' # compare the states to the sign of expected correlations in the sigma matrix
##' sigma_matrix_TFGBeta_Smad_inhib <- make_sigma_mat_dist_graph(TGFBeta_Smad_graph,
##' cor = 0.8,
##' absolute = FALSE)
##' # visualise matrix
##' heatmap.2(sigma_matrix_TFGBeta_Smad_inhib,
##' scale = "none", trace = "none",
##' dendrogram = "none", Rowv = FALSE, Colv = FALSE,
##' col = colorpanel(50, "blue", "white", "red"))
##'
##' # compare the states to the sign of final correlations in the simulated matrix
##' TFGBeta_Smad_data <- generate_expression(100, TGFBeta_Smad_graph, cor = 0.8)
##' heatmap.2(cor(t(TFGBeta_Smad_data)), scale = "none", trace = "none",
##' dendrogram = "none", Rowv = FALSE, Colv = FALSE,
##' col = colorpanel(50, "blue", "white", "red"))
##'
##'
##' @return An integer matrix indicating the resolved state
##' (activating or inhibiting for each edge or path between nodes)
##' @export
make_state_matrix <- function(graph, state = NULL){
if(!is.null(get.edge.attribute(graph, "state"))){
state <- get.edge.attribute(graph, "state")
} else {
# add default state if not specified
if(is.null(state)){
state <- "activating"
}
}
if(is.null(get.edge.attribute(graph, "state"))){
E(graph)$state <- state
}
# this could also be done with igraph::simplify(graph, remove.multiple = TRUE, remove.loops = TRUE, edge.attr.comb = igraph_opt("edge.attr.comb"))
## to do: migrate to this
graph <- as.undirected(graph, mode = "collapse", edge.attr.comb = function(x) ifelse(any(x %in% list(-1, 2, "inhibiting", "inhibition")), -1, 1))
# remove duplicate edges (and corresponding states)
graph <- as.directed(graph, mode = "arbitrary")
state <- get.edge.attribute(graph, "state")
if(length(state) == 1) state <- rep(state, length(E(graph)))
state[state == 2] <- -1
state[state == 1 | state == 0] <- 1
if(is.numeric(state)) state[!(state %in% -1:2)] <- sign(state[!(state %in% -1:2)])
if(is.character(state)){
state <- as.list(state)
state[grep("activating", state)] <- 1
state[grep("activation", state)] <- 1
state[grep("activate", state)] <- 1
state[grep("active", state)] <- 1
state[grep("positive", state)] <- 1
state[grep("inhibiting", state)] <- -1
state[grep("inhibition", state)] <- -1
state[grep("inhibitory", state)] <- -1
state[grep("inhibit", state)] <- -1
state[grep("negative", state)] <- -1
if(is.character(state)){
warning("Please give state as a scalar or vector of length(E(graph)): input must be 'activating', 'inhibiting' or an integer")
}
state <- unlist(as.numeric(state))
}
if(!all(state %in% -1:2)){
state <- sign(state) # coerce to vector or 1 and -1 if not already
warning("State inferred from non-integer weighted edges: Please give numeric states as integers: 0 or 1 for activating, -1 or 2 for inhibiting")
}
#define starting states
state_mat <- matrix(1, length(V(graph)), length(V(graph)))
rownames(state_mat) <- colnames(state_mat) <- names(V(graph))
# check for inhibiting edges
if(all(state == 1)){
#return without resolving for activations
return(state_mat)
} else if (length(V(graph)) == 1){
state_mat <- 1
return(state_mat)
} else if (length(V(graph)) == 2){
if(length(state) == 1){
state_mat[row(state_mat) != col(state_mat)] <- state # one edge
return(state_mat)
} else{
warning("only one edge state expected")
}
} else {
#resolve inhibitions downstream
# define shortest paths to all possible nodes
#check if connected or use connected subgraphs
compute_paths <- function(graph){
if(is.connected(graph)){
#define paths from 1st node
paths <- shortest_paths(as.undirected(graph), V(graph)$name[1])$vpath
} else {
subgraphs <- decompose(graph)
nodes <- sapply(subgraphs, function(subgraph) V(subgraph)$name[1])
paths <- as.list(rep(NA, length(V(graph))))
jj <- 0
for(ii in 1:length(nodes)){
subpaths <- shortest_paths(as.undirected(subgraphs[[ii]]), V(subgraphs[[ii]])$name[1])$vpath
paths[jj+1:length(subpaths)] <- subpaths
jj <- jj + length(subpaths)
}
}
paths
}
paths <- compute_paths(graph)
#check for cycles
is.cyclic <- function(paths){
cylic_paths <- sapply(paths, function(path){
if(length(path) <= 1){
FALSE
} else if (path[1] == path[length(path)]){
TRUE
} else {
FALSE
}
})
any(cylic_paths)
}
if(is.cyclic(paths)){
warning("Graph contains a cycle, computing minimal spanning tree. This may result in unresolved inhibitions.")
tree <- minimum.spanning.tree(graph)
paths <- compute_paths(tree)
if(is.cyclic(paths)){
warning("Graph contains cycles and cannot compute a minimal spanning tree. This will result in unresolved inhibitions.")
stop()
}
}
#sort paths into longest to shortest
paths <- paths[order(sapply(paths, length), decreasing = TRUE)]
#remove subpaths
for(ii in 2:length(paths)){
remove <- FALSE
for(jj in 1:(ii-1)){
# test if all nodes in path in a larger path above
if(all(names(paths[[ii]]) %in% names(paths[[jj]]))){
remove <- TRUE
}
}
if(remove){
paths[[ii]] <- NA
}
}
na.omit.list <- function(y) { return(y[!sapply(y, function(x) all(is.na(x)))]) }
paths <- na.omit.list(paths)
# check inhibition state through graph paths
edges <- as.matrix(get.edgelist(graph))
for(ii in 1:length(paths)){
if(length(paths[[ii]]) > 1){
state_path <- c(1, rep(NA, length(paths[[ii]])-1))
for(jj in 2:length(paths[[ii]])){
#find edges
kk <- which(
edges[,1] == names(paths[[ii]])[jj-1] & edges[,2] == names(paths[[ii]])[jj] |
edges[,2] == names(paths[[ii]])[jj-1] & edges[,1] == names(paths[[ii]])[jj]
)
#check for multiple edges
if(length(kk) > 1){
#check for same edges
if(!all(apply(edges[kk,], 2, function(x) all(x == x[1])))){
#reverse edge if no match
edges[kk,][edges[kk[1],1] != edges[kk,1],] <- edges[kk,][edges[kk[1],1] != edges[kk,1],2:1]
}
#check is still no match
if(all(apply(edges[kk,], 2, function(x) all(x == x[1])))){
#check for same state
if(all(state[kk] == state[kk][1])){
state[kk] <- state[kk][1]
kk <- kk[1]
} else {
warning(paste("Conflicting state for redundant edges", edges[kk[1],1], edges[kk[1],2], "with state", state[kk]))
Mode <- function(x) unique(x)[which.max(tabulate(match(x, unique(x))))][1]
state[kk] <- Mode(state[kk])
kk <- kk[1]
}
}
}
#state for edges
if(length(state_path[jj]) == length(state[kk])){
state_path[jj] <- state[kk]
} else {
print(ii)
print(jj)
stop("State not compatible with graph paths")
}
}
#cumulative product changes state along a path (*-1)
state_path <- cumprod(state_path)
#skip if activating path
if(any(state_path != 1)){
#cross-product produces a matrix form
state_cross <- state_path %*% t(state_path)
#add to state_matrix (all adjusted to start of path)
state_mat[names(paths[[ii]]), names(paths[[ii]])] <- state_cross
}
}
}
#check for negative and positive links not in paths (against the 1st node)
neg_nodes <- colnames(state_mat)[which(state_mat[1,] == -1)]
pos_nodes <- colnames(state_mat)[2:ncol(state_mat)][which(state_mat[1,2:ncol(state_mat)] == 1)] # can assume at least 3 nodes due to checks above
#check for if not computed in shorted paths already
##paths <- compute_paths(tree) # not minimal spanning tree
for(aa in neg_nodes){
for(bb in pos_nodes){
aa_in_path <- sapply(paths, function(path) aa %in% names(path))
bb_in_path <- sapply(paths, function(path) bb %in% names(path))
#check if in ANY same path
if(!any(aa_in_path & bb_in_path)){
#if not invert
state_mat[aa, bb] <- state_mat[bb, aa] <- -1
}
}
}
# check if same sign (to do: allow disconnected)
for(aa in pos_nodes){
for(bb in pos_nodes){
aa_in_path <- sapply(paths, function(path) aa %in% names(path))
bb_in_path <- sapply(paths, function(path) bb %in% names(path))
#check if in ANY same path
if(!any(aa_in_path & bb_in_path)){
#if not invert
state_mat[aa, bb] <- state_mat[bb, aa] <- 1
}
}
}
#ensure symmetric matrix
state_mat[t(state_mat) == -1] <- -1
return(state_mat)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.