wgcna WGCNA分析，简单全面的最新教程

WGCNA基本概念

加权相关网络分析(WGCNA)是一种用于描述不同样本间基因关联模式的系统生物学方法，可用于识别高度协同的基因集，并根据基因集的互连性和基因集与表型的相关性识别候选生物标记基因或治疗靶点。

与只关注差异表达基因相比，WGCNA利用数千或近万个变化最大的基因或全部基因的信息来识别感兴趣的基因集，并与表型进行显著相关分析。一是充分利用信息；第二，将数千个基因和表型之间的关联转化为几个基因集和表型之间的关联，从而避免了多重假设检验和校正的问题。

要了解WGCNA，我们需要了解WGCNA中的以下术语及其定义。

共表达网络：定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对相关性值进行冥次运算(冥次的值也就是软阈值 (power,pickSoftThreshold这个函数所做的就是确定合适的power))。无向网络的边属性计算方式为abs(cor(genex, geney)) ^ power；有向网络的边属性计算方式为(1+cor(genex, geney)/2) ^ power; signhybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0。这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。如果没有合适的power，一般是由于部分样品与其它样品因为某种原因差别太大导致的，可根据具体问题移除部分样品或查看后面的经验值。Module(模块)：高度內连的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可以对每个模块进行三个层次的分析：1. 功能富集分析查看其功能特征是否与研究目的相符；2. 模块与性状进行关联分析，找出与关注性状相关度最高的模块；3. 模块与样本进行关联分析，找到样品特异高表达的模块。基因富集相关文章去东方，最好用的在线GO富集分析工具；GO、GSEA富集分析一网打进；GSEA富集分析-界面操作。其它关联后面都会提及。Connectivity (连接度)：类似于网络中 “度”(degree)的概念。每个基因的连接度是与其相连的基因的边属性之和。Module eigengene E:给定模型的第一主成分，代表整个模型的基因表达谱。这个是个很巧妙的梳理，我们之前讲过PCA分析的降维作用，之前主要是拿来做可视化，现在用到这个地方，很好的用一个向量代替了一个矩阵，方便后期计算。(降维除了PCA，还可以看看tSNE)Intramodular connectivity:给定基因与给定模型内其他基因的关联度，判断基因所属关系。Module membership: 给定基因表达谱与给定模型的eigengene的相关性。Hub gene: 关键基因 (连接度最多或连接多个模块的基因)。Adjacency matrix(邻接矩阵)：基因和基因之间的加权相关性值构成的矩阵。TOM (Topological overlapmatrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。

基本分析过程

构建基因共表达网络：使用加权的表达相关性。识别基因集：基于加权相关性，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示。如果有表型信息，计算基因模块与表型的相关性，鉴定性状相关的模块。研究模型之间的关系，从系统层面查看不同模型的互作网络。从关键模型中选择感兴趣的驱动基因，或根据模型中已知基因的功能推测未知基因的功能。导出TOM矩阵，绘制相关性图。

WGCNA包实战

R-package WGCNA是一组计算各种加权关联分析的函数，可用于网络构建、基因筛选、基因簇识别、拓扑特征计算、数据模拟和可视化等。

输入数据和参数选择

WGCNA本质是基于相关系数的网络分析方法，适用于多样品数据模式，一般要求样本数多于15个。样本数多于20时效果更好，样本越多，结果越稳定。基因表达矩阵:常规表达矩阵即可，即基因在行，样品在列，进入分析前做一个转置。RPKM、FPKM或其它标准化方法影响不大，推荐使用Deseq2的varianceStabilizingTransformation或log2(x+1)对标准化后的数据做个转换。如果数据来自不同的批次，需要先移除批次效应(记得上次转录组培训课讲过如何操作)。如果数据存在系统偏移，需要做下quantile normalization。性状矩阵：用于关联分析的性状必须是数值型特征(如下面示例中的Height, Weight,Diameter)。如果是区域或分类变量，需要转换为0-1矩阵的形式(1表示属于此组或有此属性，0表示不属于此组或无此属性，如样品分组信息WT,KO, OE)。 ID WT KO OE Height Weight Diameter samp1 1 0 0 1 2 3 samp2 1 0 0 2 4 6 samp3 0 1 0 10 20 50 samp4 0 1 0 15 30 80 samp5 0 0 1 NA 9 8 samp6 0 0 1 4 8 7推荐使用Signed network和Robust correlation (bicor)。(这个根据自己的需要，看看上面写的每个网络怎么计算的，更知道怎么选择)无向网络在power小于15或有向网络power小于30内，没有一个power值可以使无标度网络图谱结构R^2达到0.8或平均连接度降到100以下，可能是由于部分样品与其他样品差别太大造成的。这可能由批次效应、样品异质性或实验条件对表达影响太大等造成,可以通过绘制样品聚类查看分组信息、关联批次信息、处理信息和有无异常样品(可以使用之前讲过的热图简化，增加行或列属性)。如果这确实是由有意义的生物变化引起的，也可以使用后面程序中的经验power值。

安装WGCNA

WGCNA依赖很多包装，生物导体上的包装需要自己安装，而起重机上的包装可以自动安装。一般WGCNA的安装可以通过在r中运行以下四个语句来完成。

建议在编译安装r时加入-with-blas-with-lapack，以提高矩阵运算速度。详见r和Rstudio安装。

source(“https://bioconductor.org/biocLite.R“) biocLite(c(“AnnotationDbi”, “impute”,”GO.db”, “preprocessCore”)) site=”https://mirrors.tuna.tsinghua.edu.cn/CRAN“ install.packages(c(“WGCNA”, “stringr”, “reshape2”), repos=site)

WGCNA实战

实战中采用政府提供的清洗矩阵，原始矩阵信息量太大，容易误导，后台回复WGCNA获取数据。

数据读入

图书馆

正在加载所需的包:dynamicTreeCut正在加载所需的包:fastcluster正在附加包:“fastcluster”以下对象被“pa Package:stats”屏蔽:h CLuster = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =* *重要说明:您的系统似乎支持多线程，*但在r中的WGCNA内未启用。*要在WGCNA内允许所有可用内核的多线程，请在r中使用* * allowgcnathreads()* *。如有必要，请使用disableWGCNAThreads()禁用线程。*或者，在您的系统上设置以下环境变量:* ALLOW _ WGCNA _ THREADS = * *例如* * ALLOW _ WGCNA _ THREADS = 48 * *要在linux bash shell中设置环境变量，请在运行r之前键入* * export ALLOW _ WGCNA _ THREES = 48 * *。其他操作系统或shell将*具有类似的命令来实现相同的目的。* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =附加包:“WGCNA”以下对象被“包:统计信息”屏蔽:cor

库(reshape2)库(stringr)

选项(字符串因子=假)

打开多线程

enableWGCNAThreads()

允许并行执行多达47个工作流程。

常规表达矩阵，log2转换后或 Deseq2的varianceStabilizingTransformation转换的数据如果有批次效应，需要事先移除，可使用removeBatchEffect 如果有系统偏移(可用boxplot查看基因表达分布是否一致)，需要quantile normalization

exprMat <。- "WGCNA/LiverFemaleClean.txt "

官方推荐 “signed” 或 “signed hybrid” 为与原文档一致，故未修改

type = "unsigned "

相关性计算官方推荐 biweight mid-correlation & bicor corType: pearson or bicor 为与原文档一致，故未修改

corType = "pearson "

corFnc = if else(Cortype = = " Pearson "，cor，bicor)

对二元变量，如样本性状信息计算相关性时，或基因表达严重依赖于疾病状态时，需设置下面参数

maxPOutliers = if else(Cortype = = " Pearson "，1，0.05)

关联样品性状的二元变量时，设置

robustY = if else(Cortype = = " Pearson "，T，F)

输入数据

dataExpr <。- read.table(exprMat，sep='t '，row.names=1，header=T，quote= "，comment= "，check.names=F)

dim(dataExpr)

[1] 3600 134

head(dataExpr)[，1:8]

F2 _ 2 F2 _ 3 F2 _ 14 F2 _ 15 F2 _ 19 F2 _ 20 MMT 000000044-0.01810 0.0642 6.44 e-05-0.05800 0.04830-0.15197410 MMT 0000046-0.07730-0.0297 1.12 e-01-0.05890 0 0.04 430

数据筛选

筛选中位数绝对偏差在75%之前的基因，至少MAD大于0.01，会减少计算量，丢失部分信息或者不筛选，使MAD大于0

m.mad & lt-应用(dataExpr，1，mad)DataExprvar & lt；-m . mad = " " gt。max(分位数(m.mad，probs=seq(0，1，0.25))[2]，0.01))，]

转换成一个矩阵，行中有样本，列中有基因

dataExpr <。- as.data.frame(t(dataExprVar))

缺失值的检测

gsg = goodSamplesGenes(dataExpr，verbose = 3)

标记缺失值过多的基因和样本…..步骤1

if(！gsg$allOK){

#或者，打印被删除的基因和样本名称:if (sum(！gsg$goodGenes)>;0)print flush(paste(" remaining genes:"，paste(names(dataExpr)[！gsg$goodGenes]，collapse = "，")；if (sum(！gsg$goodSamples)>。0) printFlush(粘贴("移除样本:"，粘贴(行名(dataExpr)[！gsg$goodSamples]，collapse = "，")；#从数据中删除违规基因和样本:数据表达式=数据表达式[gsg$goodSamples，gsg $ goodsgenes]

}

nGenes = ncol(DataExpr)NSAmples = nrow(DataExpr)

dim(dataExpr)

[1] 134 2697

head(dataExpr)[，1:8]

MMT 00000051 MMT 000000080 MMT 00000102 MMT 00000149 MMT 00000159 F2 _ 2-0.02260000-0.04870000 0 0.17600000-0.07680000-0.148000000 F2 _ 3 0.061700000

样本树= h样本(dist(dataExpr)，method = "average ")图(样本树，main = "样本聚类以检测异常值"，sub= "，xlab= " ")

软阈值的筛选原则是使构建的网络更符合无标度网络特性。

powers = c(c(1:10)，seq(from = 12，to=30，by = 2))sft = PickSoftthreshold(DataExpr，powerVector=powers，networkType=type，verbose=5)

pickSoftThreshold:将使用块大小2697。选择软阈值:计算给定功率的连通性…..正在研究SFT 2697人中的基因1到2697。R.sq斜率被截断。R.sq均值. k .中位数. k .最大值. k . 1 0.1370 0.825 0.412 587.000 5.95 e+02 922.0 2 2 0.0416-0.332 0.630 206.000 2.02 e+02 443.0 3 0.2280-0.747 0.920 91.500 8.43 e+01 247。 0.855 4.780 2.01 e+00 60.1 11 12 0.8020-1.410 0.866 3.160 9.61 e-01 53.2 12 14 0.8470-1.340 0.909 2.280 4.84 e-01 47。 7 13 16 0.8850-1.250 0.932 1.750 2.64 e-01 43.1 14 18 0.8830-1.210 0.922 1.400 1.46 e-01 39.1 15 20 0 0.9110-1.180 0.926 1.150 8.35 e-02 35.6 16 22 0.916 0-1.22

par(mfrow = c(1，2)) cex1 = 0.9

横轴是Soft threshold (power)，纵轴是无标度网络的评估参数，数值越高，网络越符合无标度特征 (non-scale)

标绘(SFTTFITINDICes[，3])SFTTFITINDICes[，1]，-符号(SFT $ FITINDICs[，3])SFT $ FITINDICs[，2]，标签=幂，cex=cex1，col="red ")

筛选标准。R-square=0.85

abline(h=0.85，col= "红色")

Soft threshold与平均连通性

绘图(SFtfitidices[，5]，xlab= "软阈值(功率)"，ylab= "平均连通性"，type="n "，main = paste("平均连通性"))文本(SFtfitidices[，5]，标签=功率，cex=cex1，col="red ")

功率= sft $功率估计功率

[1] 6经验功率(在没有合格功率时选择)

无向网络在power小于15或有向网络power小于30内，没有一个power值可以使无标度网络图谱结构R^2达到0.8，平均连接度较高如在100以上，可能是由于部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。

if(is . na(power)){ power = if else(NSA simples & lt；20，ifelse(type == "unsigned "，9，18)，if else(NSAmples & lt；30，ifelse(type == "unsigned "，8，16)，if else(NSAmples & lt；40，ifelse(type == "unsigned "，7，14)，ifelse(type == "unsigned "，6，12)))}

网络构建一步网络构建和模块检测

power: 上一步计算的软阈值 maxBlockSize: 计算机能处理的最大模块的基因数量 (默认5000)； 4G内存电脑可处理8000-10000个，16G内存电脑可以处理2万个，32G内存电脑可以处理3万个计算资源允许的情况下最好放在一个block里面。 corType: pearson or bicor numericLabels: 返回数字而不是颜色作为模块的名字，后面可以再转换为颜色 saveTOMs：最耗费时间的计算，存储起来，供后续使用 mergeCutHeight: 合并模块的阈值，越大模块越少

net = blockwiseModules(dataExpr，power = power，maxBlockSize = nGenes，TOMType = type，minModuleSize = 30，reassignThreshold = 0，mergeCutHeight = 0.25，numericLabels = TRUE，pamRespectsDendro = FALSE，saveTOMs=TRUE，corType = corType，maxPOutliers=maxPOutliers，loadTOMs=TRUE，saveTOMFileBase = paste0(exprMat，")。tom”)，verbose = 3)

从所有基因分块计算模块特征基因标记缺失值过多的基因和样本…..步骤1..在第一区工作。汤姆计算:邻接....将使用47个平行螺纹。慢速计算分数:0.000000..连通性....矩阵乘法....正常化....完成。..正在将块1的TOM保存到文件WGCNA/LiverFemleclean . txt . TOM-block . 1 . rdata…聚类..….检测模块..….计算模块特征基因..….检查模块中的kME....从模块1中删除3个基因，因为它们的KME太低。..从模块12中移除5个基因，因为它们的KME太低。..从模块14中删除1个基因，因为它们的KME太低。..合并过于接近的模块..合并关闭模块:合并距离小于0.25的模块计算新的MEs…

根据模块中基因数目的多少，降序排列，依次编号为 1-最大模块数。 0 (grey)表示未分入任何模块的基因。

表格(净价$colors)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 135 472 356 333 307 303 177 158 102 94 69 66 63 62分级聚类树显示每个模块的灰色基因，这些基因没有被分类到模块中。

Convert labels to colors for plotting

module labels = net $ colors module colors = labels 2 colors(module labels)

Plot the dendrogram and the module colors underneath 如果对结果不满意，还可以recutBlockwiseTrees，节省计算时间

绘图树和颜色(网络块基因[[1]]]，“模块颜色”，树标签=假，悬挂= 0.03，添加引导=真，引导悬挂= 0.05)

画出模块之间的关联图

module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示

MEs =净$MEs

不需要重新计算，只需更改以下名称即可。官方教程重新计算，不用一开始就费心了

MeS _ col = MeS col names(MeS _ col)= paste 0(" ME "，labels 2 colors(as . numeric(str _ replace _ all(col names(MeS)，" ME "，" "))))MEs_col = orderMEs(MEs_col)

根据基因间表达量进行聚类所得到的各模块间的相关性图 marDendro/marHeatmap 设置下、左、上、右的边距

绘图特征基因网络(MEs_col，“特征基因邻接热图”，marDendro = c(3，3，2，4)，marHeatmap = c(3，4，2，2)，绘图树图= T，xLabelsAngle = 90)

如果有表型数据，也可以和ME数据放在一起，一起作图

MEs_colpheno = orderMEs(cbind(MEs_col, traitData)) plotEigengeneNetworks(MEs_colpheno, “Eigengene adjacency heatmap”, marDendro = c(3,3,2,4), marHeatmap = c(3,4,2,2), plotDendrograms = T, xLabelsAngle = 90)

视觉基因网络

如果采用分步计算，或设置的blocksize>=总基因数，直接load计算好的TOM结果否则需要再计算一遍，比较耗费时间 TOM = TOMsimilarityFromExpr(dataExpr, power=power, corType=corType, networkType=type)

加载(net$TOMFiles[1]，verbose=T)

加载对象:TOM

汤姆& lt矩阵

dissTOM = 1-TOM

Transform dissTOM with a power to make moderately strong connections more visible in the heatmap

普罗托姆= dissTOM^7

Set diagonal to NA for a nicer plot

diag(绘图)= NA

Call the plot function 这一部分特别耗时，行列同时做层级聚类

基因图(网络热图，所有基因)

细胞景观出口网络

细胞角绘制的网络图见我们更新的视频教程或https://bioinfo.ke.qq.com/。

探针=列名(DataExpr)dim names(TOM)& lt；-列表(探针、探针)

Export the network into edge and node list files Cytoscape can read threshold 默认为0.5, 可以根据自己的需要调整，也可以都导出后在 cytoscape中再调整

cyt = exportNetworkToCytoscape(TOM，edgeFile = paste(exprMat，" . edges.txt "，sep= ")，nodeFile = paste(exprMat，" . nodes.txt "，sep= ")，weighted = TRUE，threshold = 0，nodeNames =探测器，nodeAttr = moduleColors)