R/functions.R
In abhik1368/netpredictor: Algorithms to search for missing link in a graph.
## Performing jaccard similarity between two entities
## Performing jaccard similarity between two entities
## Using parallel to make cluster
jaccard.sim <- function(df){
cl <- makeCluster(detectCores())
a <- suppressWarnings(snow::parMM(cl,df,t(df)))
b <- suppressWarnings(snow::parMM(cl,df,(1-t(df))))
c <- suppressWarnings(snow::parMM(cl,(1-df),t(df)))
sim <- a/(a + b + c)
# changing diagonals to 1 and removing NaN's and NA
diag(sim) <- 1
s = as.matrix(sim)
s[is.nan(s)] <- 0
return (s)
stopCluster(cl)
}
## Normalizing the matrix based on matrix columns
colNorm <- function(PTmatrix){
if (ncol(PTmatrix) > 1){
col_sum <- apply(PTmatrix, 2, sum)
col_sum_matrix <- matrix(rep(col_sum, nrow(PTmatrix)), ncol=ncol(PTmatrix), nrow=nrow(PTmatrix), byrow =T)
res <- PTmatrix/col_sum_matrix
res <- as.matrix(res)
res[is.na(res)] <- 0
return(res)
}else {
res <- PTmatrix[,1]/sum(PTmatrix[,1])
return (as.matrix(res))
}
}
## Get the transition matrix for a Bipartite Graph. Input a sequence similarity matrix (s1),
## chemical similarity matrix (s2) and drug target adjacency matrix (g1) where rows are
## protein targets and columns as Drugs.
## @export
tMat <- function(g1,s1,s2,normalise="chen"){
cl <- makeCluster(detectCores())
g1 <- t(g1)
seq<- as.matrix(s1) ## sequence similairty matrix normalised between 0 and 1
drugProt <- as.matrix(g1) ## drug target matrix
csim <- as.matrix(s2) ## drug similarity matrix normalised between 0 and 1
new.drug_drug <- suppressWarnings(snow::parMM(cl,drugProt,t(drugProt)))
new.prot_prot <- suppressWarnings(snow::parMM(cl,t(drugProt),drugProt))
#calculate drug-drug similarity based on shared proteins based on jaccard similarity
norm_drug <- jaccard.sim(drugProt)
# Jaccard similarity of two proteins based on shared compounds
norm_prot <- jaccard.sim(t(drugProt))
# Normalizing the matrices with equal weights
drug.similarity.final <- 0.5*(csim)+0.5*(norm_drug)
prot.similarity.final <- 0.5*(seq)+0.5*(norm_prot)
if(normalise == "laplace"){
D1 <- diag(x=(rowSums(drugProt))^(-0.5))
D2 <- diag(x=(colSums(drugProt))^(-0.5))
MTD1 <- suppressWarnings(snow::parMM(cl,D1,g1))
MTD <- suppressWarnings(snow::parMM(cl,MTD1,D2))
D3 <- diag(x=(rowSums(prot.similarity.final))^(-0.5))
MTT1 <- suppressWarnings(snow::parMM(cl,D3,seq))
MTT <- suppressWarnings(snow::parMM(cl,MTT1,D3))
D4 <- diag(x=(rowSums(drug.similarity.final))^(-0.5))
MDD1 <- suppressWarnings(snow::parMM(cl,D4,csim))
MDD <- suppressWarnings(snow::parMM(cl,MDD1,D4))
M1<-cbind(MTT,t(MTD))
M2<-cbind(MTD,MDD)
M <- rbind(M1,M2)
M <- as.matrix(M)
M[is.na(M)]<-0
n =c(colnames(g1),rownames(g1))
rownames(M) <- n
colnames(M) <- n
# Returning the final matrix
return(as.matrix(M))
}
if(normalise == "chen"){
ADT<-colSums(drugProt)
ATD<-rowSums(drugProt)
Sd<-rowSums(drug.similarity.final)
St<-rowSums(prot.similarity.final)
MDT<-mat.or.vec(nrow(drugProt),nrow(seq))
MTD<-mat.or.vec(nrow(seq),nrow(drugProt))
print (dim(MDT))
print (dim(MTD))
A<-t(drugProt)
D1 <- diag(x=(rowSums(drugProt))^(-0.5))
D2 <- diag(x=(colSums(drugProt))^(-0.5))
#MTD1 <- suppressWarnings(snow::parMM(cl,D1,g1))
##MTD <- suppressWarnings(snow::parMM(cl,MTD1,D2))
#MDT <- t(MTD)
D3 <- diag(x=(rowSums(prot.similarity.final))^(-0.5))
MTT1 <- suppressWarnings(snow::parMM(cl,D3,seq))
MTT <- suppressWarnings(snow::parMM(cl,MTT1,D3))
D4 <- diag(x=(rowSums(drug.similarity.final))^(-0.5))
MDD1 <- suppressWarnings(snow::parMM(cl,D4,csim))
MDD <- suppressWarnings(snow::parMM(cl,MDD1,D4))
for (i in 1:nrow(seq)){
for (j in 1:nrow(seq)){
if (ATD[i]==0){
MTT[i,j]<-prot.similarity.final[i,j]/St[i]
}
else{
MTT[i,j]<-(0.8*(prot.similarity.final[i,j])/St[i])
}
}
}
for (i in 1:ncol(drugProt)){
for (j in 1:ncol(drugProt)){
if (ADT[i]==0){
MDD[i,j]<-drug.similarity.final[i,j]/Sd[i]
}
else{
MDD[i,j]<-(0.8*(drug.similarity.final[i,j])/Sd[i])
}
}
}
for (i in 1:ncol(drugProt)){
for (j in 1:nrow(seq)){
if (ADT[i]!=0){
MDT[i,j]<- (0.2*A[i,j])/ADT[i]
}
else{
MDT[i,j]<-0
}
}
}
for (i in 1:nrow(seq)){
for (j in 1:ncol(drugProt)){
if (ATD[i]!=0){
MTD[i,j]<- (0.2*A[j,i])/ATD[i]
}
else{
MTD[i,j]<-0
}
}
}
M1<-cbind(MTT,MTD)
M2<-cbind(MDT,MDD)
M <- rbind(M1,M2)
M <- as.matrix(M)
M[is.na(M)]<-0
n =c(colnames(g1),rownames(g1))
rownames(M) <- n
colnames(M) <- n
# Returning the final matrix
return(as.matrix(M))
}
if(normalise == "none"){
MDD <- drug.similarity.final
MTT <- prot.similarity.final
MTD <- drugProt
M1<-cbind(MTT,t(MTD))
M2<-cbind(MTD,MDD)
M <- rbind(M1,M2)
M <- as.matrix(M)
M[is.na(M)]<-0
n =c(colnames(g1),rownames(g1))
rownames(M) <- n
colnames(M) <- n
M <- colNorm(M)
# Returning the final matrix
return(as.matrix(M))
}
on.exit(stopCluster(cl))
}
## performing random walk with restart both in parallel and non-parallel way.
## It takes the transition matrix W the initial matrix P0 matrix parameter for parallization
## and the number of cores to run parallelisation on. Also a restart parameter is available as r
## to get better results with different datasets one needs to tune restart parameter r.
rwr <- function(W,P0matrix,r=0.9){
stop_step <- 50
library(foreach)
library(doParallel)
registerDoParallel(detectCores())
message(paste(c("Executing RWR in parallel way .. \n")))
PTmatrix <- matrix(0, nrow=nrow(P0matrix), ncol=ncol(P0matrix))
pb <- txtProgressBar(min = 1, max = ncol(P0matrix), style = 3)
#flag = dCheckParallel(multicores = 4,verbose = T)
#if (flag){
j = 1
PTmatrix <- foreach::foreach(j=1:ncol(P0matrix), .inorder=T, .combine='cbind') %dopar% {
setTxtProgressBar(pb, j)
P0 <- P0matrix[,j]
## Initializing variables
step <- 0
PT <- P0
## Iterative update till convergence (delta<=1e-7)
while (step<=stop_step){
PX <- (1-r) * t(W) %*% PT + r * P0
delta <- sum(abs(PX-PT))
PT <- PX
step <- step+1
}
as.matrix(PT)
}
return(PTmatrix)
}
#' Get the communities in a graph
#' @title Get communities from a graph
#' @name get.Communities
#' @description Get the communities in a graph. It First decomposes the graph into components and then uses igraph's
#' various community detection algorithms to detect communities in a graph.
#' @param g igraph object
#' @param num.nodes Number of nodes to keep in the community.
#' @param calgo The community algorithm to use to find the communities.
#' @usage get.Communities(g,num.nodes = 2,calgo = walktrap.community)
#' @export
get.Communities <- function(g,num.nodes = 2,calgo = walktrap.community){
if (class(g)!= "igraph"){
stop("The function must apply to 'igraph' object.\n")
}
## decompose the graph
gps <- decompose.graph(g,min.vertices=num.nodes)
# Create a list to store list of graphs for community identification
result <- list()
nc<- sapply(gps, vcount)
# Get the indexes which has more than num.nodes
indx <- which(nc > num.nodes)
if (length(indx) < 1){
stop("No communities found less than num.nodes")
}
#g = Z$cgraph
community.significance.test <- function(graph, vs, ...) {
if (is.directed(graph)) stop("This method requires an undirected graph")
graph = g
vs <- membership(comm)
subgraph <- induced.subgraph(graph, vs)
in.degrees <- degree(subgraph)
out.degrees <- degree(graph, vs) - in.degrees
wilcox.test(in.degrees, out.degrees, ...)
}
j <- 1
for (i in indx){
comm <- cluster_walktrap(gps[[i]],weights=E(gps[[i]])$weight)
mem <- as.data.frame(as.matrix(membership(comm)))
colnames(mem)[1] <- "Members"
comm.sig <- community.significance.test(g,membership(comm))
result[[j]] = list(community = comm,
community.sig = comm.sig,
cgraph = gps[[i]],
members = mem)
j <- j+1
}
class(result) <- "community"
return (result)
}
#' plotting Communities
#' @title plot_Community
#' @description This uses an object of getCommuntiy class and extracts parameters to plot current communities.
#' @name plot_Community
#' @param gc receives the results from \code{get.Communities}
#' @param cols number of columns to show in the plot
#' @examples
#' \dontrun{
#' ## Run the Bipartite Random walk with Restart
#' Q = biNetwalk(g1,S,S1,normalise="laplace",parallel=FALSE,verbose=T)
#' ## Get the significant vertices which are in columns
#' Z = sig.net(data=A,g=g1,Amatrix=Q,num.permutation=100,adjp.cutoff=0.01,p.adjust.method="BH",parallel=FALSE)
#' ## Get the graph for plotting commnuties
#' g <- Z$cgraph
#' gp <- get.Communities(g)
#' ## Total number of communities
#' length(gp)
#' ## Plot the communities with 5 columns
#' plot_Community(gp,cols=5)
#' }
#' @export
plot_Community <- function(gc,cols=3) {
if (class(gc) != "community"){
stop("The function applies to community object.\n")
}
numPlots = length(gc)
r = ceiling(numPlots/cols)
if (numPlots==1) {
plot(gc[[1]]$community, gc[[1]]$cgraph)
} else {
par(mfrow = c(ceiling(numPlots/cols), cols))
par(cex = 0.6)
par(mar = c(4, 4, 2, 0.5), oma = c(2, 2, 2, 2))
for (i in 1:numPlots){
plot(gc[[i]]$community, gc[[i]]$cgraph,
main = paste(c("Community",i)),
layout=layout.fruchterman.reingold,
vertex.label.family="sans",
vertex.color= rainbow(10, .8, .8, alpha=.8))
}
}
par(mfrow = c(ceiling(numPlots/cols), cols))
mtext("Significant Communties", cex=1,font=2, col="Black", outer=TRUE)
}
#' Compute Degree Centrality of the Bipartite Graph
#' @title Degree centrality for Bipartite Graphs.
#' @name get.biDegreeCentrality
#' @description Measures graph degree centrality for bipartite graphs as well as Unipartite Graphs.The degree centrality of a vertex can be defined as fraction of incident edges over the number of all possible edges.
#' This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @param loops Either TRUE or FALSE keep loops or remove it.
#' @param SM Either TRUE or FALSE . If TRUE returns average centrality for rows and columns. If false returns degree centrality for all the nodes.
#' @return It returns a list of score of degree centrality of nodes at rows and nodes at columns.
#' @usage get.biDegreeCentrality(A,loops=FALSE,SM=TRUE)
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of the degree centrality the nodes in the rows are scored first followed by column vertices.
#' By deafult \code{SM=TRUE} which return the betweeness centrality for two different sets of nodes.
#' @export
get.biDegreeCentrality <-function(A,loops=FALSE,SM=TRUE){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (SM==TRUE){
if (!is.null(V(g)$type)){
V(g)$degree <- degree(g)
# determine maximum degrees for each vertex set
max_d_TRUE <- max(V(g)[type==TRUE]$degree)
max_d_FALSE <- max(V(g)[type==FALSE]$degree)
# determine centralization for TRUE vertex subset
g$sm.degree.centralization.columns <-sum(max_d_TRUE - V(g)[type==TRUE]$degree)/((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1))
# determine centralization for FALSE vertex subset
g$sm.degree.centralization.rows <-sum(max_d_FALSE - V(g)[type==FALSE]$degree)/((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1))
# return both values as list
return(list("SM.bi.degreeCentrality.columns"=g$sm.degree.centralization.columns,"SM.bi.degreeCentrality.rows"=g$sm.degree.centralization.rows))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
} else if (!is.null(V(g)$type)){
for (i in V(g)){
V(g)[i]$degreeCentrality <- degree(g,v=i)/length(V(g)[type==!V(g)[i]$type])
}
# return value vector as list item
return(list("bi.degreeCentrality"=V(g)$degreeCentrality))
}
else {
for (i in V(g)){
if (!loops){
V(g)[i]$degreeCentrality <- degree(g,v=i,loops=FALSE)/(length(V(g))-1)
}
else{
V(g)[i]$degreeCentrality <- degree(g,v=i,loops=TRUE)/(length(V(g)))
}
}
# return value vector as list item
return(list("bi.degreeCentrality"=V(g)$degreeCentrality))
}
}
#' Compute Graph density Bipartite Graph
#' @title Graph Density for Bipartite graphs.
#' @name get.biDensity
#' @description Measures graph density of bipartite graphs . The density captures the fraction of actual present edges over the number of all possible edges, given that multiple edges are not allowed.
#' This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @return It returns a list of score of degree centrality of nodes at rows and nodes at columns.
#' @usage get.biDensity(A)
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of the degree Density the nodes in the rows are scored first followed by column vertices.
#' @export
get.biDensity <-function(A){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (!is.null(V(g)$type)){
# return value as list item
if (!is.directed(g)){
return(list("bi.Density"=length(E(g))/(length(which(V(g)$type))*length(which(!V(g)$type)))))
}
else{
return(list("bi.Density"=length(E(g))/(2*length(which(V(g)$type))*length(which(!V(g)$type)))))
}
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
}
bipartite.closeness.centrality <- function(g){
if (class(g) != "igraph"){
stop("The function must apply to igraph object.\n")
}
if (!is.null(V(g)$type)){
# determine maximal raw scores for both vertex subsets
mrs_TRUE <- length(V(g)[type==FALSE]) + 2*length(V(g)[type==TRUE]) - 2
mrs_FALSE <- length(V(g)[type==TRUE]) + 2*length(V(g)[type==FALSE]) - 2
# get sum of all geodesic paths for each vertex
sp <- shortest.paths(g)
sp[sp==Inf] <- length(V(g))
rowsums_shortest_paths <- rowSums(sp)
# "bipartite" normalization of scores
for (i in V(g)){
if (V(g)[i]$type==TRUE){
V(g)[i]$closeness.centrality <- mrs_TRUE/rowsums_shortest_paths[i+1]
}
else{
V(g)[i]$closeness.centrality <- mrs_FALSE/rowsums_shortest_paths[i+1]
}
}
# return value as list
return(list("Bipartite.Closeness.Centrality"=V(g)$closeness.centrality))
}
else {
# boolean vertex attribute 'type' is required
cat("vertex attribute <type> is missing")
}
}
#' Closeness centrality for bipartite graphs
#' @title Closeness centrality for bipartite graphs.
#' @name get.biClosenessCentralit
#' @description Measures Closeness centrality of bipartite graphs . The closeness centrality of a vertex is inversely proportional to the total geodesic distance to all other vertices.
#' It also has way to calculate single mode betweeness centrality. Single mode centralization for bipartite graphs measures the extent to which vertices in one vertex subset are central relative only to other vertices in the same subset.
#' This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @param SM Either TRUE or FALSE . If TRUE returns average centrality for rows and columns. If false returns degree centrality for all the nodes.
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of closeness centrality of the nodes.By deafult \code{SM=TRUE} which return the closeness centrality for two different sets of nodes.
#' @export
get.biClosenessCentrality <-function(A,SM=TRUE){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (SM==TRUE){
if (!is.null(V(g)$type)){
# determine bipartite closeness centrality scores
d <- bipartite.closeness.centrality(g)[[1]]
d[is.na(d)] <- 0
V(g)$Bipartite.Closeness.Centrality <- d
# get maximum scores
max_c_TRUE <- max(V(g)[type==TRUE]$Bipartite.Closeness.Centrality)
max_c_FALSE <- max(V(g)[type==FALSE]$Bipartite.Closeness.Centrality)
# determine denominators for both vertex subsets
if (length(V(g)[type==FALSE])<length(V(g)[type==TRUE])){
denom_TRUE <- ((length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-2))/(2*length(V(g)[type==TRUE])-3) + ((length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-length(V(g)[type==FALSE])))/(length(V(g)[type==TRUE])+length(V(g)[type==FALSE])-2)
}
else{
denom_TRUE <- ((length(V(g)[type==TRUE])-2)*(length(V(g)[type==TRUE])-1))/(2*length(V(g)[type==TRUE])-3)
}
if (length(V(g)[type==TRUE])<length(V(g)[type==FALSE])){
denom_FALSE <- ((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-2))/(2*length(V(g)[type==FALSE])-3) + ((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-length(V(g)[type==TRUE])))/(length(V(g)[type==FALSE])+length(V(g)[type==TRUE])-2)
}
else{
denom_FALSE <- ((length(V(g)[type==FALSE])-2)*(length(V(g)[type==FALSE])-1))/(2*length(V(g)[type==FALSE])-3)
}
# determine centralization for TRUE vertex subset
g$sm.closeness.centralization.columns <-sum(max_c_TRUE - V(g)[type==TRUE]$Bipartite.Closeness.Centrality)/denom_TRUE
# determine centralization for FALSE vertex subset
g$sm.closeness.centralization.rows <-sum(max_c_FALSE - V(g)[type==FALSE]$Bipartite.Closeness.Centrality)/denom_FALSE
# return both values as list
return(list("SM.bi.ClosenessCentrality.columns"=g$sm.closeness.centralization.columns,"SM.bi.ClosenessCentrality.rows"=g$single.mode.closeness.centralization.rows))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
} else if (!is.null(V(g)$type)){
# determine maximal raw scores for both vertex subsets
mrs_TRUE <- length(V(g)[type==FALSE]) + 2*length(V(g)[type==TRUE]) - 2
mrs_FALSE <- length(V(g)[type==TRUE]) + 2*length(V(g)[type==FALSE]) - 2
# get sum of all geodesic paths for each vertex
sp <- shortest.paths(g)
sp[sp==Inf] <- length(V(g))
rowsums_shortest_paths <- rowSums(sp)
# "bipartite" normalization of scores
for (i in V(g)){
if (V(g)[i]$type==TRUE){
V(g)[i]$closeness.centrality.columns <- mrs_TRUE/rowsums_shortest_paths[i+1]
}
else{
V(g)[i]$closeness.centrality.rows <- mrs_FALSE/rowsums_shortest_paths[i+1]
}
}
# return values as list
return(list("bi.ClosenessCentrality.columns"=na.omit(V(g)$closeness.centrality.columns),"bi.ClosenessCentrality.rows"=na.omit(V(g)$closeness.centrality.rows)))
} else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
}
#' Betweeness centrality for bipartite graphs
#' @title Betweeness centrality for bipartite graphs.
#' @name get.biBetweenessCentrality
#' @description Measures Betweeness centrality of bipartite graphs. The betweenness centrality of a vertex may be roughly defined as the number of geodesic paths that pass through a given vertex,
#' weighted inversely by the total number of equivalent paths between the same two vertices, including those that do not pass through the given vertex.
#' It also has way to calculate single mode betweeness centrality. Single mode centralization for bipartite graphs measures the extent to which vertices in one vertex subset are central
#' relative only to other vertices in the same subset.This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @param SM Either TRUE or FALSE . If TRUE returns average centrality for rows and columns. If false returns betweeness centrality for all the nodes.
#' @usage get.biBetweenessCentrality(A,SM=TRUE)
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of closeness centrality of the nodes. By deafult \code{SM=TRUE} which return the betweeness centrality for two different sets of nodes.
#' @export
get.biBetweenessCentrality <- function(A,SM=TRUE){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (SM==TRUE){
if (!is.null(V(g)$type)){
# determine denominators for both vertex subsets
if (length(V(g)[type==FALSE])<length(V(g)[type==TRUE])){
denom_TRUE <- 2*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1)
}
else{
denom_TRUE <- (length(V(g)[type==TRUE])-1) * ( 0.5*length(V(g)[type==FALSE])*(length(V(g)[type==FALSE])-1) + 0.5*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==TRUE])-2) + (length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1) )
}
if (length(V(g)[type==TRUE])<length(V(g)[type==FALSE])){
denom_FALSE <- 2*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1)
}
else{
denom_FALSE <- (length(V(g)[type==FALSE])-1) * ( 0.5*length(V(g)[type==TRUE])*(length(V(g)[type==TRUE])-1) + 0.5*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==FALSE])-2) + (length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1) )
}
# determine raw betweenness centrality scores from igraph
V(g)$betweenness.raw <- betweenness(g,directed=FALSE)
# get maximum scores
max_c_TRUE <- max(V(g)[type==TRUE]$betweenness.raw)
max_c_FALSE <- max(V(g)[type==FALSE]$betweenness.raw)
# determine centralization for TRUE vertex subset
g$sm.betweenness.centralization.columns <-sum(max_c_TRUE - V(g)[type==TRUE]$betweenness.raw)/denom_TRUE
# determine centralization for FALSE vertex subset
g$sm.betweenness.centralization.rows <-sum(max_c_FALSE - V(g)[type==FALSE]$betweenness.raw)/denom_FALSE
# return both values as list
return(list("SM.bi.BetweenessCentrality.columns"=g$sm.betweenness.centralization.columns,"SM.bi.BetweenessCentrality.rows"=g$sm.betweenness.centralization.rows))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
} else if (!is.null(V(g)$type)){
# determine maximal raw scores for both vertex subsets
if (length(V(g)[type==FALSE])<length(V(g)[type==TRUE])){
mrs_TRUE <- 2*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1)
}
else{
mrs_TRUE <- 0.5*length(V(g)[type==FALSE])*(length(V(g)[type==FALSE])-1)+0.5*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==TRUE])-2)+(length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1)
}
if (length(V(g)[type==TRUE])<length(V(g)[type==FALSE])){
mrs_FALSE <- 0.5*length(V(g)[type==TRUE])*(length(V(g)[type==TRUE])-1)+0.5*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==FALSE])-2)+(length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1)
}
else{
mrs_FALSE <- 2*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1)
}
# get raw betweenness centrality scores from igraph
betweenness_rs <- betweenness(g,directed=FALSE)
# "bipartite" normalization of scores
for (i in V(g)){
if (V(g)[i]$type==TRUE){
V(g)[i]$betweenness.centrality.columns <- betweenness_rs[i+1]/mrs_TRUE
}
else{
V(g)[i]$betweenness.centrality.rows <- betweenness_rs[i+1]/mrs_FALSE
}
}
# return value as list
return(list("bi.BetweenessCentrality.columns"=na.omit(V(g)$betweenness.centrality.columns),"bi.BetweenessCentrality.rows"=na.omit(V(g)$betweenness.centrality.rows)))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
}
#' Projecting Bipartite Networks
#' @title Projecting Bipartite Networks.
#' @name get.biWeightedProjection
#' @description Projecting the Bipartite Network based on the vertex type (\code{TRUE} or \code{FALSE}).The mode \code{"shared-neighbours"} adds the number of shared neighbours in the bipartite graph to the edge linking two neighbours in the monopartite projection.
#' \code{"newman"} adopts the weighting scheme presented by Newman (2001), that weighs the contributions of shared neighbours by the size of their linkage profiles minus one.
#' @param A Biparite matrix object
#' @param vertex Either TRUE of FALSE where if FALSE then off-diagonal elements represent edge weights, diagonal elements vertex weights
#' @param mode \code{'shared-neighbours'} and \code{'newman'}
#' @usage get.biWeightedProjection(A,vertex=FALSE,mode='shared-neighbours')
#' @references
#' \itemize{
#' \item Newman M (2001) Scientific collaboration networks. II. Shortest paths, weighted networks, and centrality. Physical Review E 64.
#' }
#' @return Returns a igraph object of the projected matrix with edge weights.
#' @export
get.biWeightedProjection <- function(A,vertex=FALSE,mode='shared-neighbours'){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (!bipartite.mapping(g)$res){
stop("The function applies to bipartite graphs.\n")
}
if (is.bipartite(g)){
V(g)$type <- bipartite.mapping(g)$type
if (vertex) proj <- bipartite.projection(g)[[2]]
else proj <- bipartite.projection(g)[[1]]
V(proj)$name <- V(g)[type==vertex]$name
## number of shared vertices as edge weights
if (mode=='shared-neighbours'){
incidenceMatrix <- get.incidence(g) ## TODO: might need some more testing if type and vType match
if (!vertex){
## off-diagonal elements represent edge weights, diagonal elements vertex weights (currently unused)
m <- incidenceMatrix %*% t(incidenceMatrix)
}
else {
m <- t(incidenceMatrix) %*% incidenceMatrix
}
## assign edge weights from matrix m
if (length(E(proj))>0){
for (i in 1:(length(E(proj)))){
E(proj)[i]$weight <- m[V(proj)[get.edge(proj,i)[[1]]]$name, V(proj)[get.edge(proj,i)[[2]]]$name]
}
}
return(proj)
}
## Newman2001 Formula
if (mode=='newman'){
incidenceMatrix <- get.incidence(g)
if (!vertex){
incidenceMatrix <- t(incidenceMatrix)
if (length(c(which(rowSums(incidenceMatrix)==1)))!=0){
m2 <- incidenceMatrix[c(which(rowSums(incidenceMatrix)==1))*(-1),]
}
else{
m2 <- incidenceMatrix
}
m <- t(m2) %*% (m2*(1/(rowSums(m2)-1)))
}
else {
if (length(c(which(rowSums(incidenceMatrix)==1)))!=0){
m2 <- incidenceMatrix[c(which(rowSums(incidenceMatrix)==1))*(-1),]
}
else{
m2 <- incidenceMatrix
}
m <- t(m2) %*% (m2*(1/(rowSums(m2)-1)))
}
if (length(E(proj))>0){
for (i in 1:(length(E(proj)))){
E(proj)[i]$weight <- m[V(proj)[get.edge(proj,i)[[1]]]$name, V(proj)[get.edge(proj,i)[[2]]]$name]
}
}
return(proj)
}
}
}
abhik1368/netpredictor documentation built on May 10, 2019, 4:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## Performing jaccard similarity between two entities
## Performing jaccard similarity between two entities
## Using parallel to make cluster
jaccard.sim <- function(df){
cl <- makeCluster(detectCores())
a <- suppressWarnings(snow::parMM(cl,df,t(df)))
b <- suppressWarnings(snow::parMM(cl,df,(1-t(df))))
c <- suppressWarnings(snow::parMM(cl,(1-df),t(df)))
sim <- a/(a + b + c)
# changing diagonals to 1 and removing NaN's and NA
diag(sim) <- 1
s = as.matrix(sim)
s[is.nan(s)] <- 0
return (s)
stopCluster(cl)
}
## Normalizing the matrix based on matrix columns
colNorm <- function(PTmatrix){
if (ncol(PTmatrix) > 1){
col_sum <- apply(PTmatrix, 2, sum)
col_sum_matrix <- matrix(rep(col_sum, nrow(PTmatrix)), ncol=ncol(PTmatrix), nrow=nrow(PTmatrix), byrow =T)
res <- PTmatrix/col_sum_matrix
res <- as.matrix(res)
res[is.na(res)] <- 0
return(res)
}else {
res <- PTmatrix[,1]/sum(PTmatrix[,1])
return (as.matrix(res))
}
}
## Get the transition matrix for a Bipartite Graph. Input a sequence similarity matrix (s1),
## chemical similarity matrix (s2) and drug target adjacency matrix (g1) where rows are
## protein targets and columns as Drugs.
## @export
tMat <- function(g1,s1,s2,normalise="chen"){
cl <- makeCluster(detectCores())
g1 <- t(g1)
seq<- as.matrix(s1) ## sequence similairty matrix normalised between 0 and 1
drugProt <- as.matrix(g1) ## drug target matrix
csim <- as.matrix(s2) ## drug similarity matrix normalised between 0 and 1
new.drug_drug <- suppressWarnings(snow::parMM(cl,drugProt,t(drugProt)))
new.prot_prot <- suppressWarnings(snow::parMM(cl,t(drugProt),drugProt))
#calculate drug-drug similarity based on shared proteins based on jaccard similarity
norm_drug <- jaccard.sim(drugProt)
# Jaccard similarity of two proteins based on shared compounds
norm_prot <- jaccard.sim(t(drugProt))
# Normalizing the matrices with equal weights
drug.similarity.final <- 0.5*(csim)+0.5*(norm_drug)
prot.similarity.final <- 0.5*(seq)+0.5*(norm_prot)
if(normalise == "laplace"){
D1 <- diag(x=(rowSums(drugProt))^(-0.5))
D2 <- diag(x=(colSums(drugProt))^(-0.5))
MTD1 <- suppressWarnings(snow::parMM(cl,D1,g1))
MTD <- suppressWarnings(snow::parMM(cl,MTD1,D2))
D3 <- diag(x=(rowSums(prot.similarity.final))^(-0.5))
MTT1 <- suppressWarnings(snow::parMM(cl,D3,seq))
MTT <- suppressWarnings(snow::parMM(cl,MTT1,D3))
D4 <- diag(x=(rowSums(drug.similarity.final))^(-0.5))
MDD1 <- suppressWarnings(snow::parMM(cl,D4,csim))
MDD <- suppressWarnings(snow::parMM(cl,MDD1,D4))
M1<-cbind(MTT,t(MTD))
M2<-cbind(MTD,MDD)
M <- rbind(M1,M2)
M <- as.matrix(M)
M[is.na(M)]<-0
n =c(colnames(g1),rownames(g1))
rownames(M) <- n
colnames(M) <- n
# Returning the final matrix
return(as.matrix(M))
}
if(normalise == "chen"){
ADT<-colSums(drugProt)
ATD<-rowSums(drugProt)
Sd<-rowSums(drug.similarity.final)
St<-rowSums(prot.similarity.final)
MDT<-mat.or.vec(nrow(drugProt),nrow(seq))
MTD<-mat.or.vec(nrow(seq),nrow(drugProt))
print (dim(MDT))
print (dim(MTD))
A<-t(drugProt)
D1 <- diag(x=(rowSums(drugProt))^(-0.5))
D2 <- diag(x=(colSums(drugProt))^(-0.5))
#MTD1 <- suppressWarnings(snow::parMM(cl,D1,g1))
##MTD <- suppressWarnings(snow::parMM(cl,MTD1,D2))
#MDT <- t(MTD)
D3 <- diag(x=(rowSums(prot.similarity.final))^(-0.5))
MTT1 <- suppressWarnings(snow::parMM(cl,D3,seq))
MTT <- suppressWarnings(snow::parMM(cl,MTT1,D3))
D4 <- diag(x=(rowSums(drug.similarity.final))^(-0.5))
MDD1 <- suppressWarnings(snow::parMM(cl,D4,csim))
MDD <- suppressWarnings(snow::parMM(cl,MDD1,D4))
for (i in 1:nrow(seq)){
for (j in 1:nrow(seq)){
if (ATD[i]==0){
MTT[i,j]<-prot.similarity.final[i,j]/St[i]
}
else{
MTT[i,j]<-(0.8*(prot.similarity.final[i,j])/St[i])
}
}
}
for (i in 1:ncol(drugProt)){
for (j in 1:ncol(drugProt)){
if (ADT[i]==0){
MDD[i,j]<-drug.similarity.final[i,j]/Sd[i]
}
else{
MDD[i,j]<-(0.8*(drug.similarity.final[i,j])/Sd[i])
}
}
}
for (i in 1:ncol(drugProt)){
for (j in 1:nrow(seq)){
if (ADT[i]!=0){
MDT[i,j]<- (0.2*A[i,j])/ADT[i]
}
else{
MDT[i,j]<-0
}
}
}
for (i in 1:nrow(seq)){
for (j in 1:ncol(drugProt)){
if (ATD[i]!=0){
MTD[i,j]<- (0.2*A[j,i])/ATD[i]
}
else{
MTD[i,j]<-0
}
}
}
M1<-cbind(MTT,MTD)
M2<-cbind(MDT,MDD)
M <- rbind(M1,M2)
M <- as.matrix(M)
M[is.na(M)]<-0
n =c(colnames(g1),rownames(g1))
rownames(M) <- n
colnames(M) <- n
# Returning the final matrix
return(as.matrix(M))
}
if(normalise == "none"){
MDD <- drug.similarity.final
MTT <- prot.similarity.final
MTD <- drugProt
M1<-cbind(MTT,t(MTD))
M2<-cbind(MTD,MDD)
M <- rbind(M1,M2)
M <- as.matrix(M)
M[is.na(M)]<-0
n =c(colnames(g1),rownames(g1))
rownames(M) <- n
colnames(M) <- n
M <- colNorm(M)
# Returning the final matrix
return(as.matrix(M))
}
on.exit(stopCluster(cl))
}
## performing random walk with restart both in parallel and non-parallel way.
## It takes the transition matrix W the initial matrix P0 matrix parameter for parallization
## and the number of cores to run parallelisation on. Also a restart parameter is available as r
## to get better results with different datasets one needs to tune restart parameter r.
rwr <- function(W,P0matrix,r=0.9){
stop_step <- 50
library(foreach)
library(doParallel)
registerDoParallel(detectCores())
message(paste(c("Executing RWR in parallel way .. \n")))
PTmatrix <- matrix(0, nrow=nrow(P0matrix), ncol=ncol(P0matrix))
pb <- txtProgressBar(min = 1, max = ncol(P0matrix), style = 3)
#flag = dCheckParallel(multicores = 4,verbose = T)
#if (flag){
j = 1
PTmatrix <- foreach::foreach(j=1:ncol(P0matrix), .inorder=T, .combine='cbind') %dopar% {
setTxtProgressBar(pb, j)
P0 <- P0matrix[,j]
## Initializing variables
step <- 0
PT <- P0
## Iterative update till convergence (delta<=1e-7)
while (step<=stop_step){
PX <- (1-r) * t(W) %*% PT + r * P0
delta <- sum(abs(PX-PT))
PT <- PX
step <- step+1
}
as.matrix(PT)
}
return(PTmatrix)
}
#' Get the communities in a graph
#' @title Get communities from a graph
#' @name get.Communities
#' @description Get the communities in a graph. It First decomposes the graph into components and then uses igraph's
#' various community detection algorithms to detect communities in a graph.
#' @param g igraph object
#' @param num.nodes Number of nodes to keep in the community.
#' @param calgo The community algorithm to use to find the communities.
#' @usage get.Communities(g,num.nodes = 2,calgo = walktrap.community)
#' @export
get.Communities <- function(g,num.nodes = 2,calgo = walktrap.community){
if (class(g)!= "igraph"){
stop("The function must apply to 'igraph' object.\n")
}
## decompose the graph
gps <- decompose.graph(g,min.vertices=num.nodes)
# Create a list to store list of graphs for community identification
result <- list()
nc<- sapply(gps, vcount)
# Get the indexes which has more than num.nodes
indx <- which(nc > num.nodes)
if (length(indx) < 1){
stop("No communities found less than num.nodes")
}
#g = Z$cgraph
community.significance.test <- function(graph, vs, ...) {
if (is.directed(graph)) stop("This method requires an undirected graph")
graph = g
vs <- membership(comm)
subgraph <- induced.subgraph(graph, vs)
in.degrees <- degree(subgraph)
out.degrees <- degree(graph, vs) - in.degrees
wilcox.test(in.degrees, out.degrees, ...)
}
j <- 1
for (i in indx){
comm <- cluster_walktrap(gps[[i]],weights=E(gps[[i]])$weight)
mem <- as.data.frame(as.matrix(membership(comm)))
colnames(mem)[1] <- "Members"
comm.sig <- community.significance.test(g,membership(comm))
result[[j]] = list(community = comm,
community.sig = comm.sig,
cgraph = gps[[i]],
members = mem)
j <- j+1
}
class(result) <- "community"
return (result)
}
#' plotting Communities
#' @title plot_Community
#' @description This uses an object of getCommuntiy class and extracts parameters to plot current communities.
#' @name plot_Community
#' @param gc receives the results from \code{get.Communities}
#' @param cols number of columns to show in the plot
#' @examples
#' \dontrun{
#' ## Run the Bipartite Random walk with Restart
#' Q = biNetwalk(g1,S,S1,normalise="laplace",parallel=FALSE,verbose=T)
#' ## Get the significant vertices which are in columns
#' Z = sig.net(data=A,g=g1,Amatrix=Q,num.permutation=100,adjp.cutoff=0.01,p.adjust.method="BH",parallel=FALSE)
#' ## Get the graph for plotting commnuties
#' g <- Z$cgraph
#' gp <- get.Communities(g)
#' ## Total number of communities
#' length(gp)
#' ## Plot the communities with 5 columns
#' plot_Community(gp,cols=5)
#' }
#' @export
plot_Community <- function(gc,cols=3) {
if (class(gc) != "community"){
stop("The function applies to community object.\n")
}
numPlots = length(gc)
r = ceiling(numPlots/cols)
if (numPlots==1) {
plot(gc[[1]]$community, gc[[1]]$cgraph)
} else {
par(mfrow = c(ceiling(numPlots/cols), cols))
par(cex = 0.6)
par(mar = c(4, 4, 2, 0.5), oma = c(2, 2, 2, 2))
for (i in 1:numPlots){
plot(gc[[i]]$community, gc[[i]]$cgraph,
main = paste(c("Community",i)),
layout=layout.fruchterman.reingold,
vertex.label.family="sans",
vertex.color= rainbow(10, .8, .8, alpha=.8))
}
}
par(mfrow = c(ceiling(numPlots/cols), cols))
mtext("Significant Communties", cex=1,font=2, col="Black", outer=TRUE)
}
#' Compute Degree Centrality of the Bipartite Graph
#' @title Degree centrality for Bipartite Graphs.
#' @name get.biDegreeCentrality
#' @description Measures graph degree centrality for bipartite graphs as well as Unipartite Graphs.The degree centrality of a vertex can be defined as fraction of incident edges over the number of all possible edges.
#' This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @param loops Either TRUE or FALSE keep loops or remove it.
#' @param SM Either TRUE or FALSE . If TRUE returns average centrality for rows and columns. If false returns degree centrality for all the nodes.
#' @return It returns a list of score of degree centrality of nodes at rows and nodes at columns.
#' @usage get.biDegreeCentrality(A,loops=FALSE,SM=TRUE)
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of the degree centrality the nodes in the rows are scored first followed by column vertices.
#' By deafult \code{SM=TRUE} which return the betweeness centrality for two different sets of nodes.
#' @export
get.biDegreeCentrality <-function(A,loops=FALSE,SM=TRUE){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (SM==TRUE){
if (!is.null(V(g)$type)){
V(g)$degree <- degree(g)
# determine maximum degrees for each vertex set
max_d_TRUE <- max(V(g)[type==TRUE]$degree)
max_d_FALSE <- max(V(g)[type==FALSE]$degree)
# determine centralization for TRUE vertex subset
g$sm.degree.centralization.columns <-sum(max_d_TRUE - V(g)[type==TRUE]$degree)/((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1))
# determine centralization for FALSE vertex subset
g$sm.degree.centralization.rows <-sum(max_d_FALSE - V(g)[type==FALSE]$degree)/((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1))
# return both values as list
return(list("SM.bi.degreeCentrality.columns"=g$sm.degree.centralization.columns,"SM.bi.degreeCentrality.rows"=g$sm.degree.centralization.rows))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
} else if (!is.null(V(g)$type)){
for (i in V(g)){
V(g)[i]$degreeCentrality <- degree(g,v=i)/length(V(g)[type==!V(g)[i]$type])
}
# return value vector as list item
return(list("bi.degreeCentrality"=V(g)$degreeCentrality))
}
else {
for (i in V(g)){
if (!loops){
V(g)[i]$degreeCentrality <- degree(g,v=i,loops=FALSE)/(length(V(g))-1)
}
else{
V(g)[i]$degreeCentrality <- degree(g,v=i,loops=TRUE)/(length(V(g)))
}
}
# return value vector as list item
return(list("bi.degreeCentrality"=V(g)$degreeCentrality))
}
}
#' Compute Graph density Bipartite Graph
#' @title Graph Density for Bipartite graphs.
#' @name get.biDensity
#' @description Measures graph density of bipartite graphs . The density captures the fraction of actual present edges over the number of all possible edges, given that multiple edges are not allowed.
#' This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @return It returns a list of score of degree centrality of nodes at rows and nodes at columns.
#' @usage get.biDensity(A)
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of the degree Density the nodes in the rows are scored first followed by column vertices.
#' @export
get.biDensity <-function(A){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (!is.null(V(g)$type)){
# return value as list item
if (!is.directed(g)){
return(list("bi.Density"=length(E(g))/(length(which(V(g)$type))*length(which(!V(g)$type)))))
}
else{
return(list("bi.Density"=length(E(g))/(2*length(which(V(g)$type))*length(which(!V(g)$type)))))
}
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
}
bipartite.closeness.centrality <- function(g){
if (class(g) != "igraph"){
stop("The function must apply to igraph object.\n")
}
if (!is.null(V(g)$type)){
# determine maximal raw scores for both vertex subsets
mrs_TRUE <- length(V(g)[type==FALSE]) + 2*length(V(g)[type==TRUE]) - 2
mrs_FALSE <- length(V(g)[type==TRUE]) + 2*length(V(g)[type==FALSE]) - 2
# get sum of all geodesic paths for each vertex
sp <- shortest.paths(g)
sp[sp==Inf] <- length(V(g))
rowsums_shortest_paths <- rowSums(sp)
# "bipartite" normalization of scores
for (i in V(g)){
if (V(g)[i]$type==TRUE){
V(g)[i]$closeness.centrality <- mrs_TRUE/rowsums_shortest_paths[i+1]
}
else{
V(g)[i]$closeness.centrality <- mrs_FALSE/rowsums_shortest_paths[i+1]
}
}
# return value as list
return(list("Bipartite.Closeness.Centrality"=V(g)$closeness.centrality))
}
else {
# boolean vertex attribute 'type' is required
cat("vertex attribute <type> is missing")
}
}
#' Closeness centrality for bipartite graphs
#' @title Closeness centrality for bipartite graphs.
#' @name get.biClosenessCentralit
#' @description Measures Closeness centrality of bipartite graphs . The closeness centrality of a vertex is inversely proportional to the total geodesic distance to all other vertices.
#' It also has way to calculate single mode betweeness centrality. Single mode centralization for bipartite graphs measures the extent to which vertices in one vertex subset are central relative only to other vertices in the same subset.
#' This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @param SM Either TRUE or FALSE . If TRUE returns average centrality for rows and columns. If false returns degree centrality for all the nodes.
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of closeness centrality of the nodes.By deafult \code{SM=TRUE} which return the closeness centrality for two different sets of nodes.
#' @export
get.biClosenessCentrality <-function(A,SM=TRUE){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (SM==TRUE){
if (!is.null(V(g)$type)){
# determine bipartite closeness centrality scores
d <- bipartite.closeness.centrality(g)[[1]]
d[is.na(d)] <- 0
V(g)$Bipartite.Closeness.Centrality <- d
# get maximum scores
max_c_TRUE <- max(V(g)[type==TRUE]$Bipartite.Closeness.Centrality)
max_c_FALSE <- max(V(g)[type==FALSE]$Bipartite.Closeness.Centrality)
# determine denominators for both vertex subsets
if (length(V(g)[type==FALSE])<length(V(g)[type==TRUE])){
denom_TRUE <- ((length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-2))/(2*length(V(g)[type==TRUE])-3) + ((length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-length(V(g)[type==FALSE])))/(length(V(g)[type==TRUE])+length(V(g)[type==FALSE])-2)
}
else{
denom_TRUE <- ((length(V(g)[type==TRUE])-2)*(length(V(g)[type==TRUE])-1))/(2*length(V(g)[type==TRUE])-3)
}
if (length(V(g)[type==TRUE])<length(V(g)[type==FALSE])){
denom_FALSE <- ((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-2))/(2*length(V(g)[type==FALSE])-3) + ((length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-length(V(g)[type==TRUE])))/(length(V(g)[type==FALSE])+length(V(g)[type==TRUE])-2)
}
else{
denom_FALSE <- ((length(V(g)[type==FALSE])-2)*(length(V(g)[type==FALSE])-1))/(2*length(V(g)[type==FALSE])-3)
}
# determine centralization for TRUE vertex subset
g$sm.closeness.centralization.columns <-sum(max_c_TRUE - V(g)[type==TRUE]$Bipartite.Closeness.Centrality)/denom_TRUE
# determine centralization for FALSE vertex subset
g$sm.closeness.centralization.rows <-sum(max_c_FALSE - V(g)[type==FALSE]$Bipartite.Closeness.Centrality)/denom_FALSE
# return both values as list
return(list("SM.bi.ClosenessCentrality.columns"=g$sm.closeness.centralization.columns,"SM.bi.ClosenessCentrality.rows"=g$single.mode.closeness.centralization.rows))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
} else if (!is.null(V(g)$type)){
# determine maximal raw scores for both vertex subsets
mrs_TRUE <- length(V(g)[type==FALSE]) + 2*length(V(g)[type==TRUE]) - 2
mrs_FALSE <- length(V(g)[type==TRUE]) + 2*length(V(g)[type==FALSE]) - 2
# get sum of all geodesic paths for each vertex
sp <- shortest.paths(g)
sp[sp==Inf] <- length(V(g))
rowsums_shortest_paths <- rowSums(sp)
# "bipartite" normalization of scores
for (i in V(g)){
if (V(g)[i]$type==TRUE){
V(g)[i]$closeness.centrality.columns <- mrs_TRUE/rowsums_shortest_paths[i+1]
}
else{
V(g)[i]$closeness.centrality.rows <- mrs_FALSE/rowsums_shortest_paths[i+1]
}
}
# return values as list
return(list("bi.ClosenessCentrality.columns"=na.omit(V(g)$closeness.centrality.columns),"bi.ClosenessCentrality.rows"=na.omit(V(g)$closeness.centrality.rows)))
} else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
}
#' Betweeness centrality for bipartite graphs
#' @title Betweeness centrality for bipartite graphs.
#' @name get.biBetweenessCentrality
#' @description Measures Betweeness centrality of bipartite graphs. The betweenness centrality of a vertex may be roughly defined as the number of geodesic paths that pass through a given vertex,
#' weighted inversely by the total number of equivalent paths between the same two vertices, including those that do not pass through the given vertex.
#' It also has way to calculate single mode betweeness centrality. Single mode centralization for bipartite graphs measures the extent to which vertices in one vertex subset are central
#' relative only to other vertices in the same subset.This function is the adaption to bipartite graphs as presented in Borgatti and Everett (1997).
#' @param A Biparite matrix object
#' @param SM Either TRUE or FALSE . If TRUE returns average centrality for rows and columns. If false returns betweeness centrality for all the nodes.
#' @usage get.biBetweenessCentrality(A,SM=TRUE)
#' @references
#' \itemize{
#' \item Borgatti, S. P. and Everett, M. G. (1997) Network analysis of 2--mode data. Social Networks \bold{19}, 243--269.
#' }
#' @return Returns a list scores of closeness centrality of the nodes. By deafult \code{SM=TRUE} which return the betweeness centrality for two different sets of nodes.
#' @export
get.biBetweenessCentrality <- function(A,SM=TRUE){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (SM==TRUE){
if (!is.null(V(g)$type)){
# determine denominators for both vertex subsets
if (length(V(g)[type==FALSE])<length(V(g)[type==TRUE])){
denom_TRUE <- 2*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1)
}
else{
denom_TRUE <- (length(V(g)[type==TRUE])-1) * ( 0.5*length(V(g)[type==FALSE])*(length(V(g)[type==FALSE])-1) + 0.5*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==TRUE])-2) + (length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1) )
}
if (length(V(g)[type==TRUE])<length(V(g)[type==FALSE])){
denom_FALSE <- 2*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1)
}
else{
denom_FALSE <- (length(V(g)[type==FALSE])-1) * ( 0.5*length(V(g)[type==TRUE])*(length(V(g)[type==TRUE])-1) + 0.5*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==FALSE])-2) + (length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1) )
}
# determine raw betweenness centrality scores from igraph
V(g)$betweenness.raw <- betweenness(g,directed=FALSE)
# get maximum scores
max_c_TRUE <- max(V(g)[type==TRUE]$betweenness.raw)
max_c_FALSE <- max(V(g)[type==FALSE]$betweenness.raw)
# determine centralization for TRUE vertex subset
g$sm.betweenness.centralization.columns <-sum(max_c_TRUE - V(g)[type==TRUE]$betweenness.raw)/denom_TRUE
# determine centralization for FALSE vertex subset
g$sm.betweenness.centralization.rows <-sum(max_c_FALSE - V(g)[type==FALSE]$betweenness.raw)/denom_FALSE
# return both values as list
return(list("SM.bi.BetweenessCentrality.columns"=g$sm.betweenness.centralization.columns,"SM.bi.BetweenessCentrality.rows"=g$sm.betweenness.centralization.rows))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
} else if (!is.null(V(g)$type)){
# determine maximal raw scores for both vertex subsets
if (length(V(g)[type==FALSE])<length(V(g)[type==TRUE])){
mrs_TRUE <- 2*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1)
}
else{
mrs_TRUE <- 0.5*length(V(g)[type==FALSE])*(length(V(g)[type==FALSE])-1)+0.5*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==TRUE])-2)+(length(V(g)[type==FALSE])-1)*(length(V(g)[type==TRUE])-1)
}
if (length(V(g)[type==TRUE])<length(V(g)[type==FALSE])){
mrs_FALSE <- 0.5*length(V(g)[type==TRUE])*(length(V(g)[type==TRUE])-1)+0.5*(length(V(g)[type==FALSE])-1)*(length(V(g)[type==FALSE])-2)+(length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1)
}
else{
mrs_FALSE <- 2*(length(V(g)[type==TRUE])-1)*(length(V(g)[type==FALSE])-1)
}
# get raw betweenness centrality scores from igraph
betweenness_rs <- betweenness(g,directed=FALSE)
# "bipartite" normalization of scores
for (i in V(g)){
if (V(g)[i]$type==TRUE){
V(g)[i]$betweenness.centrality.columns <- betweenness_rs[i+1]/mrs_TRUE
}
else{
V(g)[i]$betweenness.centrality.rows <- betweenness_rs[i+1]/mrs_FALSE
}
}
# return value as list
return(list("bi.BetweenessCentrality.columns"=na.omit(V(g)$betweenness.centrality.columns),"bi.BetweenessCentrality.rows"=na.omit(V(g)$betweenness.centrality.rows)))
}
else {
# boolean vertex attribute 'type' is required
warning("vertex attribute <type> is missing")
}
}
#' Projecting Bipartite Networks
#' @title Projecting Bipartite Networks.
#' @name get.biWeightedProjection
#' @description Projecting the Bipartite Network based on the vertex type (\code{TRUE} or \code{FALSE}).The mode \code{"shared-neighbours"} adds the number of shared neighbours in the bipartite graph to the edge linking two neighbours in the monopartite projection.
#' \code{"newman"} adopts the weighting scheme presented by Newman (2001), that weighs the contributions of shared neighbours by the size of their linkage profiles minus one.
#' @param A Biparite matrix object
#' @param vertex Either TRUE of FALSE where if FALSE then off-diagonal elements represent edge weights, diagonal elements vertex weights
#' @param mode \code{'shared-neighbours'} and \code{'newman'}
#' @usage get.biWeightedProjection(A,vertex=FALSE,mode='shared-neighbours')
#' @references
#' \itemize{
#' \item Newman M (2001) Scientific collaboration networks. II. Shortest paths, weighted networks, and centrality. Physical Review E 64.
#' }
#' @return Returns a igraph object of the projected matrix with edge weights.
#' @export
get.biWeightedProjection <- function(A,vertex=FALSE,mode='shared-neighbours'){
if (class(A) != "matrix"){
stop("The function must apply to matrix object.\n")
}
g <- graph.incidence(A,mode = 'all')
if (!bipartite.mapping(g)$res){
stop("The function applies to bipartite graphs.\n")
}
if (is.bipartite(g)){
V(g)$type <- bipartite.mapping(g)$type
if (vertex) proj <- bipartite.projection(g)[[2]]
else proj <- bipartite.projection(g)[[1]]
V(proj)$name <- V(g)[type==vertex]$name
## number of shared vertices as edge weights
if (mode=='shared-neighbours'){
incidenceMatrix <- get.incidence(g) ## TODO: might need some more testing if type and vType match
if (!vertex){
## off-diagonal elements represent edge weights, diagonal elements vertex weights (currently unused)
m <- incidenceMatrix %*% t(incidenceMatrix)
}
else {
m <- t(incidenceMatrix) %*% incidenceMatrix
}
## assign edge weights from matrix m
if (length(E(proj))>0){
for (i in 1:(length(E(proj)))){
E(proj)[i]$weight <- m[V(proj)[get.edge(proj,i)[[1]]]$name, V(proj)[get.edge(proj,i)[[2]]]$name]
}
}
return(proj)
}
## Newman2001 Formula
if (mode=='newman'){
incidenceMatrix <- get.incidence(g)
if (!vertex){
incidenceMatrix <- t(incidenceMatrix)
if (length(c(which(rowSums(incidenceMatrix)==1)))!=0){
m2 <- incidenceMatrix[c(which(rowSums(incidenceMatrix)==1))*(-1),]
}
else{
m2 <- incidenceMatrix
}
m <- t(m2) %*% (m2*(1/(rowSums(m2)-1)))
}
else {
if (length(c(which(rowSums(incidenceMatrix)==1)))!=0){
m2 <- incidenceMatrix[c(which(rowSums(incidenceMatrix)==1))*(-1),]
}
else{
m2 <- incidenceMatrix
}
m <- t(m2) %*% (m2*(1/(rowSums(m2)-1)))
}
if (length(E(proj))>0){
for (i in 1:(length(E(proj)))){
E(proj)[i]$weight <- m[V(proj)[get.edge(proj,i)[[1]]]$name, V(proj)[get.edge(proj,i)[[2]]]$name]
}
}
return(proj)
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.