生信分析学习笔记—之初探WGCNA

什么是WGCNA

Weighted gene coexpression network analysis是分析不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集, 并根据基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。简单来说，该分析方法旨在寻找协同表达的基因模块(module)，并探索基因网络与关注的表型之间的关联关系，以及网络中的核心基因。

相比较于 Coexpression Network而言，WGCNA不仅仅说明基因间相关与否，还记录着它们的相关性数值，数值就是基因之间的联系的权重(相关性)。

什么时候用WGCNA

WGCNA适用于复杂的多样本转录组数据，常应用于研究不同器官/组织类型和不同阶段的发育调控、生物和非生物胁迫的不同时间点响应机制等。

1、时间序列样品，两个材料相同处理，可以一起做WGCNA分析，也可以分别做然后比较，相同材料不同处理之间也一样

2、同一物种，不同来源的转录组数据，可以放在一起做WGCNA，也可以分开来比较

3、一般建议15个样品以上进行WGCNA分析有生物学意义，可以是5个时间点三个重复的15个样品

4、表型数据，建议收集可以统计量化的性状数据，可将模块和表型数据关联分析，有助于筛选关键基因模块和基因来解释相关表型

WGCNA分析流程

共表达网络（加权基因网络）：点代表基因，边代表基因表达相关性。加权表示对相关性值进行幂次运算（幂次值=软阈值）。

Module（模块，高度内连的基因集）：在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可对模块进行三个层次的分析：第一层次，功能富集分析查看其功能特征是否与研究目的相符；第二层次，模块与性状进行关联分析，找出与关注性状相关度最高的模块；第三层次，模块与样本进行关联分析，找到样品特异高表达的模块。

connectivity（连接度）：类似于网络中“degree“的概念。每个基因的连接度是与其相连的基因的边属性之和。

module membership：给定基因表达谱与给定模型的eigengene的相关性

module eigengene：给定模型的第一主成分，代表整个模型的基因表达谱。

hub gene：连接度最多或连接多个模块的基因。

adjacency matrix：邻接矩阵，基因与基因之间的加权相关性值构成的矩阵。

TOM：把邻接矩阵转换成拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可用于构建网络或绘制TOM图。

WGCNA实际操作

#样本聚类检查离群值

## 查看是否有离群样品,聚类是针对样本用hclust
sampleTree = hclust(dist(t(fpkm)), method = 'average')
plot(sampleTree, main = 'Sample clustering to detect outliers', sub='', xlab='')

看所有样本的基因表达数据，然后进行聚类：比如是先下降后上升，还是先上升后下降，把表达模式相同的聚类在一起。

#软阈值筛选

为什么要进行软阈值的筛选？要解释这个问题就要谈到无尺度分布，我们的基因网络通常不是随意连接的，它是一个有hub的网络，一个强中介可以连接很多很多的基因，类似于你想找一个人的电话号码，实际上你打114就可以了，这个114的号码就是一个强中介，一个hub。你是没有办法来衡量你和你想要找的那个人的距离远近的，这种就叫作无尺度分布，该分布又叫做幂律分布。

在网络中，节点表示基因，若两个基因间的相关系数大于指定的阈值，则会以边连接以表示两者间可能存在相互作用。这里我们引用了加权的概念，强中介的基因的权重一定就是比较大的。这样就会根据基因间相关性的强弱来选择阈值。这样的阈值就是软阈值。如下图所示，基本在26时，拟合幂律分布的结果是比较好的；同时节点的凭据连通性也趋于稳定了。之后会使用这个参数建立网络。

##选择合适“软阀值（soft thresholding power）”beta
powers = c(c(1:20), seq(from = 12, to=30, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 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')

#要求R2大于0.70。
abline(h=0.70,col='green')
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')

左图纵轴：相关系数的平方，越高说明该网络越逼近无网络尺度的分布。右图纵轴：Connectivity类似于度的概念，基因模块中所有基因邻近函数的均值，每个基因的连接度是与其相连的基因的边属性之和。

# 层级聚类树展示各个模块

基于基因在所有样品中表达量数据计算的两两之间的相关性，将相似表达模式的基因鉴定为一个模块（module）。从而将上千个差异基因归为十几或者几十个基因模块。

动态剪切树法切割模块,选择动态混合切割法，可自行设置模块大小。

#用power=26建立邻接矩阵。
softPower = 26
adjacency = adjacency(datExpr, power = softPower)
#将邻接矩阵转换为拓扑矩阵。
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
#以TOM为基础聚类。
geneTree = hclust(as.dist(dissTOM), method = 'average')
sizeGrWindow( 12,9)
plot(geneTree, xlab='',sub='',main = 'Gene clustering on TOM-based dissimilarity',labels = FALSE, hang = 0.04)

#动态剪切树法切割模块,选择动态混合切割法。
minModuleSize = 30 
dynamicMods = cutreeDynamic(dendro = geneTree,distM = dissTOM,deepSplit = 2,pamRespectsDendro = FALSE,
 minClusterSize = minModuleSize)
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
sizeGrWindow( 10,6)
plotDendroAndColors(geneTree,dynamicColors, 'Dynamic Tree Cut',dendroLabels = FALSE,
 hang = 0.03,addGuide = TRUE,guideHang = 0.05,main = 'Gene dendrogram and module colors')

不同的颜色代表不同的模块，每种模块的表达模式可能略有差别。可能基因都是先上升后下降，但是有的模块是上升很快，下降很慢，总之是可以有不同的模式的。

#合并相似的共表达网络

这个过程相对于把上一步筛选得到的modules再进行聚类,使用函数moduleEigengenes(),确定修剪高度后,使用mergeCloseModules()函数合并modules。接下来需要判断有无非常相似的模块，进行合并操作，根据模块的特征值，计算不同模块的相似性，进行层次聚类。

MEList = moduleEigengenes(datExpr,colors = dynamicColors)
MEs = MEList$eigengenes
MEDiss = 1-cor(MEs)
METree = hclust(as.dist(MEDiss),method = 'average')
sizeGrWindow(7, 6)
plot(METree,main = 'Clustering of module eigengenes',xlab = '',sub = '')
#选择相关系数大于0.75的进行合并。
MEDissThres = 0.25
abline(h=MEDissThres, col = 'red')

选择相关系数大于0.75的进行合并

merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres,verbose = 3)
mergedColors = merge$colors
mergedMEs = merge$newMEs
sizeGrWindow( 12, 9)
plotDendroAndColors(geneTree,cbind(dynamicColors,mergedColors),c('Dynamic tree cut','Merged dynamic'),
 dendroLabels = FALSE,hang=0.03,addGuide = TRUE,guideHang = 0.05)

#共表达模块之间的相关性

根据基因间表达量进行聚类所得到的各模块间的相关性图

#根据基因间表达量进行聚类所得到的各模块间的相关性图
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(MEs)
sizeGrWindow(7, 6) 
plotEigengeneNetworks(MET, 'Eigengene adjacency heatmap',marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)

选择需要导出的模块颜色，保存模块赋值和模块特征基因信息，以供后续进一步分析。

本站是提供个人知识管理的网络存储空间，所有内容均由用户发布，不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息，谨防诈骗。如发现有害或侵权内容，请点击一键举报。