#' Edge Density Threshold
#'
#'@description
#'Get the adjacency matrix considering only the highest absolute value weights
#'up to the desired link density.
#'
#'@param weights correlation matrix
#'@param edge.density threshold value
#'
#'@export
density_link_threshold <- function(weights, edge.density=.05){
#Checks
if(nrow(weights)!=ncol(weights)) stop("matrix must be square")
if(any(rownames(weights)!=colnames(weights))) stop("Weights rownames and colnames must be the same")
if(edge.density<0 || edge.density>1) stop("edge.density must be in range [0,1]")
x <- weights
diag(x) <- 0
# Take correlation upper triangular values as vector
xv <- x[upper.tri(x)]
# Sort Aboslute values
xv.sort <- sort(abs(xv),decreasing=TRUE)
# Take the Absolute Min Values
th.value <- xv.sort[round(length(xv.sort)*edge.density)]
adj <- x * (abs(x)>=th.value)
return(adj)
}
#------------------------------------------------------------------------------#
#' Signed-Weighted Graph Layout
#'
#' @description It elaborates the coordinates for the representation of
#' the vertices of the graph considering only the links with a positive sign.
#'
#'@importFrom igraph layout.fruchterman.reingold subgraph.edges E
#'@export
layout_signed <- function(graph, seed=123){
graph.sub <- subgraph.edges(graph=graph,
eids=which(E(graph)$weight>0),
delete.vertices=FALSE)
set.seed(seed)
layout <- layout.fruchterman.reingold(graph.sub)
return(layout)
}
#------------------------------------------------------------------------------#
#' Get Communities of a Signed Weighted Graph
#'
#'@description Adaptation in R of the algorithm for detention of weighted communities with sign of a graph developed by Sergio Gomez.
#' https://deim.urv.cat/~sergio.gomez/radatools.php
#'
#'@importFrom igraph graph_from_adjacency_matrix write_graph make_clusters is.weighted is.directed as_adjacency_matrix make_clusters
#'@importFrom stringr str_split
#'
#'@export
community_detection <- function(graph){
#Check Graph
if(!is.weighted(graph)) stop("graph must be weighted graph")
if(is.directed(graph)) stop("graph must be undirected")
adj <- as_adjacency_matrix(graph, attr="weight", sparse=FALSE)
#End Checks
#get path of executable Communities_Detection.exe
path <- system.file("exec", package="MetaGenomicR", mustWork=TRUE)
path <- paste(path,"/",sep="")
#write graph in pajek format as request from executable
write_graph(graph = graph_from_adjacency_matrix(adjmatrix=adj, mode="undirected", weighted=TRUE),
file = paste(path,"graph.net",sep=""),
format ="pajek")
#path for graph and results files
path.graph <- paste(path,"graph.net",sep="")
path.result <- paste(path, "res.txt", sep="")
#save command for the execution
cmd <- paste(path,"Communities_Detection.exe none WS l 1",sep="")
cmd <- paste(cmd, path.graph, path.result)
#make communities detection
system(cmd, ignore.stdout=TRUE, ignore.stderr=TRUE)
#open file of results
res.file = file(path.result, "r")
#read and store results
file.lines <- readLines(res.file)
info <- file.lines[2]
modularity <- as.numeric(str_split(file.lines[3]," ")[[1]][3])
vertex.num <- as.numeric(str_split(file.lines[5]," ")[[1]][4])
comm.num <- as.numeric(str_split(file.lines[6]," ")[[1]][4])
comm.vert <- vector(mode="list",length=comm.num)
for(line in 8:(7+comm.num)){
i <- line-7
#comm[[i]]
tmp <- unlist(str_split(file.lines[line]," "))
tmp <- as.numeric(tmp[2:length(tmp)])
comm.vert[[i]] <- tmp
}
comm <- vector(mode="integer", length=vertex.num)
for(c in 1:comm.num){comm[comm.vert[[c]]] <- c}
#close connection
close(res.file)
#remove files
file.remove(path.graph)
file.remove(path.result)
#return the results as communities structure of igraph package
return(make_clusters(graph=graph,
membership=comm,
algorithm="signed weights (radatools)",
modularity=modularity))
}
# Get Weights
#------------------------------------------------------------------------------#
#------------------------------------------------------------------------------#
#' Spiec-Easi MB weights extraction
#'@description Extract the weight matrix of the Meinshausen and Bühlmann criterion
#'in the Spiec-Easi algorithm.
#'@importFrom SpiecEasi symBeta getOptBeta
#'
#'@param res.mb Spiec-Easi algorithm results with mb criteria
#'
#'@export
get_mb_weights <- function(res.mb){
#CHECK ARGUMENTS
if(class(res.mb)!="pulsar.refit") stop("class must be pulsar.refit")
if(res.mb$est$method!="mb") stop("must be mb method")
#END CHECK
# get weights matrix W
W <- as.matrix(symBeta(as.matrix(getOptBeta(res.mb)), mode='maxabs'))
colnames(res.mb$est$data) -> colnames(W) -> rownames(W)
W[W>0] <- 1
W[W<0] <- -1
return(W)
}
#' Spiec-Easi GLASSO weights extraction
#'@description Extract the weight matrix of the GLASSO criterion
#'in the Spiec-Easi algorithm.
#'@importFrom SpiecEasi getOptCov getRefit
#'@importFrom stats cov2cor
#'
#'@export
get_gl_weights <- function(res.gl, Refit=TRUE){
#CHECK ARGUMENTS
if(class(res.gl)!="pulsar.refit") stop("class must be pulsar.refit")
if(res.gl$est$method!="glasso") stop("must be gl method")
#END CHECK
# get weights matrix W
W <- cov2cor(as.matrix(getOptCov(res.gl)))
if(Refit){W <- W * as.matrix(getRefit(res.gl))}
colnames(res.gl$est$data) -> colnames(W) -> rownames(W)
return(W)
}
#'@title communities_colormap
#'@author Alessandro Fuschi
#'@importFrom qualpalr qualpal
#'@export
communities_colormap <- function(n=20){
assert((round(n)==n && n>0),"n must be an integer positive number")
set.seed(1)
colormap <- qualpal(n=n,colorspace=list(h=c(0,360),s=c(.25,1),l=c(.25,.75)))
colormap <- rownames(colormap$RGB)
colormap <- c(rgb(1,1,1,.8),colormap)
names(colormap) <- as.character(0:n)
return(colormap)
}
#'@title taxonomy_colormap
#'@author Alessandro Fuschi
#'@importFrom qualpalr qualpal
#'@export
taxonomy_colormap <- function(taxaID){
assert(is.character(taxaID),"taxaID must be character vector")
assert(identical(taxaID,unique(taxaID)),"find at least a duplicate in taxaID")
n <- length(taxaID)
set.seed(1)
colormap <- qualpal(n=n,colorspace=list(h=c(0,360),s=c(.25,1),l=c(.5,.9)))
colormap <- rownames(colormap$RGB)
names(colormap) <- taxaID
return(colormap)
}
#'@title remove.small.communities
#'@author Alessandro Fuschi
#'@importFrom igraph vcount induced.subgraph V
#'@export
setGeneric("remove.small.communities", function(obj,vertex.num.threshold,keep.discarded) standardGeneric("remove.small.communities"))
setMethod("remove.small.communities", c("MetaGenomicGraph","numeric","logical"), function(obj, vertex.num.threshold,keep.discarded){
#Checks Arguments
assert(vertex.num.threshold>0,"vertex.num.threshold must be integer greater than 0")
#End Checks
graph <- obj@netw
comm <- obj@comm
keep.comm.names <- as.numeric(names(sizes(comm)[sizes(comm)>=vertex.num.threshold]))
keep.comm.vids <- which(comm$membership %in% keep.comm.names)
if(keep.discarded){
graph.sub <- graph
comm.sub <- comm
comm.sub$membership[setdiff(1:comm$vcount, keep.comm.vids)] <- 0
comm.sub$modularity <- NA
} else {
graph.sub <- induced.subgraph(graph, V(graph)[keep.comm.vids])
comm.sub <- comm
comm.sub$membership <- comm$membership[keep.comm.vids]
comm.sub$vcount <- length(comm.sub$membership)
comm.sub$modularity <- NA
}
return(new("MetaGenomicGraph",
data=obj@data, meta=obj@meta, taxa=obj@taxa,
netw=graph.sub,comm=comm.sub))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.