linux系统如何进行聚类树,WGCNA分析使用教程

加权基因共表达网络分析 (WGCNA, Weighted correlation network analysis)是用来描述不同样品之间基因关联模式的系统生物学方法，可以用来鉴定高度协同变化的基因集, 并根据基因集的内连性和基因集与表型之间的关联鉴定候补生物标记基因或治疗靶点。

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

理解WGCNA，需要先理解下面几个术语和它们在WGCNA中的定义。共表达网络：定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对相关性值进行冥次运算 (冥次的值也就是软阈值 (power, pickSoftThreshold这个函数所做的就是确定合适的power))。无向网络的边属性计算方式为 abs(cor(genex, geney)) ^ power；有向网络的边属性计算方式为 (1+cor(genex, geney)/2) ^ power; sign hybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0。这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。如果没有合适的power，一般是由于部分样品与其它样品因为某种原因差别太大导致的，可根据具体问题移除部分样品或查看后面的经验值。

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

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

Module eigengene E: 给定模型的第一主成分，代表整个模型的基因表达谱。这个是个很巧妙的梳理，我们之前讲过PCA分析的降维作用，之前主要是拿来做可视化，现在用到这个地方，很好的用一个向量代替了一个矩阵，方便后期计算。(降维除了PCA，还可以看看tSNE)

Intramodular connectivity: 给定基因与给定模型内其他基因的关联度，判断基因所属关系。

Module membership: 给定基因表达谱与给定模型的eigengene的相关性。

Hub gene: 关键基因 (连接度最多或连接多个模块的基因)。

Adjacency matrix (邻接矩阵)：基因和基因之间的加权相关性值构成的矩阵。

TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。

基本分析流程

构建基因共表达网络：使用加权的表达相关性。

识别基因集：基于加权相关性，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示。

如果有表型信息，计算基因模块与表型的相关性，鉴定性状相关的模块。

研究模型之间的关系，从系统层面查看不同模型的互作网络。

从关键模型中选择感兴趣的驱动基因，或根据模型中已知基因的功能推测未知基因的功能。

导出TOM矩阵，绘制相关性图。

WGCNA包实战

R包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依赖的包比较多，bioconductor上的包需要自己安装，cran上依赖的包可以自动安装。通常在R中运行下面4条语句就可以完成WGCNA的安装。

建议在编译安装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 获取数据。

数据读入library(WGCNA)

## Loading required package: dynamicTreeCut

## Loading required package: fastcluster

##

## Attaching package: 'fastcluster'

## The following object is masked from 'package:stats':

##

## hclust

## ==========================================================================

## *

## * Package WGCNA 1.63 loaded.

## *

## * Important note: It appears that your system supports multi-threading,

## * but it is not enabled within WGCNA in R.

## * To allow multi-threading within WGCNA with all available cores, use

## *

## * allowWGCNAThreads()

## *

## * within R. Use disableWGCNAThreads() to disable threading if necessary.

## * Alternatively, set the following environment variable on your system:

##