WGCNA 分析

---

2021.9.28 参考了生信技能树的教程
第一次学习WGCNA 分析(这个属于单细胞高级分析)
还有细节未探究,WGCNA+GO+KEGG功能注释

参考：
https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/
https://cloud.tencent.com/developer/article/1516749
https://zhuanlan.zhihu.com/p/388190968
https://www.jianshu.com/p/817d492d4904
https://mp.weixin.qq.com/s?__biz=MzAxMDkxODM1Ng==&mid=2247491248&idx=1&sn=afb11345a1e89cd57ace7dcf827dddbf&scene=21#wechat_redirect

#应用场景：
1.找到模块的核心基因

2.利用关系预测基因功能

#一般流程：

得到模块之后的分析有：
1.模块的功能富集

2.模块与性状之间的相关性

3.模块与样本间的相关系数

#主要函数：

### 计算邻接矩阵
adjacency = adjacency(dataExpr, power = power)

### 把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得距离矩阵。
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM

### 层级聚类计算基因之间的距离树 
geneTree = hclust(as.dist(dissTOM), method = "average")

### 模块合并
# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)

### 通过计算模块的代表性模式和模块之间的定量相似性评估，合并表达图谱相似
的模块
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs)
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")
MEDissThres = 0.25

# Call an automatic merging function
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged

#基本概念

• Module
  模块(module)：表达模式相似的基因分为一类，这样的一类基因成为模块；
• Eigengene
  Eigengene（eigen + gene）：基因和样本构成的矩阵，https://en.wiktionary.org/wiki/eigengene
• Adjacency Matrix
  邻近矩阵：是图的一种存储形式，用一个一维数组存放图中所有顶点数据；用一个二维数组存放顶点间关系（边或弧）的数据，这个二维数组称为邻接矩阵；在WGCNA分析里面指的是基因与基因之间的相关性系数矩阵。 如果用了阈值来判断基因相关与否，那么这个邻近矩阵就是0/1矩阵，只记录基因相关与否。但是WGCNA没有用阈值来卡基因的相关性，而是记录了所有基因之间的相关性。
• Topological Overlap Matrix (TOM)
  利用上面的邻近矩阵，又计算了一个新的邻近矩阵。一般来说，TOM就是WGCNA分析的最终结果，后续的只是对TOM的下游注释。

#原理解析：

#Soft Threshold ？Hard Threshold？
网络中基因与基因之间是否连边取决于这两个基因是否发生显著共表达：nodes represent genes and nodes are connected if the corresponding genes are significantly co-expressed across appropriately chosen tissue samples

皮尔森相关系数(Pearson) correlation coefficient：It is standard to use the (Pearson) correlation coefficient as a co-expression measure, e.g., the absolute value of Pearson correlation is often used in a gene expression cluster analysis.

那么现在问题来了：我得到了两个基因的相关系数值之后，如何决定这两个基因在构建网络时是否连边呢？

1，定义一个阈值，比如有统计学意义的r>0.9，那么就连边；否则，就不连边。这个办法就是‘hard’ threshold，也即不加权（unweighted）。
2，真实网络与随机网络的特征差异
真实的生物学网络比如：

human disease network，
gene co-expression networks，
protein-protein interaction networks，
cell-cell interaction networks，
the world wide web and social interaction networks。

• 真实网络服从幂率分布（power law） p(k)~ k−γ，即无标度网络(Scale-free networks)
  

• 随机网络一般服从泊松(Poisson)分布
  

这两个网络分布的特点是：真实网络中绝大部分点的度都很低，只有少数的点（即常说的hub节点）度很高，而随机网络中绝大部分的点的度都处于上图中的峰值处，即度相似。

无标度网络节点的度分布特征使得它有一个很大的特点，那就是稳健性：随机去除网络中的一个节点，网络还能依然保持（绝大多数节点的度很低）。但是它也有脆弱的一面：那就是去除网络中的hub节点，网络就散了。但是hub节点只是网络中极少数的几个节点，被攻击的概率非常小。

基于以上，即‘hard’ threshold的缺点和应用于真实网络造成的信息丢失（参考文献中有说明）等原因，官网教程中的作者提出了一种方法即用‘soft’ thresholding 来权衡网络中的连边。

就是将相关系数的值[0,1]通过一个换算映射到[0,1]，也即加权（weighted）的思想：

那么这里的参数β就是软阈值。

---

#如何选取Soft Threshold？

左图的纵坐标scale-free fit index，即signed R2，代表对应的网络中log(k)与log(p(k))相关系数的平方乘以一个方向向量，由slope决定（The sign of the scale-free model fitting index R2 is determined by minus the sign of the slope），拟合的线性方程为下图，来源于WNCGCNA包中的源代码：
相关系数的平方越高，说明该网络越逼近无标度网络的分布。相关参考文献中有大量数据证明当signed R2 大于0.85时，网络就已经符合无标度网络的分布。
因此，WGCNA包中计算SoftThreshold的函数pickSoftThreshold中RsquaredCut 默认值为 0.85，最佳的powers值保存在sft$powerEstimate。上图中，作者重新设定了一个阈值cex1=0.9，因此你会看到图片作图中有一条红色的横线，表示第一个signed R2达到这条红线时的最佳powers值，此图中是6。

右图的纵轴代表对应的网络中所有基因连接数（即节点的度）的均值。

这张图有三个部分：

1，聚类树：那么你对聚类方法又了解多少？
2，Dynamic Tree Cut：根据某一height，聚类树分成了多少类（也即模块）
3，Merged dynamic：合并后的聚类模块，那么又是根据什么指标要将Dynamic Tree Cut中的类别进行合并得到新的聚类模块呢？
根据最佳软阈值和表达谱数据首先计算得到一个邻接矩阵adjacency，然后将邻接矩阵转换为TOM矩阵再计算一个距离矩阵dissTOM，这个距离矩阵就是后面用来画聚类树的图。

• 层次聚类

We now use hierarchical clustering to produce a hierarchical clustering tree (dendrogram) of genes. Note that we use the function hclust that provides a much faster hierarchical clustering routine than the standard hclust function.

geneTree = hclust(as.dist(dissTOM), method = "average");

采用的聚类方法为平均距离（method = “average”），

图片中，每一个树枝代表一个gene，纵坐标的Heigth为聚类距离。

设定一个聚类距离阈值，就可以将gene聚成module，这里又有很多方法设定这个阈值。而官网教程成给定的是dynamicTreeCut包中的函数cutreeDynamic，并且设定了最小模块中包含的基因个数为30个基因。

# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30;

# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                deepSplit = 2, pamRespectsDendro = FALSE,
                minClusterSize = minModuleSize);

table(dynamicMods)

为了定量每一个模块的共表达相似性，这里计算了每一个模块的eigengenes，即每个模块的特征值。

# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs