WGCNA分析生物信息学R语言源码生信分析流程

WGCNA构建共表达网络四

2020-02-24  本文已影响0人  多啦A梦的时光机_648d

一:提取指定模块的基因名

提取基因信息,进行下游分析包括GO/KEGG等功能数据库的注释。

# Select module
module = "turquoise";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule];

二:模块的导出

主要模块里面的基因直接的相互作用关系信息可以导出到cytoscape,VisANT等网络可视化软件。

# Recalculate topological overlap
TOM = TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate); 
# Select module
module = "turquoise";
# Select module probes
probes = colnames(datExpr) ## 我们例子里面的probe就是基因名
inModule = (moduleColors==module);
modProbes = probes[inModule]; 
## 也是提取指定模块的基因名
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)

三:模块对应的基因关系矩阵

1.首先是导出到VisANT

vis = exportNetworkToVisANT(modTOM,
                            file = paste("VisANTInput-", module, ".txt", sep=""),
                            weighted = TRUE,
                            threshold = 0)

2.然后是导出到cytoscape

cyt = exportNetworkToCytoscape(
  modTOM,
  edgeFile = paste("CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
  nodeFile = paste("CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
  weighted = TRUE,
  threshold = 0.02,
  nodeNames = modProbes, 
  nodeAttr = moduleColors[inModule]
);

四:模块内的分析—— 提取hub genes

hub genes指模块中连通性(connectivity)较高的基因,如设定排名前30或前10%)。
高连通性的Hub基因通常为调控因子(调控网络中处于偏上游的位置),而低连通性的基因通常为调控网络中偏下游的基因(例如,转运蛋白、催化酶等)。
HubGene: |kME| >=阈值(0.8)

4.1 计算连通性

### Intramodular connectivity, module membership, and screening for intramodular hub genes

# (1) Intramodular connectivity

# moduleColors <- labels2colors(net$colors)
connet=abs(cor(datExpr,use="p"))^6
Alldegrees1=intramodularConnectivity(connet, moduleColors)
head(Alldegrees1)

4.2 模块内的连通性与gene显著性的关系

# (2) Relationship between gene significance and intramodular connectivity
which.module="blue"
EB= as.data.frame(datTraits[,1]); # change specific 
names(EB) = "EB"
GS1 = as.numeric(cor(EB,datExpr, use="p"))
GeneSignificance=abs(GS1)
colorlevels=unique(moduleColors)
sizeGrWindow(9,6)
par(mfrow=c(2,as.integer(0.5+length(colorlevels)/2)))
par(mar = c(4,5,3,1))
for (i in c(1:length(colorlevels)))
{
  whichmodule=colorlevels[[i]];
  restrict1 = (moduleColors==whichmodule);
  verboseScatterplot(Alldegrees1$kWithin[restrict1],
                     GeneSignificance[restrict1], col=moduleColors[restrict1],
                     main=whichmodule,
                     xlab = "Connectivity", ylab = "Gene Significance", abline = TRUE)
}
结果

4.3 计算模块内所有基因的连通性, 筛选hub genes

abs(GS1)> .9 可以根据实际情况调整参数
abs(datKME$MM.black)>.8 至少大于 >0.8

###(3) Generalizing intramodular connectivity for all genes on the array
datKME=signedKME(datExpr, MEs, outputColumnName="MM.")
# Display the first few rows of the data frame
head(datKME)
##Finding genes with high gene significance and high intramodular connectivity in
# interesting modules
# abs(GS1)> .9 可以根据实际情况调整参数
# abs(datKME$MM.black)>.8 至少大于 >0.8
FilterGenes= abs(GS1)> .9 & abs(datKME$MM.black)>.8
table(FilterGenes)

4.4 another plot for realtionship between module eigengenes

displaying module heatmap and the eigengene

sizeGrWindow(8,7);
which.module="blue"
ME=MEs[, paste("ME",which.module, sep="")]
par(mfrow=c(2,1), mar=c(0.3, 5.5, 3, 2))
plotMat(t(scale(datExpr[,moduleColors==which.module ]) ),
        nrgcols=30,rlabels=F,rcols=which.module,
        main=which.module, cex.main=2)
par(mar=c(5, 4.2, 0, 0.7))
barplot(ME, col=which.module, main="", cex.main=2,
        ylab="eigengene expression",xlab="MPP")
结果
上一篇 下一篇

猜你喜欢

热点阅读