WGCNA有了相关性矩阵为什么还要计算拓扑矩阵?

我和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的分析。

生物信息学

经验分享:转录组测序(比对篇)

2020-6-22 21:16:29

生物信息学

利用TCGA数据库发表的文章都是哪些杂志?

2020-6-23 20:13:05

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索