我和WGCNA的纠葛
WGCNA的我忍了他四年了。
第一次听到是一位朋友需要学习WGCNA，公司给他报的价格是3万。
然后他就硬着头皮自学，然后发了文章毕业了。
这大概是四年前的事情。
从那个时候，我就知道了WGCNA的名字，但是一直想学但就是没有行动。
因为那时候我的R语言没有入门，看到别人limma包都写了本书来说明，就觉得所有R包的学习都需要很长时间。
直到有一天Jimmy跟我说，他学习一个新的R包就跟喝水一样，彻底把我震撼了。
我说我想看一下你学习R包的过程。
他就开了一个网络讲座，现场直播WGCNA的学习过程，在视频里面他说是要给一个朋友演示过程，那个朋友就是果子。
但是，不争气的是，那个时刻我见识到他从头到尾学会了WGCNA，后来整个公众号的帖子都是以他的流程作为模板来分析，但是，我还是觉得太难。
现在想来，是因为基础太差，数据的增删改查都不能独自完成。
怎么样做才能不辜负Jimmy的这番好意呢？
假装会了。
我确实也是这么做的，总不能说自己不会吧。
一开始真的很心慌，因为怕抽查，到了后面就坦然了，因为我觉得随时都可以学会。
到了这个月，组上的样本堆在我面前，我不能再这么下去了，就开始学习。
学习也是直接看的官方教程，然后就开始拆解他的代码。
拆解的第一部分弄清楚了，什么叫无尺度网络分布
挑选合适的软阈值是为了让基因间的相关性符合无尺度网络分布。
接下来就是聚类，找出模块。
拓扑重叠矩阵如何计算出来
但是WGCNA在找模块之前多做了一步，就是把相关性矩阵（Adjacency Matrix）变成了拓扑重叠矩阵(Topology Overlap Matrix)。
这个操作是依赖TOMsimilarity函数来完成的。
## 加载R包
library(WGCNA)
## 确定软阈值
softPower = 6
## 得到邻接矩阵
adjacency = adjacency(datExpr, power = softPower)
## 转换为拓扑矩阵
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
上面的代码中，计算邻接矩阵的这一步
adjacency = adjacency(datExpr, power = softPower)
实际上等同于之前讲的相关性分析
adjacency=abs(cor(datExpr,use = 'p'))^power
是这个样子的
就这个邻接矩阵，已经可以通过聚类来找出模块了
hierADJ <- hclust(as.dist(1-adjacency),method="average")
plot(hierADJ,labels=F, hang = 0.04)
但是作者觉得，仅仅通过表达的相关性矩阵，不足于反应体内真实的情况。
假如基因1和基因2的表达相关性是0.8，基因2和基因3的表达相关性也是0.8
你是不能说这两个关系是相等的，因为基因1和基因2之间还有很多共同的小伙伴跟着两个基因建立联系。
也就是说，我们评价两个基因相关性的时候，不能仅仅看这个两个基因的表达相关性，还要考虑跟其他基因之间的互作。
那么其他基因之间的互作怎么表示呢？
用和其他基因相关性的乘积之和来表示。
举个例子：
假如现在矩阵A是这个样子的，选取A21，也就是第二行第一列的数是2
经过这个式子之后是多少呢？
是
A21*A11+A22*A21+A23*A31+A24*A41
也就是
2*1+6*2+10*3+14*4=100
在R语言里面有个符号可以方便的做
A %*% A
在计算好的值后面再加上一个本身，就是100+2最终就是102
A %*% A+A
那方程式下面是干什么呢？
我解释不了，作者说的是，从数学上让整个式子的值在0到1之间，就这个作用。
我们现在从头来算一次。
adjacency = adjacency(datExpr, power = softPower)
### 现在分布操作
adjmat1=adjacency
### 对角线为0
diag(adjmat1)=0
### NA的变为0
adjmat1[is.na(adjmat1)]=0
adjmat1现在是这样的
现在就是表达式的上面部分,
numTOM=as.dist(adjmat1 %*% adjmat1 +adjmat1)
然后计算分母
### 获取每个基因的连通性
kk=apply(adjmat1,2,sum)-1
## 新的矩阵，行列跟基因一样
Dhelp1=matrix(kk,ncol=length(kk),nrow=length(kk))
denomTOM= pmin(as.dist(Dhelp1),as.dist(t(Dhelp1)))+as.dist(1-adjmat1)
最终算出TOM matrix
TOMmatrix=as.matrix(numTOM/denomTOM)
算出来的是相关性，用1减去这个值，就变成不相信，便于画图展示
## 变成dissimilarity
dissTOM = 1-TOMmatrix
## 对角线变成1
diag(dissTOM) = 1
以上步骤实际上被写成了一个函数，可以从表达量直接计算
## 加载R包
library(WGCNA)
## 确定软阈值
softPower = 6
## 得到邻接矩阵
adjacency = adjacency(datExpr, power = softPower)
## 转换为拓扑矩阵
TOM = TOMsimilarity(adjacency)
dissTOM = 1-TOM
如果仔细观察，发现值不怎么一样，这是怎么回事呢？
是这个内置函数处理数据的问题，他在计算连通性的时候，没有减去自身的1
kk=apply(adjmat1,2,sum)
而我们是这样计算的,更加准确
kk=apply(adjmat1,2,sum)-1如果我们改成跟他一样，那么算出来肯定也是一样的。
说明我们的想法完全正确。
现在相关性矩阵（Adjacency Matrix）变成了拓扑重叠矩阵(Topology Overlap Matrix)。
只要聚类就可以发现有很多branch，分支，这些就是模块了
geneTree = hclust(as.dist(dissTOM), method = "average")
plot(geneTree,labels = FALSE, hang = 0.04)
下次我们看看他如何确定模块以及合并模块的，然后再讲讲每个模块的eigengene是如何计算出来的。
在结束的那一刻，我们就搞懂了WGCNA所有的细节。
也意味着，不加载WGCNA的包，我们也可以自己进行WGCNA的分析。
