lncRNA实战项目-第六步-WGCNA相关性分析

响应生信技能树的号召：lncRNA数据分析传送门, 一起来一个lncRNA数据分析实战！

因为样本数量比较可观，所以可以进行WGCNA分析。这里是并不需要选取所有的基因来做WGCNA分析，挑选的标准可以是top变异程度大的基因集合，或者显著差异表达的基因集合等等。

这里可以参考：https://github.com/jmzeng1314/my_WGCN

WGCNA将lncRNA分成18个模块（3635个lncRNA），空间模块中lncRNA表达呈现明显的组织区域特异性，如：CB (M1, 794个lncRNAs),DG/CA1 (M2, 443个lncRNAs), CA1 (M4, 369个lncRNAs),neocortex (M7, 123个lncRNAs)和OC (M10,57个lncRNAs)。时间模块中lncRNA表达与年龄有关，而与组织区域不明显；性别模块中lncRNA表达与性别和年龄都相关。每个模块就必须做pathway/go等数据库的注释分析咯！

资料收集：

google搜索或在生信技能树和生信菜鸟团搜索WGCNA ，能找到很多教程，下面列出几个中文教程和英文教程，强烈推荐中文教程1和英文教程3。

一文学会WGCNA分析

学习WGCNA总结

Tutorials for the WGCNA package

WGCNA Background and glossary

Data input and cleaning

Network construction and module detection

Relating modules to external information and identifying important genes

Interfacing network analysis with other data such as functional annotation and gene ontology

Network visualization using WGCNA functions

Exporting a gene network to external visualization software

背景知识

基本概念:

WGCNA(weighted correlation network analysis)加权基因共表达网络分析, 用于提取与性状或临床特征相关的基因模块，解析基础代谢途径，转录调控途径、翻译水平调控等生物学过程。WGCNA适合于复杂的数据模式，推荐5组以上的数据，如：

不同器官、组织类型发育调控；

同一组织不同时期发育调控；

非生物胁迫不同时间点应答；

病原物侵染后不同时间点应答。

基本步骤：

WGCNA分为表达量聚类分析和表型关联两部分，具体步骤包括基因之间相关系数的计算，共表达网络的构建，筛选特定模块，模块与性状关联，核心基因的筛选。

Overview of a typical WGCNA analysis.

术语：

Co-expression weighted network:是一个无向有权重（undirected, weighted）的网络。“无权重（unweighted network）”,基因与基因之间的相关度只能是0或者1,0表示两个基因没有联系，而1表示有。“有权重（weighted network）”基因之间不仅仅是相关与否，还记录着它们的相关性数值，数值就是基因之间的联系的权重(相关性)。

Module:（模块）指表达模式相似的基因聚为一类，这样的一类基因称为模块。

Connectivity（连接度）:连接度是指某个基因与该网络中其他基因的相关程度，常为“相关度”之和。

Eigengene（eigen- +‎ gene）：基因和样本构成的矩阵。

Hub gene: 模块中的高度相关的核心基因。

Gene signicance ，GS:基因显著性参数，为非负数字信息，比如基于样本的临床数据（clinical trails）和基于每个基因的-log(p-value)等。

分析流程

WGCNA输入文件需要一个表达矩阵，最好是RPKM或其他归一化好的表达量，还需要一个矩阵关于临床信息或者其它关于样本属性的信息。

STEP1: 输入数据的准备

表达矩阵可以从作者GitHub上下载（https://github.com/DChenABLife/RhesusLncRNA），我这里只尝试对lncRNA的表达矩阵(https://github.com/DChenABLife/RhesusLncRNA/blob/master/All_sample_LncRNA_exp_RPKM.xlsx)做后续分析。下载的表达矩阵文件是Excel格式的，需要转为csv格式方便后续用R处理，可以直接打开这个excel文件，然后另存为csv格式即可。

读入原始表达数据

原始数据包含64个样本，9904个lncRNA表达量，这时的矩阵行为基因，列为样本信息，其中第一列是lncRNA ID，(这里的lncRNA ID是cufflinks 组装时给的自由的ID, 是需要和已有的ID对应, 对于新的转录本再通过nr/nt等数据库注释）， 第66列是作者给出的注释（我查了几个注释有的也查不出来是什么意思）。

rawdata

重命名数据列表，行名和列名

fpkm

准备表型信息

这里有64个样本，包含猕猴脑不同空间区域，不同发育时期，以及性别，因为每个样本都交叉包含着三种不同的信息，如果选择全部表型信息，我试了试，后续的模块和性状完全看不清关系，所以我这里仅选择脑不同区域的表型信息,包括CB、DG、PFC、PCC、CA1、OC、PC、TC。

sample.Info

下载并载入WGCNA包

行列转置

WGCNA针对的是基因进行聚类，而一般我们的聚类是针对样本用hclust即可，也就是说要转置为行表示样本, 列表示基因。

datExpr

确定临床表型与样本名字

datExpr和datTraits准备好后，接下来就是构建基因网络，鉴定模块。网络构建有三种方法：1）一步法构建网络；2）多步法构建网络；3）block-wise构建网络（主要针对大数据集）。下面的介绍的步骤是一步法构建网络。

STEP2：确定最佳soft-thresholding power

选择合适的“软阀值（soft thresholding power）”beta, 用到的函数是pickSoftThreshold，pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)，powerVector可以是一系列数值，从而选择最优值。这个函数返回一个列表，第一项是powerEstimate是估计的最优power；第二项是fitIndices是详细的矩阵数据，其中第五列是mean.k表示平均“连接度（connectivity）”。

Constructing a weighted gene network entails the choice of the soft thresholding power to which co-expression similarity is raised to calculate adjacency.

最佳beta值是3。

Soft Threshold

STEP3: 一步法构建共表达矩阵

一步法构建网络，使用函数blockwiseModules(), 这个函数包含很多参数，其中power=sft$powerEstimate=3 即上一步得到的最佳软阈值；maxBlockSize 默认为5000， 表示在这个数值内的基因将整体被计算，如果调大需要更多的内存；numerricLabels 默认为返回颜色，设置为TRUE则返回数字；mergeCutHeight是合并模块阈值的一个参数。

STEP4:模块鉴定及可视化

模块鉴定

上一步的返回结果是一个列表，可以用table()函数查看，0表示没有任何module接受。table(net$colors) 可以看总共有多少模块，每个模块的大小，这里共有9个模块，从1-9每个模块的大小是递减的，从2254-115，0表示这些基因不在所有模块内。

可视化

dendrograms表示在一个block中所有基因的进化树图, 使用函数plotDendroAndColors()查看系统发生树；blockGenes是一个block中所有的基因。

cluster Dendrogram

Sample dendrogram and trait heatmap

STEP5:模块和性状的关系

MEs是一个关于modules的特征量矩阵，行数等于筛选的modules数，列数等于样本数；

图中第二列第五行，即CB/turquoise相关性最强（0.97），pvalue=1e-41,后续分析可以挑选这个模块。

每一列对应的样本特征可以从design这里查看。

STEP6:感兴趣性状的模块的具体基因分析

选定CB/turquoise这个模块后，接下来看看模块与样本特征和基因的相关性，1. 首先计算全模块与基因的相关性矩阵，2. 再计算性状与基因的相关性矩阵，3. 使用verboseScatterplot()绘制相关性图：

首先计算模块与基因的相关性矩阵

再计算性状与基因的相关性矩阵

最后把两个相关性矩阵联合起来,指定感兴趣模块进行分析

上图可以看出基因跟其对应的性状高度相关，可以导出做个GO/KEGG注释，看看这些基因的具体功能。

STEP7:网络的可视化

Topological Overlap Matrix (TOM)图是指基因和筛选模型的可视化表示。WGCNA认为基因之间的简单的相关性不足以计算共表达，所以它利用邻近矩阵，又计算了一个新的邻近矩阵。一般来说，TOM就是WGCNA分析的最终结果，后续的只是对TOM的下游注释。

生成TOM图的步骤：1. 生成全基因不相似TOM矩阵，比如dissTOM=1-TOMsimilarityFromExpr(datExpr, power = 6)，可以把得到的不相关矩阵加幂，这样绘制的TOM图色彩差异会比较明显。2. 使用TOMplot()绘制TOM图。

图上和图左是全基因系统发育树，不同颜色亮带表示不同的module，每一个亮点表示每个基因与其他基因的相关性强弱（越亮表示相关性越强）。

STEP8:提取指定模块的基因名

提取基因信息,可以做GO/KEGG等分析，进而解释这些module的生物学意义。这里的lncRNA ID转换着有点麻烦，这一步先略过，之后再看看。

STEP9: 模块的导出

主要模块里面的基因直接的相互作用关系信息可以导出到cytoscape,VisANT等网络可视化软件。

首先是导出到VisANT

也可以输出模块里的部分基因：

然后是导出到cytoscape

编辑：jimmy

后续分析，请大家持续关注

~