WGCNA：（加权共表达网络分析）
文章目录
- 1 介绍
- 2 名词解释
- Co-expression network
- Module
- Connectivity
- Intramodular connectivity kIMk_{IM}kIM
- Module eigengene E
- Eigengene significance
- Module Membership kMEk_{ME}kME
- Hub gene
- Gene significance GS
- Module significance
- 3 流程
- 3.1 数据输入、清洗、预处理
- 3.1.1 输入表达数据
- 3.1.2 检查基因和样品是否有过多缺失值，以及异常样品
- 3.1.3 输入临床特征数据
- 3.2 一步构建共表达网络和划分模块
- 3.2.1 筛选软阀值 β\betaβ
- 3.2.2 构建网络和划分模块
- 3.3 模块与性状关联分析
- 3.3.1 量化模块特征关联
- 3.3.2 基因与性状和重要模块的关系：基因重要性和模块成员
- 3.3.3 模内分析：鉴定具有高GS和MM的基因
- 4 总结
- 参考资料
1 介绍
WGCNA相比于差异表达基因可以获得跟多信息，通过考虑测得的转录本之间的关系可以更完整地表示微列阵数据，这可以通过基因表达谱之间的成对相关性进行评估。WGCNA从数千个基因的层次开始，确定临床上感兴趣的基因模块，最后使用模块内连接性，基因重要性来识别疾病途径中的关键基因，以进行进一步验证。与其将成千上万的基因与微阵列样品特征相关联，不如着重于几个（通常少于10个）模块与样品特征之间的关系。为此，它计算每个模块的特征基因显着性（样本特征与特征基因之间的相关性）和相应的p值。 模块定义不使用先验定义的基因集。取而代之的是，通过使用层次聚类从表达式数据构造模块。 尽管建议将生成的模块与基因本体信息相关联以评估其生物学合理性，但这不是必需的。由于模块可能对应于生物学途径，因此将分析重点放在模块内Hub基因（或模块特征基因）上就等于一种出于生物学目的数据降维。由于模块内hub基因的表达谱高度相关，通常会产生数十种候选生物标记。尽管这些候选者在统计学上是等效的，但它们可能在生物学上的合理性或临床效用之间存在差异。基因本体信息可用于进一步确定模块内Hub基因的优先级。图一显示了WGCNA分析流程。
- 构建共表达网络
- 划分模块
- 模块与性状关联分析
- 模块之间关联分析
- 模块中核心基因鉴定
2 名词解释
Co-expression network
将共表达网络定义为无方向的加权基因网络。 这种网络的节点对应于基因表达谱，基因之间的连线由基因表达之间的成对相关性决定。通过将相关系数的绝对值提高到幂 β>1\beta >1β>1 （软阈值）， 加权基因共表达网络的构建强调了高相关性，却以低相关性为代价。具体来说 aij=∣cor(xi,xj)∣βa_{ij}=|cor(x_i,x_j)|^\betaaij=∣cor(xi,xj)∣β 表示无符号（unsigned）网络的连接，另外，也可以指定有符号的加权共表达网络，如 aij=∣(1+cor(xi,xj))/2∣βa_{ij}=|(1+cor(x_i,x_j))/2|^\betaaij=∣(1+cor(xi,xj))/2∣β。
Module
模块是高度互连的基因簇。 在无符号共表达网络中，模块对应于具有高度绝对相关性的基因簇。 在有符号网络中，模块对应于正相关的基因。
Connectivity
对于每个基因，连通性定义为与其他网络基因的连接强度之和 ki=Σμ≠iaμik_i=\Sigma_{\mu \neq i}a_{\mu i}ki=Σμ=iaμi，在共表达网络中，连通性可衡量基因与所有其他网络基因之间的相关性。
Intramodular connectivity kIMk_{IM}kIM
模内连通性测量给定基因相对于特定模块的基因如何连接或共表达。 模内连通性可以解释为模块成员资格的量度。
Module eigengene E
模块特征向量E被定义为给定模块的主成分一。 可以认为是模块中基因表达谱的代表。
Eigengene significance
当可获得微阵列样品特征y（例如病例对照状态或体重）时，可以将模块特征基因与此结果相关联。 相关系数称为特征基因显著性
Module Membership kMEk_{ME}kME
特征基因连通性，对于每个基因，我们通过将其基因表达谱与给定模块的模块本征基因相关联来定义模块成员的“模糊”度量。如，MMblue(i)=Kcor,iblue=cor(xi,Eblue)MM^{blue}(i)=K^{blue}_{cor,i}=cor(x_i,E^{blue})MMblue(i)=Kcor,iblue=cor(xi,Eblue) 用来测量基因 iii 与蓝色模块特征基因的相关程度。如果 MMblue(i)MM^{blue}(i)MMblue(i) 接近0，则第 iii 个基因不属于blue模块的一部分。 另一方面，如果 MMblue(i)MM^{blue}(i)MMblue(i) 接近于1或 -1，则它与蓝色模块基因高度相关。模块成员的符号编码该基因与蓝色模块本征基因是正相关还是负相关， 可以为所有输入基因定义模块成员度量（无论其原始模块成员如何）。 事实证明，模块成员资格度量与模块内连接性 kIMk_{IM}kIM 高度相关。 高度连接的模块内Hub基因倾向于对各个模块具有较高的模块成员值。
Hub gene
这个松散定义的术语被用作“高度连接的基因”的缩写。通过定义，共表达模块内部的基因往往具有高度的连通性。
Gene significance GS
为了将外部信息整合到共表达网络中，我们利用了基因显著性方法。抽象地说，GSiGS_iGSi 的绝对值越高，第 iii 个基因的生物学意义就越大。例如，GSiGS_iGSi 可以编码通路成员（例如，如果该基因是已知的凋亡基因，则为1，否则为0），敲除必需性或与外部微阵列样品性状的相关性。基因显着性度量也可以通过减去p值的对数来定义。 唯一的要求是，基因显着性0表示该基因对于所关注的生物学问题不重要。 基因显著性可以取正值或负值。
Module significance
模块显著性被定义为给定模块中所有基因的平均绝对基因显著性的度量。 当将基因显着性定义为基因表达与外部性状y的相关性时，此度量往往与模块特征基因与y的相关性高度相关。
3 流程
3.1 数据输入、清洗、预处理
表达数据和性状数据下载
3.1.1 输入表达数据
library(WGCNA)
options(stringsAsFactors = FALSE)
femData = read.csv("LiverFemale3600.csv")
datExpr0 = as.data.frame(t(femData[, -c(1:8)]))
names(datExpr0) = femData$substanceBXH
3.1.2 检查基因和样品是否有过多缺失值，以及异常样品
检测缺失值
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK
[1] TRUE
如果输出语句为TRUE，不需要删除基因和样本数据；如果没有返回TURE，需要从数据中删除潜在的基因和样本。
if (!gsg$allOK) {if (sum(!gsg$goodGenes)>0)printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))if (sum(!gsg$goodSamples)>0)printFlush(paste("Removing samples:", paste(names(datExpr0)[!gsg$goodSamples], collapse = ", ")))datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}
对样本进行聚类，看是否有异常样本
sampleTree = hclust(dist(datExpr0), method = "average")
sizeGrWindow(12,9)
par(cex = 0.6)
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
abline(h = 15, col = 'red')
从上图可以看出，样品F2_221样品疑似异样，考虑去除。
clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust)
clust0 1 1 134
去除异常样品后的表达矩阵。
datExpr = datExpr0[clust == 1, ]
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
3.1.3 输入临床特征数据
traitData = read.csv("ClinicalTraits.csv")
dim(traitData)
[1] 361 38
names(traitData)
[1] "X" "Mice" "Number" "Mouse_ID" [5] "Strain" "sex" "DOB" "parents" [9] "Western_Diet" "Sac_Date" "weight_g" "length_cm"
[13] "ab_fat" "other_fat" "total_fat" "comments"
[17] "X100xfat_weight" "Trigly" "Total_Chol" "HDL_Chol"
[21] "UC" "FFA" "Glucose" "LDL_plus_VLDL"
[25] "MCP_1_phys" "Insulin_ug_l" "Glucose_Insulin" "Leptin_pg_ml"
[29] "Adiponectin" "Aortic.lesions" "Note" "Aneurysm"
[33] "Aortic_cal_M" "Aortic_cal_L" "CoronaryArtery_Cal" "Myocardial_cal"
[37] "BMD_all_limbs" "BMD_femurs_only"
去除不需要的临床特征和样品
allTraits = traitData[, -c(31, 16)]
allTraits = allTraits[, c(2, 11:36)]
datTraits = allTraits[match(rownames(datExpr),allTraits$Mice),]
rownames(datTraits) = datTraits[,1]
datTraits = datTraits[,-1]
对表达矩阵再次聚类
sampleTree2 = hclust(dist(datExpr), method = "average")
将临床特征大小用颜色深浅表示：白色表示数值低，红色表示数值高，灰色表示缺少输入
traitColors = numbers2colors(datTraits, signed = FALSE)
样本聚类热图
plotDendroAndColors(sampleTree2, traitColors, groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap")
3.2 一步构建共表达网络和划分模块
3.2.1 筛选软阀值 β\betaβ
设置一个软阀值向量
powers = c(c(1:10), seq(from = 12, to=20, by=2))
使用网络拓扑分析函数
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
软阈值 β\betaβ 与无尺度网络评价系数 R2R^2R2 的关系，以及软阈值 β\betaβ 与平均连通性的关系
sizeGrWindow(9, 5)
par(mfrow = c(1,2))
cex1 = 0.9
plot(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=cex1,col="red")
abline(h=0.90,col="red") plot(sft$fitIndices[,1], sft$fitIndices[,5], xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n", main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
R2R^2R2 为无尺度网络评价系数，一般设置为0.9，β\betaβ 取值标准：R2R^2R2 第一次到达0.9时对应的 β\betaβ 值，由上图可知，β\betaβ 选为6
3.2.2 构建网络和划分模块
构建共表达网络、划分模块、合并相似模块。
net = blockwiseModules(datExpr,power = 6, # 软阀值选为6TOMType = "unsigned", # 构建无尺度网络minModuleSize = 30, # 最小模块基因数为30reassignThreshold = 0, mergeCutHeight = 0.25, # 模块合并阀值numericLabels = F, # 模块颜色标签pamRespectsDendro = FALSE, saveTOMs = TRUE, # 保存TOM矩阵saveTOMFileBase = "femaleMouseTOM",verbose = 3,maxBlockSize = 5000) # 可以处理的数据集的最大基因数，默认5000
# 所有模块个数和各个模块中基因数量
table(net$colors)
black blue brown cyan green greenyellow grey 211 460 409 77 312 100 99 grey60 lightcyan lightgreen magenta midnightblue pink purple 47 58 34 123 76 157 106 red salmon tan turquoise yellow 221 91 94 609 316
结果显示有18个模块，模块大小为34至609个基因，模块grey表示所有模块外部的基因。
moduleColors = net$colors # 模块颜色标签
MEs = net$MEs # 模块特征向量
TOM矩阵层次聚类结果，其中dissTom = 1-Tom，后面进一步介绍TOM矩阵。
geneTree = net$dendrograms[[1]]
geneTree
Call:
fastcluster::hclust(d = as.dist(dissTom), method = "average")
Cluster method : average
Number of objects: 3600
TOM矩阵的层次聚类以及基因所属模块可视化
sizeGrWindow(12, 9)
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]],moduleColors,"Module colors",dendroLabels = FALSE,hang = 0.03,addGuide = TRUE,guideHang = 0.05)
TOM矩阵介绍：
首先由表达矩阵计算邻接矩阵 aij=∣cor(xi,xj)∣βa_{ij}=|cor(x_i,x_j)|^\betaaij=∣cor(xi,xj)∣β ，xix_ixi 和 xjx_jxj表示任意两个基因 ，为了最大程度地减少噪声和虚假关联的影响，需要将邻接矩阵转换为拓扑重叠矩阵（TOM矩阵），常用TOMsimilarity()函数将邻接矩阵转为TOM矩阵。
3.3 模块与性状关联分析
3.3.1 量化模块特征关联
计算模块特征向量和临床性状之间相关系数矩阵，并对相关系数进行检验。
MEs = orderMEs(MEs)
moduleTraitCor = cor(MEs, datTraits, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
可视化
pdf(file = 'labeledHeatmap.pdf',width = 12,height = 8)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits),yLabels = names(MEs), ySymbols = names(MEs),colorLabels = FALSE, colors = greenWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.7, zlim = c(-1,1), cex.lab.y = 0.6,cex.lab.x = 0.7,yColorWidth = strheight("M")/2,main = paste("Module-trait relationships"))
dev.off()
3.3.2 基因与性状和重要模块的关系：基因重要性和模块成员
通过将基因显着性GS定义为基因和性状之间的相关性的绝对值，我们可以量化单个基因与我们感兴趣的性状（体重）的关联。 对于每个模块，我们还定义了模块成员MM的定量度量，作为模块特征基因与基因表达谱的相关性。 这样能够量化微阵列上所有基因与每个模块的相似性。
weight = as.data.frame(datTraits$weight_g)
names(weight) = "weight"
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", names(weight), sep="")
names(GSPvalue) = paste("p.GS.", names(weight), sep="")
3.3.3 模内分析：鉴定具有高GS和MM的基因
module = "brown"
column = match(module, modNames)
moduleGenes = moduleColors==module
table(moduleGenes)
sizeGrWindow(7, 7)
par(mfrow = c(1,1))
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),abs(geneTraitSignificance[moduleGenes, 1]),xlab = paste("Module Membership in", module, "module"),ylab = "Gene significance for body weight",main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
4 总结
加权共表达网络最主要的两个步骤：
- 共表达网络构建和模块划分
- 模块和性状关联分析
通过上面两个步骤找到与目标性状显着相关的模块，对模块里面的基因进行后续分析，如结合差异表达分析。
参考资料
WGCNA参考手册
本博客内容将同步更新到个人微信公众号：生信玩家。欢迎大家关注~~~
WGCNA：（加权共表达网络分析）相关推荐
- WGCNA加权基因共表达网络分析（1）简介、原理
WGCNA简介 WGCNA(Weighted Gene Co-Expression Network Analysis, 加权基因共表达网络分析),鉴定表达模式相似的基因集合(module),解析基因集 ...
- 基因共表达网络分析java,WGCNA:加权基因共表达网络分析
加权基因表达网络分析(Weighted gene co-expression network analysis, WGCNA),又叫权重基因共表达网络分析,其根本思想是根据基因表达模式的不同,挖掘出相 ...
- WGCNA（加权基因共表达网络分析）
WGCNA(加权基因共表达网络分析) 序章 这个工具现在很火,高分文章用到很多. 加权基因共表达网络分析(WGCNA,Weighted gene co-expression network analy ...
- RNA 13. SCI 文章中加权基因共表达网络分析之 WGCNA
WGCNA 分析流程 2008 年发表在 BMC 之后的影响力还是很高的,先后在各大期刊都能看到,但是就其分析的过程来看,还是需要有一定 R 语言的基础才能完整的复现出来文章中的结果,这期就搞出来供大 ...
- 基因共表达网络分析java,基因共表达——基因共表达网络分析
Gene co-expression(基因共表达)是一种使用大量基因表达数据构建基因间的相关性,从而挖掘基因功能的一类分析方法. 在很多情况下,有着相似行为/变化的物质,会存在着一定的联系.在生物中, ...
- 基因共表达网络分析java,好用的基因共表达网络分析工具
原标题:好用的基因共表达网络分析工具 基因共表达网络(GeneCo-expreesion Network)是用来展现基因间相互作用关系的一种手段,是基于基因间表达数据而构建调控网络图.今天推荐一个查询 ...
- 基因共表达网络分析java,基因共表达网络分析-WGCNA
今天推荐给大家一个R包WGCNA,针对我们的表达谱数据进行分析. 简单介绍:WGCNA首先假定基因网络服从无尺度分布,并定义基因共表达相关矩阵.基因网络形成的邻接函数,然后计算不同节点的相异系数,并据 ...
- 基因共表达网络分析java,RNA-seq数据的基因共表达网络分析
目录 1 背景 BackgroundTypes of biological networks Motivation for using co-expression networks Network i ...
- WGCNA将共表达基因与表型数据相关联
欢迎关注微信公众号<生信修炼手册>! 单纯的共表达基因集合的结果并不能与我们的实验设计相关联,对于识别到的几十个共表达基因集合,一一进行富集分析去挖掘其功能,看上去如此的盲目,没有目的性, ...
- 基因共表达网络分析java,WGCNA 加权基因共表达网络分析教程（1）
学是很多做科研的同学都想学的,包括我在内,现阶段正在学习这个,深夜整理材料不易,请多多关照支持! 先安装Rstudio以及以及加载WGCNA数据包 Rstudio软件下载地址: WGCNA加载按照网站 ...
最新文章
- 【系统缓慢、CPU 100%、频繁Full GC问题】的定位排查思路！
- 如何在自己开发的日程管理页面插入提醒功能_微信中6个藏得很深但却很有用的功能...
- html外链式css运行不出来div,html+css外链式
- linux shell命令行及脚本编程实例详解_Linux高手必看的10本经典书籍
- 计算机水平考试改革,浅析全国计算机等级考试改革及应对策略
- 无限极评论怎么删除php,TP5 无限极评论回复
- C# datagridview 删除行（转 学会、放弃博客）
- 消格子时一个很深的bug的修复纪录
- 爆款 | Medium上6900个赞的AI学习路线图，让你快速上手机器学习
- springboot+vue网络课程教学网站系统java源码介绍
- 十八个经典问答，讲透了RS485接口！-小白收藏
- 基于vue的手机阅读小说类webapp
- 数据安全技术专利态势分析
- 最新国民人均年薪出炉，你有没有拉国家的后腿？
- 什么是迭代（迭代法）
- java中controller层是干嘛的？
- 2016流行这2种色彩！附优秀网页设计案例
- 【转】地心历险记 2：神秘岛 迅雷 下载 地址|神秘岛 高清 下载地址
- OSChina 周五乱弹 —— 到底哪个更重要
- 如何评价B端产品经理的能力