計算網(wǎng)絡(luò)節(jié)點模塊內(nèi)連通度(within modular degree)和模塊間連通度(between modular degree)_網(wǎng)絡(luò)的連通度-CSDN博客

1.為什么要計算網(wǎng)絡(luò)節(jié)點模塊內(nèi)連通度和模塊間連通度
2.網(wǎng)絡(luò)節(jié)點模塊內(nèi)連通度（within modular degree）和模塊間連通度(among modular degree)計算公式
3.計算代碼及注釋
1.為什么要計算網(wǎng)絡(luò)節(jié)點模塊內(nèi)連通度和模塊間連通度
??我們可以通過了解節(jié)點的模塊內(nèi)連通度和模塊間連通度來對節(jié)點的重要性進行度量，具有較低的模塊內(nèi)連通度和模塊間連通度的節(jié)點的重要性則相對較低,具體劃分如下表（其中Zi為節(jié)點i的模塊內(nèi)連通度的z-scores，Pi為根據(jù)節(jié)點i的模塊間連通度計算的參與系數(shù)，participation coefficient）：

定義

計算公式

#該代碼參考自：https://github.com/cwatson/brainGraph/blob/master/R/vertex_roles.R
#g表示graph文件，其中g(shù)raph的節(jié)點含有modular屬性；
#也就是說，使用下面的函數(shù)計算前請先計算各節(jié)點的模塊屬性，并將其賦值給節(jié)點
#wtc <- cluster_louvain(igraph,NA)
#modularity(wtc)
#V(igraph)$module<-membership(wtc)

#計算模塊內(nèi)連接度的z-scores
within_module_deg_z_score <- function(g, A=NULL, weighted=FALSE) {
  stopifnot(is_igraph(g))
  if (is.null(A)) {
    if (isTRUE(weighted)) {
      A <- as_adj(g, sparse=FALSE, names=TRUE, attr='weight')
    } else {
      A <- as_adj(g, sparse=FALSE, names=TRUE)
    }
  }
  memb <- vertex_attr(g, "module")
  N <- max(memb)
  nS <- tabulate(memb)
  z <- Ki <- rep.int(0, dim(A)[1L])
  Ksi <- sigKsi <- rep.int(0, N)
  names(z) <- names(Ki) <- rownames(A)
  for (S in seq_len(N)) {
    x <- rowSums(A[memb == S, memb == S])
    Ki[memb == S] <- x
    Ksi[S] <- sum(x) / nS[S]
    sigKsi[S] <- sqrt(sum((x - Ksi[S])^2) / (nS[S]-1))
  }
  z <- (Ki - Ksi[memb]) / sigKsi[memb]
  z[is.infinite(z)] <- 0
  df <- data.frame(Ki,z,row.names = names(Ki))
  return(df)
}

#計算節(jié)點的參與系數(shù)
part_coeff <- function(g, A=NULL, weighted=FALSE) {
  stopifnot(is_igraph(g))
  if (is.null(A)) {
    if (isTRUE(weighted)) {
      A <- as_adj(g, sparse=FALSE, attr='weight')
    } else {
      A <- as_adj(g, sparse=FALSE)
    }
  }
  memb <- vertex_attr(g, "module")
  Ki <- colSums(A)
  N <- max(memb)
  Kis <- t(rowsum(A, memb))
  pi <- 1 - ((1 / Ki^2) * rowSums(Kis^2))
  names(pi) <- rownames(A)
  return(pi)
}

還有其它博客也介紹了Zi和Pi的計算，例如：
使用ggClusterNet一條代碼計算網(wǎng)絡(luò)模塊內(nèi)連通度（Zi）和模塊間連通度（Pi）
加載了ggClusterNet包后，只需要一行函數(shù)即可計算Zi和Pi，但仔細(xì)閱讀了其計算的函數(shù)之后，對一些地方存在一些疑問（在代碼中注釋了），此時該函數(shù)計算的結(jié)果也與上面函數(shù)within_module_deg_z_score()計算結(jié)果不同，經(jīng)過修改后的函數(shù)within_module_degree2()計算的結(jié)果與上面函數(shù)within_module_deg_z_score()計算的結(jié)果計算結(jié)果相同。
因此，在使用R包時，特別是一些未經(jīng)嚴(yán)格檢驗的R包時，需要了解更多的細(xì)節(jié)，以免計算錯誤，得到錯誤的結(jié)果。

#參考自：https://github.com/taowenmicro/ggClusterNet/blob/master/R/module.roles.R
network_degree <- function(comm_graph){
  ki_total <-NULL
  net_degree <- degree(comm_graph)
  for(i in 1:length(V(comm_graph))){    
    ki <- net_degree[i]    
    tmp <- data.frame(taxa=names(ki), total_links=ki)   
    if(is.null(ki_total)){ki_total<-tmp} else{ki_total <- rbind(ki_total, tmp)}   
  } 
  return(ki_total) 
}

#compute within-module degree for each of the features

within_module_degree <- function(comm_graph){
  mods <- vertex_attr(comm_graph, "module")
  vs <- as.list(V(comm_graph))
  modvs <- data.frame("taxon"= names(vs), "mod"=mods)
  #不理解此處為什么要用decompose.graph(),難道不應(yīng)該直接用模塊內(nèi)的節(jié)點去取子圖嗎？我在within_module_degree2()中對我的想法進行了實現(xiàn)
  sg1 <- decompose.graph(comm_graph,mode="strong")
  df <- data.frame()
  for(mod in unique(modvs$mod)){ 
    mod_nodes <- subset(modvs$taxon,modvs$mod==mod) 
    neighverts <- unique(unlist(sapply(sg1,FUN=function(s){if(any(V(s)$name %in% mod_nodes)) V(s)$name else NULL})))  
    g3 <- induced.subgraph(graph=comm_graph,vids=neighverts)   
    mod_degree <- degree(g3) 
    for(i in mod_nodes){      
      ki <- mod_degree[which(names(mod_degree)==i)]     
      tmp <- data.frame(module=mod, taxa=names(ki), mod_links=ki)      
      df <- rbind(df,tmp)      
    }   
  } 
  return(df) 
}

within_module_degree2 <- function(comm_graph){ 
  mods <- vertex_attr(comm_graph, "module")
  vs <- as.list(V(comm_graph)) 
  modvs <- data.frame("taxon"= names(vs), "mod"=mods)
  df <- data.frame()
  for(mod in unique(modvs$mod)){  
    mod_nodes <- subset(modvs$taxon,modvs$mod==mod) 
    g3 <- induced.subgraph(graph=comm_graph,vids=mod_nodes)  
    mod_degree <- degree(g3) 
    for(i in mod_nodes){      
      ki <- mod_degree[which(names(mod_degree)==i)]      
      tmp <- data.frame(module=mod, taxa=names(ki), mod_links=ki)      
      df <- rbind(df,tmp)    
    }   
  } 
  return(df)
}
#calculate the degree (links) of each node to nodes in other modules.
among_module_connectivity <- function(comm_graph){ 
  mods <- vertex_attr(comm_graph, "module")  
  vs <- as.list(V(comm_graph)) 
  modvs <- data.frame("taxa"= names(vs), "mod"=mods)
  df <- data.frame() 
  for(i in modvs$taxa){   
    for(j in modvs$taxa){     
      if(are_adjacent(graph=comm_graph, v1=i , v2=j)){      
        mod <- subset(modvs$mod, modvs$taxa==j)       
        tmp <- data.frame(taxa=i, taxa2=j, deg=1, mod_links=mod)       
        df <- rbind(df, tmp)       
      }     
    }    
  }  
  out <- aggregate(list(mod_links=df$deg), by=list(taxa=df$taxa, module=df$mod), FUN=sum)  
  return(out)  
}

#compute within-module degree z-score which
#measures how well-connected a node is to other nodes in the module.

zscore <- function(mod.degree){  
  ksi_bar <- aggregate(mod_links ~ module, data=mod.degree, FUN = mean)  
  ksi_sigma <- aggregate(mod_links ~ module, data=mod.degree, FUN = sd) 
  z <- NULL 
  for(i in 1:dim(mod.degree)[1]){   
    mod_mean <- ksi_bar$mod_links[which(ksi_bar$module == mod.degree$module[i])]    
    mod_sig <- ksi_sigma$mod_links[which(ksi_bar$module == mod.degree$module[i])]    
    z[i] <- (mod.degree$mod_links[i] - mod_mean)/mod_sig   
  }  
  z <- data.frame(row.names=rownames(mod.degree), z, module=mod.degree$module)  
  return(z) 
}

#The participation coefficient of a node measures how well a  node is distributed
# in the entire network. It is close to 1 if its links are uniformly
#distributed among all the modules and 0 if all its links are within its own module.

participation_coeffiecient <- function(mod.degree, total.degree){ 
  p <- NULL  
  for(i in total.degree$taxa){   
    ki <- subset(total.degree$total_links, total.degree$taxa==i)   
    taxa.mod.degree <- subset(mod.degree$mod_links, mod.degree$taxa==i)   
    p[i] <- 1 - (sum((taxa.mod.degree)**2)/ki**2)   
  } 
  p <- as.data.frame(p)
  return(p)
}