加权基因共表达网络资料集锦：

1.https://mp.weixin.qq.com/s/LvlBt_gCLEyxrpGkKkkIZA

写在前面：大家黏贴复制代码的时候一定注意格式，有可能一个小短线或者逗号、点之类的格式不对就有可能报错，实在眼看不出来就自己手打一遍，这是小白的我最近跑数据入的最多的坑，大家要注意。

来源：生信修炼手册，我看的第一篇WGCNA的推文，写的很详细，很赞，但是代码对于初学者来说还得再细细体会，于是有了后面我推荐的各种推文。文末附。

WGCNA是目前非常火热的一项研究内容，其全称为weighted correlation network analysis, 直译就是加权基因相关性网络分析。通过这项分析，可以鉴定共表达的基因集合，这样的集合称之为modules, 而且可以将modules与表型数据进行关联分析，挖掘潜在的mark 基因。

这个高大上的分析内容的第一步就是构建基因之间的共表达网络，共表达是常用的分析相关性的一种策略，直接通过线性相关函数来计算相关性，比如pearson， spearman等相关系数，每两个基因之间可以计算出一个相关系数，那么如何构建出相关性网络呢？

在基因的相关性网络中，每个节点代表一个基因，节点之间的连线用来表示两个基因的相关性。在传统的相关性分析中，通常会给定一个阈值，比如相关系数的绝对值必须大于0.9，才认为这两个基因间存在相关性。对应的公式如下

S表示两个基因间的相关系数的绝对值，公式如下

注意是绝对值，因为协同变化的基因可以是正相关，也可以是负相关。给定一个阈值，如果两个基因之间的相关系数大于该阈值，则认为这两个基因存在相关性，在网络图中就用一条线将这两个基因连接起来；如果小于该阈值，则不存在相关性。

通过阈值筛选，将两个基因间的相关系数转换为0和1,0代表没有相关性，1代表有相关性，所有基因之间的关系可以用以下矩阵来表示

geneA geneB geneC

geneA 0 1 1

geneB 1 0 1

geneC 1 0 0

这样的矩阵称之为邻接矩阵，通过这个矩阵可以直观的表示一个网络，数值为1的点对应的两个基因在网络图中有连线。

用上述方法构建出的网络，称之为非加权的共表达网络，对于两个基因而言，其相关性是有强弱的，是一个在0到1 分为内波动的值，采用上述一刀切的方法，缺失了原本的变化趋势，所以非加权的共表达网络丢失了很多信息。

WGCNA的开发团队提出了加权基因共表达网络的概念，怎么加权呢，公式如下

在计算邻接矩阵中两个基因的值时，将原本的相关系数的绝对值做一个乘方运算。乘方运算强化了相关系数的变化层次，比如原本系数相差，乘方运算后其差距会被拉大，这样使得数据可以区分的更开，有利于后续聚类识别modules, 同时乘方运算也保证了相关性关系的不变性，公式如下

取log之后，二者是一个线性关系，采用乘方运算来计算基因间的邻接矩阵，用该矩阵构建共表达网络时，两个基因之间的连线不在是有无的关系，而有对应的数字的加权了，这个就是加权基因共表达网络。

在构建加权基因共表达网络时有几点注意事项，最关键的一点是样本数目，当样本太少时，简单线性相关系数并不能有效识别基因间的相关性，会出现很多基因间的相关系数完全一致的情况，这样的数据就很难进一步挖掘，官方推荐至少20个样本，另外就是基因表达谱的数据的预处理，在计算相关性时，表达量数值很低的基因容易造成干扰，会发现它与非常多的基因都存在相关性，所以可以指定一个阈值，将表达量很低的基因去除。

WGCNA如何从module中挖掘关键基因

识别到与表型数据相关的modules之后，还可以在该modules中进一步筛选基因，为了方便筛选，对于每个基因定义了以下三个统计量

1. connectivity

连接度，在之前的文章中，简单提过这个概念，类似于网络中节点的degree的概念，只不过在加权共表达网络中，由于每条边代表两个基因间的相关性的大小，对应一个数值，所以一个基因在共表达网络中的连接度定义为与该基因相连的所有边的数值之和。

另外，根据相连的基因是否和该基因位于同一个module, 又可以将边分为两类，和该基因位于同一个module内，定义为within，位于不同的modules, 定义为out。在WGCNA中，可以通过intramodularConnectivity函数计算连接度，用法如下

# 计算基因间的邻接值

ADJ1=abs(cor(datExpr,use="p"))^6

#计算连接度

Alldegrees1=intramodularConnectivity(ADJ1, colorh1)

计算的结果如下

> head(Alldegrees1)

kTotal kWithin kOut kDiff

Gene1 31.80186 28.37595 3.425906 24.95005

Gene2 28.88249 26.47896 2.403522 24.07544

Gene3 25.38600 23.11852 2.267486 20.85103

KTotal代表该基因的所有边的连接度，KWithin代表和该基因位于同一个module下的边的连接度，KOut代表和该基因位于不同module下的边的连接度，所以KTotal是KWithin和KOut之和，KDiff代表KWithin和KOut的差值。

在module中，会存在hub gene的概念，所谓的hub gene, 就是该module下连接度最大的基因，注意此时只考虑位于该module下的边，就是上文的KWithin。

2. module member-ship

简称MM, 将该基因的表达量与module的第一主成分，即module eigengene进行相关性分析就可以得到MM值，所以MM值本质上是一个相关系数，如果基因和某个module的MM值为0，说明二者根本不相关，该基因不属于这个module; 如果MM的绝对值接近1，说明基因与该module相关性很高。

在WGCNA中，计算基因与module之间的MM值的代码如下

datKME = signedKME(

datExpr,

datME,

outputColumnName="MM.")

第一个参数为基因表达量，第二个参数为Module Eigengene值，结果如下

> head(datKME)

MM.blue MM.brown MM.green MM.grey MM.turquoise MM.yellow

Gene1 0.6830511 0.11547756 -0.007124794 0.2840109 0.9481457 0.09588170

Gene2 0.6342657 0.02257975 0.080277091 0.3029967 0.9356343 0.06889483

Gene3 -0.6198067 -0.12531203 0.008372054 -0.2776929 -0.9121710 -0.17852211

Gene4 0.5966736 0.06469079 0.049862112 0.2671967 0.9052030 0.11707603

Gene5 0.6642214 0.14369720 -0.017975774 0.2442237 0.9017972 -0.01038067

Gene6 -0.6018161 -0.15167072 0.006667131 -0.2053897 -0.9192597 -0.17138960

3. gene signigicancer

简称GS, 将该基因的表达量与对应的表型数值进行相关性分析，最终的相关系数的值就是GS, GS反映出基因表达量与表型数据的相关性，计算GS的代码如下

GS1=as.numeric(cor(y,datExpr, use="p"))

通过以上三个量化指标，可以方便对module下的基因进行筛选。通常情况下，通过module和表型数据的相关性分析，我们可以筛选得到和感兴趣的某一表型相关的具体的modules，在该module下面深入挖掘基因时，可以通过MM和GS两个指标作为过滤手段，示例如下

FilterGenes= abs(GS1)> .2 & abs(datKME$MM.brown)>.8

假设brown是我们找到的和表型高度相关的module, 其中的关键基因可以定义为和brown这一表型的GS值大于0.2，而且MM值大于0.8的基因。筛选出关键基因后，可以通过功能富集分析进一步挖掘其功能。

WGCNA实战练习

本文采用WGCNA官网的Tutirial 1的数据，对加权基因共表达网络分析和后续的数据挖掘的具体操作进行梳理，数据可以从官网下载，示意图如下

整个分析流程可以分为以下几个步骤

1. 数据预处理

这部分包括以下4个内容

读取基因表达量数据

对样本和基因进行过滤

读取样本表型数据

可视化样本聚类树和表型数据

官方的示例数据是一个小鼠的芯片表达谱数据，包含了135个雌性小鼠的数据，在提供的表达谱数据中，除了探针ID和样本表达量之外，还有额外的探针注释信息，在读取原始数据时，需要把多余注释信息去除，代码如下

# 读取文件

options(stringsAsFactors = FALSE)

femData = read.csv("LiverFemale3600.csv")

# 去除多余的注释信息列

datExpr0 = as.data.frame(t(femData[, -c(1:8)]))

names(datExpr0) = femData$substanceBXH

rownames(datExpr0) = names(femData)[-c(1:8)]

对于基因的表达量数据，需要进行过滤，对于基因而言，可以过滤缺失值或者低表达的基因，对于样本而言，如果该样本中基因缺失值很多，也需要过滤，WGCNA内置了一个检验基因和样本的函数，通过该函数可以进行一个基本过滤，代码如下

gsg = goodSamplesGenes(datExpr0)

if (!gsg$allOK) {

datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]

}

goodSamples和goodGenes就是需要保留的基因和样本。基础过滤之后，还可以看下是否存在离群值的样本，通过样本的聚类树进行判断，代码如下

plot(sampleTree,

main = "Sample clustering to detect outliers",

sub="", xlab="", cex.lab = 1.5,

cex.axis = 1.5, cex.main = 2)

生成的图片如下

从图上可以看出，F2_221 这个样本和其他样本差距很大，可以将该样本过滤掉。代码如下

clust = cutreeStatic(

sampleTree,

cutHeight = 15,

minSize = 10)

keepSamples = (clust==1)

datExpr = datExpr0[keepSamples, ]

nGenes = ncol(datExpr)

nSamples = nrow(datExpr)

表型数据中也包含了不需要的列，而且其样本比表达谱的样本多，需要根据表达谱的样本提取对应的表型数据，代码如下

# 读取文件

traitData = read.csv("ClinicalTraits.csv")

# 删除多余的列

allTraits = traitData[, -c(31, 16)]

allTraits = allTraits[, c(2, 11:36) ]

# 样本和表达谱的样本保持一致

femaleSamples = rownames(datExpr)

traitRows = match(femaleSamples, allTraits$Mice)

datTraits = allTraits[traitRows, -1]

rownames(datTraits) = allTraits[traitRows, 1]

表达谱数据和表型数据准备好之后，可以绘制样本聚类树和表型的热图，代码如下

# 由于去除了样本，重新对剩余样本聚类

sampleTree2 = hclust(dist(datExpr), method = "average")

traitColors = numbers2colors(datTraits, signed = FALSE)

plotDendroAndColors(

sampleTree2,

traitColors,

groupLabels = names(datTraits),

main = "Sample dendrogram and trait heatmap")

生成的图片如下

上半部分为样本的聚类树，下班部分为样本对应的表型的热图，顺序和聚类树中的顺序一致，表达量从低到高，颜色从白色过渡到红色，灰色代表缺失值。

2. 构建共表达网络，识别modules

在构建共表达网络时，将基因间的相关系数进行乘方运算来表征其相关性，首先需要确定乘方的值，代码如下

# 设定一些列power梯度

powers = c(c(1:10), seq(from = 12, to=20, by=2))

sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

在sft这个对象中保存了每个power值计算出来的网络的特征，结构如下

> str(sft)

List of 2

$ powerEstimate: num 6

$ fitIndices :'data.frame': 15 obs. of 7 variables:

..$ Power : num [1:15] 1 2 3 4 5 6 7 8 9 10 ...

..$ SFT.R.sq : num [1:15] 0.0278 0.1264 0.3404 0.5062 0.6807 ...

..$ slope : num [1:15] 0.345 -0.597 -1.03 -1.422 -1.716 ...

..$ truncated.R.sq: num [1:15] 0.456 0.843 0.972 0.973 0.94 ...

..$ mean.k. : num [1:15] 747 254.5 111 56.5 32.2 ...

..$ median.k. : num [1:15] 761.7 250.8 101.7 47.2 25.1 ...

..$ max.k. : num [1:15] 1206 574 324 202 134 ...

其中powerEstimate就是最佳的power值，fitIndices保存了每个power对应的网络的特征。我们可以对不同power取值下，网络的特征进行可视化，代码如下

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=0.9,

col="red"

)

abline(h=0.90, col="red")

生成的图片如下

sft$fitIndices保存了每个power构建的相关性网络中的连接度的统计值，k就是连接度值，可以看到，对于每个power值，提供了max, median, max3种连接度的统计量，这里对连接度的均值进行可视化，代码如下

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"

)

生成的图片如下

确定好power值之后，可以直接构建相关性网络

net = blockwiseModules(

datExpr,

power = sft$powerEstimate,

TOMType = "unsigned",

minModuleSize = 30,

reassignThreshold = 0,

mergeCutHeight = 0.25,

numericLabels = TRUE,

pamRespectsDendro = FALSE,

saveTOMs = TRUE,

saveTOMFileBase = "femaleMouseTOM",

verbose = 3)

net对象中保存了所有相关性网络和module的结果，可以将基因的聚类树和对应的module进行可视化，代码如下

mergedColors = labels2colors(net$colors)

plotDendroAndColors(

net$dendrograms[[1]],

mergedColors[net$blockGenes[[1]]],

"Module colors",

dendroLabels = FALSE,

hang = 0.03,

addGuide = TRUE,

guideHang = 0.05

)

生成的图片如下

上方为基因的聚类树，聚类时的距离为1-TOM值，下方为基因对应的modules。类似的，还可以结合基因间的距离，即1-TOM值，用热图展示，代码如下

geneTree = net$dendrograms[[1]]

moduleColors = labels2colors(net$colors)

dissTOM = 1 - TOMsimilarityFromExpr(

datExpr,

power = sft$powerEstimate)

plotTOM = dissTOM ^ 7

diag(plotTOM) = NA

TOMplot(

plotTOM,

geneTree,

moduleColors,

main = "Network heatmap plot, all genes"

)

生成的图片如下

在前面文章我们提到过，在识别module的过程中共，首先用dynamicTreeCut识别modules, 然后根据Module eigengene间的相关性合并modules,net`这个对象中保存了合并前和合并后的modules, 可以将二者画在同一张图上，可视化代码如下

unmergedColors = labels2colors(net$unmergedColors)

mergedColors = labels2colors(net$colors)

plotDendroAndColors(

net$dendrograms[[1]],

cbind(unmergedColors[net$blockGenes[[1]]], mergedColors[net$blockGenes[[1]]]),

c("Dynamic Tree Cut" , "Merged colors"),

dendroLabels = FALSE,

hang = 0.03,

addGuide = TRUE,

guideHang = 0.05

)

生成的图片如下

对于合并前的modules, 其相关性分析的结果可视化如下

unmergedColors = labels2colors(net$unmergedColors)

MEList = moduleEigengenes(datExpr, colors = unmergedColors)

MEs = MEList$eigengenes

MEDiss = 1-cor(MEs)

METree = hclust(as.dist(MEDiss), method = "average")

plot(METree,

main = "Clustering of module eigengenes",

xlab = "",

sub = "")

生成的图片如下

对于每个module而言，我们希望知道该module下对应的基因，提取方式如下

> moduleColors = labels2colors(net$colors)

> unique(moduleColors)

[1] "grey" "turquoise" "grey60" "yellow" "tan"

[6] "green" "red" "black" "blue" "midnightblue"

[11] "cyan" "magenta" "salmon" "lightgreen" "brown"

[16] "purple" "pink" "greenyellow" "lightcyan"

> head(names(datExpr)[moduleColors=="red"])

[1] "MMT00000159" "MMT00000793" "MMT00000840" "MMT00001154" "MMT00001245"

[6] "MMT00001260"

同样我们也可以提取module对应的基因表达量数据，绘制热图, 代码如下

which.module="red"

plotMat(

t(scale(datExpr[,moduleColors==which.module ]) ),

nrgcols=30,

rlabels=F,

rcols=which.module,

main=which.module,

cex.main=2

)

生成的图片如下

3. 筛选与表型相关的modules

本质上是计算module的ME值与表型的相关系数，代码如下

nGenes = ncol(datExpr)

nSamples = nrow(datExpr)

MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes

MEs = orderMEs(MEs0)

moduleTraitCor = cor(

MEs,

datTraits,

use = "p"

)

moduleTraitPvalue = corPvalueStudent(

moduleTraitCor,

nSamples

)

可以对module和表型间的系数的结果进行可视化，代码如下

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 = blueWhiteRed(50),

textMatrix = textMatrix,

setStdMargins = FALSE,

cex.text = 0.5,

zlim = c(-1,1),

main = paste("Module-trait relationships")

)

生成的图片如下

指定一个我们感兴趣的表型，可以得到与其相关性最高的module, 代码如下

> which.trait <- "weight_g"

> moduleTraitCor[, which.trait]

> moduleTraitCor[, which.trait]

MEmagenta MEblack MEturquoise MEgreen MElightcyan

-0.017418109 -0.312679561 -0.272907078 0.001339804 -0.128053858

MEblue MEbrown MEred MEsalmon MEyellow

0.314323101 0.591340840 0.509942529 0.432058666 0.219900538

MElightgreen MEgreenyellow MEgrey60 MEpink MEpurple

-0.057215182 -0.022394396 -0.016705204 -0.051495573 -0.021167541

MEtan MEcyan MEmidnightblue MEgrey

0.269827166 0.181595161 0.193569095 0.089702947

以上结果中，和weight_g最相关的为module为MEred，当然也可以自己指定一个阈值，筛选出多个候选的modules。在WGCNA中，对于基因定义了GS值，表征基因和表型之间的相关性，对于module而言，也可以用所有基因GS绝对值的平均数来表征该module与表型之间的相关性，代码如下

moduleColors = labels2colors(net$colors)

which.trait <- "weight_g"

y <- datTraits[, which.trait]

GS <- as.numeric(cor(y ,datExpr, use="p"))

GeneSignificance <- abs(GS)

ModuleSignificance <- tapply(

GeneSignificance,

moduleColors, mean, na.rm=T)

plotModuleSignificance(GeneSignificance, moduleColors)

生成的图片如下

可以看到brown, red这两个模块和体重相关。对于ME和某一表型而言，还可以将数据合并，聚类展示，代码如下

weight <- datTraits[, which.trait]

MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes

MEs = orderMEs(MEs0)

MET = orderMEs(cbind(MEs, weight))

par(mar = c(4, 2, 1, 2), cex = 0.9)

plotEigengeneNetworks(

MET, "",

marDendro = c(0,4,1,2),

marHeatmap = c(3,4,1,2),

cex.lab = 0.8,

xLabelsAngle = 90

)

生成的图片如下

4. 筛选关键基因

筛选出与表型高相关的modules之后，还可以对modules下的基因进行进一步筛选，主要根据GS值和MM值，代码如下

datKME = signedKME(

datExpr,

MEs,

outputColumnName="MM.")

FilterGenes= abs(GS1)> .2 & abs(datKME$MM.brown)>.8

筛选出候选基因后，可以进行下游的功能富集分析，使用clusterProfiler等R包，进一步挖掘功能，之前的文章有详细介绍如何进行富集分析，这里就不展开了。

5. 导出module数据, 绘制网络图

可以导出指定modules对应的基因共表达网络，方便可视化，代码如下

TOM = TOMsimilarityFromExpr(datExpr, power = 6)

modules = c("brown", "red")

probes = names(datExpr)

inModule = is.finite(match(moduleColors, modules));

modProbes = probes[inModule];

modTOM = TOM[inModule, inModule];

dimnames(modTOM) = list(modProbes, modProbes)

cyt = exportNetworkToCytoscape(

modTOM,

edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),

nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),

weighted = TRUE,

threshold = 0.02,

nodeNames = modProbes,

nodeAttr = moduleColors[inModule]

)