替代WGCNA？初试MEGENA分析
提到共表达网络，大家用多了WGCNA，今天我们来试试用MEGENA来做共表达网络分析。
首先自然是安装加载R包。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("MEGENA")
library(MEGENA)然后提前设置参数
n.cores <- 2; # 要调用PCP的内核/线程数
doPar <-TRUE; # 是否要设置并行
method = "pearson"; # 相关性计算方法 pearson or spearman
FDR.cutoff = 0.05; # 定义FDR阈值
module.pval = 0.05; # 假定模块有意义的阈值。推荐的是0.05。
hub.pval = 0.05; # 由 random tetrahedral networks得到的连接度有意义的阈值
cor.perm = 10; # 用于计算所有相关对的FDRs
hub.perm = 100; # 用于计算连通显著性p值
# 要在下游执行的注释
annot.table=NULL;
id.col = 1;
symbol.col= 2;接下来我们将使用MEGENA包中自带的示例数据给大家展示。大家用自己的基因表达量矩阵也可
data(Sample_Expression) # 加载示例数据
datExpr[1:4,1:4]
ijw <- calculate.correlation(datExpr,doPerm = cor.perm,output.corTable = FALSE,output.permFDR = FALSE)在下面一步中，我们将使用平面滤波网络(PFN)计算有显著相关性的基因对。在使用不同相似度度量的情况下，可以独立地将结果格式化为3列的数据帧，列名为c(row, col, weight)，并确保权重列在0到1之间。
run.par = doPar & (getDoParWorkers() == 1)
if (run.par)
{
cl <- parallel::makeCluster(n.cores)
registerDoParallel(cl)
cat(paste("number of cores to use:",getDoParWorkers(),"\n",sep = ""))
}# 计算PFN，筛选有意义基因对
el <- calculate.PFN(ijw[,1:3],doPar = doPar,num.cores = n.cores,keep.track = FALSE)
g <- graph.data.frame(el,directed = FALSE)通过MCA聚类进行多尺度聚类分析。“MEGENA.output“”是用于下游分析总结和绘图的核心输出。
# 进行MCA聚类
MEGENA.output <- do.MEGENA(g,
mod.pval = module.pval,hub.pval = hub.pval,remove.unsig = TRUE,
min.size = 10,max.size = vcount(g)/2,
doPar = doPar,num.cores = n.cores,n.perm = hub.perm,
save.output = FALSE)# 结果总结
summary.output <- MEGENA.ModuleSummary(MEGENA.output,
mod.pvalue = module.pval,hub.pvalue = hub.pval,
min.size = 10,max.size = vcount(g)/2,
annot.table = annot.table,id.col = id.col,symbol.col = symbol.col,
output.sig = TRUE)
if (!is.null(annot.table))
{
V(g)$name <- paste(annot.table[[symbol.col]][match(V(g)$name,annot.table[[id.col]])],V(g)$name,sep = "|")
summary.output <- output[c("mapped.modules","module.table")]
names(summary.output)[1] <- "modules"
}
print)
print(summary.output$module.table)
summary.output$module.table# 画模块层次结构
module.table <- summary.output$module.table
colnames(module.table)[1] <- "id" # first column of module table must be labelled as "id".
hierarchy.obj <- plot_module_hierarchy(module.table = module.table,label.scaleFactor = 0.15,
arrow.size = 0.03,node.label.color = "blue")
print
# 画图
# 可以生成细化的模块网络图:
pnet.obj <- plot_module(output.summary = summary.output,PFN = g,subset.module = "c1_3",
layout = "kamada.kawai",label.hubs.only = TRUE,
gene.set = NULL,color.code = "grey",
output.plot = FALSE,out.dir = "modulePlot",col.names = c("magenta","green","cyan"),label.scaleFactor = 20,
hubLabel.col = "black",hubLabel.sizeProp = 1,show.topn.hubs = Inf,show.legend = TRUE)
print(pnet.obj[[1]])到此MEGENA分析就做完了，我们得到了有意义的基因对，以及这些基因分别属于哪个模块之中，为了复现“Network models of primary melanoma microenvironments identify key melanoma regulators underlying prognosis”文章中的图，我们将这些数据导出，试图通过ggplot2、Cytoscape、igraph等进行展现，最终得了下图，算是稍稍可以将就看了吧。
替代WGCNA？初试MEGENA分析
相关文章