经过软阈值确定，TOM矩阵计算，层次聚类和动态切割后，

我们获得了一些基因模块，每个模块里面都可以看成是有共同目的的团伙。

为了方便后面纷至沓来的相关性分析，尤其是跟性状的联合分析，现在需要选取出每个模块中的代表人物。

这个过程就是寻找每个模块的eigengene(读音: 爱根金)

这样，模块信息就变成了，行是eigengene，列是样本的矩阵了，十分简化。

什么是eigengene

按照作者的说法，就是每个模块基因和样本矩阵进行PCA分析后的第一主成分。

WGCNA中提供了简单的方法来计算。

就是moduleEigengenes这个函数，名字也比较好记，就是模块的Eigengene。

输入表达量数据和模块的信息就可以。

	

	

	


		

		


			moduleEigengenes(datExpr, moduleColors) 

		




	

返回的是个列表，可以把我们需要的提取出来。

	

	

	


		

		


			MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes 

		




	

行是样本，列是模块的eigengene

有了这个就可以跟性状信息求相关性，得到WGCNA最重要的热图了。

现在的问题是，这个值是怎么算出来的。

我们以第一列黑色模块来举例。

我们重复一下，是主成分分析的第1个主成分。

那就先把黑色模块提取出来，提取模块也是个重要操作。

因为后面确定了模块，肯定就需要提取模块

	

	

	


		

		


			table(moduleColors) 

		




	

黑色模块含有211个基因，提取表达矩阵。

	

	

	


		

		


			datablack <- datExpr[,moduleColors == "black"] 

		




	

行是样本，列是基因

有一些缺失值，先用knn技术把缺失值填充一下。

	

	

	


		

		


			### 缺失值填充 

		




		


			### 要求就是行是基因，列是样本 

		




		


			datModule = t(as.matrix(datablack)) 

		




		


			datModule = impute::impute.knn(datModule, k = min(10, nrow(datModule)-1)) 

		




		


			## 提取数据，行是样本，列是基因 

		




		


			datModule_knn = t(datModule$data) 

		




	

	

	

	


		

		


			pc <- prcomp(datModule_knn) 

		




	

提取第一个主成分

	

	

	


		

		


			pc$x[, 1] 

		




	

	

	

	


		

		


			[1]  0.0138038628  0.0668404876  0.0667271646 -0.0643323004  0.0637495925 -0.0008963038 

		




	

发现不怎么一样。

但是你如果计算一下相关性，会发现新的结果。

	

	

	


		

		


			cor(MEs0[,1],pc$x[, 1]) 

		




	

相关性极其高

	

	

	


		

		


			           [,1] 

		




		


			[1,] -0.9994166 

		




	

说明两种有不可告人的关系。

奇异值分解

进过阅读文献，查看源码得知，作者在计算eigengene的时候用的是奇异值分解法。

这里是让我十分沮丧的，因为奇异值分解用到了高数矩阵的知识，我看了很多资料，大脑就是缺失一环。

所以，当前我并不能深入浅出的讲清楚这个事情，搞懂是迟早的，只是当前优先度不够。

但是，我可以做出来。

首先奇异值分解，数据是需要缩放的，就用scale函数

scale函数的理解，可以看看这一篇。

z-score的标准化究竟怎么弄？

	

	

	


		

		


			datModule_scale= scale(datModule_knn) 

		




	

那么复杂的技能奇异值分解，就是一个单词的事情

	

	

	


		

		


			svd = svd(datModule_scale) 

		




	

而结果里面的左边奇异值的第一列就是我们要的结果

	

	

	


		

		


			svd$u[,1] 

		




	

数据证明是数值是对的，只是符号不一样。

而主成分得到的结果和分解奇异值得到的结果高度相关，只是数值不同

	

	

	


		

		


			cor(pc$x[, 1],svd$u[,1]) 

		




	

	

	

	


		

		


			          [,1] 

		




		


			[1,] 0.9994166 

		




	

但是现在的问题是，我们算出来的结果跟官方的有出入，正负相反。

他们处理的方法是，计算每个样本中，基因的平均值，用这个平均值跟主成分求相关性

如果相关性是负的，那么就把主成分的结果乘以负一。

	

	

	


		

		


			### 计算平均值 

		




		


			averExpr = rowMeans(datModule_scale, na.rm = TRUE)    

		




		


			### 求相关性 

		




		


			corAve = cor(averExpr,svd$u[,1], use = "p") 

		




		


			corAve  

		




	

确实是负数

	

	

	


		

		


			          [,1] 

		




		


			[1,] -0.998023 

		




	

明明可以用主成分算的事情，为什么要用奇异值分解？

又是算法层面的解释：

PCA当样本非常多的时候，计算协方差矩阵，还要进行特征值分解，计算量很大，而有些奇异值分解方法可以跳过这一步直接得到结果。

eigengene的意义在哪？

第一，性状关联。

模块里面的基因，经过压缩，变成了一个eigengene，这样跟性状信息可以求相关性。

有了相关性，我们就可以得到性状和模块的相关性热图，这也是WGCNA里面最重要的一张图。

那模块中的基因怎么提取呢？我们已经展示了。

提取出的基因可以画热图，GO分析，KEGG分析。

第二，hubgene

eigengene的存在，让找每个模块中的是hub gene变得简单。

hubgene，通俗的解释就是

买个模块中跟eigengene 相关性最强的基因。

那些cytoscape画的网络图都是次要的，因为找hub gene不需要画图，直接计算就行了。

第三，迭代WGCNA

这里面来了一个新的概念就是迭代WGCNA，我们后面会介绍。

eigengene的出现，我们筛选模块里面高价值的基因就有了标准，我们可以把模块中的每个基因都跟eigengene计算相关性，把小于0.8的全部剔除(人为定义)

每个模块都这么做，现在就到了一个新的行是基因，列是样本的矩阵。

重新来做一次WGCNA，再来筛选一次，直到无法筛选基因为止。

此时，我们就得到了都比较纯净模块。而有些文章，模块基本上都是模糊的，胆子大，心也大。