WGCNA的分析中为什么要挑选软阈值?

我们今天就搞清楚这两个问题:
1.什么是无尺度分布
2.手工计算软阈值

无尺度分布

假如我问大家,如果你想认识世界上的任何一个人,需要通过几个人来联系?

6个。

这个就是6度分隔理论,描述的是随机网络里面的情形,随机网络里面,大多数节点拥有相同的链接数。
如果用计算机模拟的。在节点完全随机的情况下,任何两个节点的平均距离却远远大于6个。
但是我会感觉6度理论是靠谱的,举两个例子:
因为你想联系上特朗普的话,应该没那么困难。你先认识地方电视台,然后通过他们联系马云,而马云跟特朗普认识。

你如果找马化腾呢?也没有那么困难,你先认识地方电视台,然后找到马云,而马云跟马化腾很熟。

你看,只要我们在这个网络中,加入一些巨人,事情就会简化很多,而这些巨人拥有的特点就是,巨人很少,单个巨人跟普通人的链接特别多。
这些巨人起到了枢纽的作用,枢纽英文叫作hub。

真实世界的网络不是随机的,而是有hub的网络。就像右边的图一样。普通节点占大多数,hub节点是少数。

随机网络是可以用一个值来衡量大多数节点之间的距离,或者是6或者是60,也就是你可以用一个尺度去衡量他,可以称为尺度网络。
而有巨人的网络,没办法来衡量两个节点之间的距离,叫作无尺度分布,该分布又叫做幂律分布,我们说的二八定律,长尾定律都是幂律分布的口头化呈现。

人和人之间的交往是这个样子,那么蛋白和蛋白之间的互作是什么情况呢?
目前的观察是,也符合幂律分布,也就是无尺度分布。一个细胞内,不是每个基因都要表达,即使表达,也不定起作用。
决定这个基因分化,功能改变,都是一些重要基因,我们把他成为hub gene。

好了,WGCNA分析就是利用了这一点规律,强行让基因间的联系符合无尺度分布。

手工计算软阈值

我们用真实的数据来演示以下这个过程

load(file = "FemaleLiver-01-dataInput.RData")

我们的数据名称是dataExpr,在第一次课的文件夹里面。
他的行是样本,134个,列是基因,3600个。

WGCNA首先通过相关性分析(cor函数),计算了任意两个基因的相关性。
如果是平常,我们就会设定一个阈值,大于0.8(假设是这个)是有联系,小于的叫没有联系。
但是正常情况并不是这样,因为0.79的基因会出来闹。

所以,作者就给这些相关性加了一个幂。

图中的绿线原来的相关性,进行幂次计算之后(8次)就变成了红色。
可以看到整体的相关性都变小了,但是里面本来就小的变得更加小。

为什么要这么做,因为这样做了之后,基因间的连通性,就开始符合幂律了。
我们以幂为10来计算一下。

power <- 10
ADJ=abs(cor(datExpr,use = 'p'))^power

现在得到了一个3600行乘以3600行的相关性矩阵

如果对每一列的基因进行求和,得到的就是这个基因跟其他基因相关性之和。

k=apply(ADJ,2,sum) -1

减去1是排除了自身。

MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080 MMT00000102 
 0.02094526  7.13906159  4.17142638  0.22300729  6.25771609  1.20155840

这个K指的就是3600基因组成的网络中,每个节点的连通性。做一个频次分布图就会有惊喜

hist(k)

横坐标是连通性依次升高,纵坐标表示该范围连通性的频次,是不是从分布上看符合了幂律呢?
当然这是靠眼睛在看,
真正的幂律是这样的,把连通性分隔,分隔内连通性的平均值取log10,跟频率的概率取log10,两者之间有线性关系。

用图展示一下,
先把k从小到大排序,切割成10份

cut1=cut(k,10)

计算每个区间的平均值

binned.k=tapply(k,cut1,mean)

结果如下

(-0.0655,6.58]    (6.58,13.2]    (13.2,19.7]    (19.7,26.3]    (26.3,32.9]    (32.9,39.5]      (39.5,46]      (46,52.6]    (52.6,59.2]    (59.2,65.8] 
      2.151281       9.145864      15.670490      22.564169      30.466627      36.649594      42.028070      49.043907      55.893530      61.864996

然后计算每个区间的频率

freq1=tapply(k,cut1,length)/length(k)

结果是这样的

(-0.0655,6.58]    (6.58,13.2]    (13.2,19.7]    (19.7,26.3]    (26.3,32.9]    (32.9,39.5]      (39.5,46]      (46,52.6]    (52.6,59.2]    (59.2,65.8] 
   0.792500000    0.130277778    0.036666667    0.011111111    0.001944444    0.003055556    0.003611111    0.008333333    0.005555556    0.006944444

此时这个均值的对数和频率的对数就是线性的

plot(log10(binned.k),log10(freq1+.000000001),xlab="log10(k)",ylab="log10(p(k))")

如果通用线性函数加线和注释就会明显一点

xx= as.vector(log10(binned.k))
lm1=lm(as.numeric(log10(freq1+.000000001))~ xx )
lines(xx,predict(lm1),col=1)
title(paste( "scale free R^2=",as.character(round(summary(lm1)$adj.r.squared,2)),", slope=", round(lm1$coefficients[[2]],2)))

R平方达到了0.81,已经很不错了。
但是,我们必须有所选择,所以我们可以几个幂次一起算,然后来选就行了。

把以上结果写一个函数

mypick <- function(powerVector,datExpr){
  power <- powerVector
  cor<-stats::cor
  ADJ=abs(cor(datExpr,use = 'p'))^power
  k=apply(ADJ,2,sum) -1
  cut1=cut(k,10)
  binned.k=tapply(k,cut1,mean)
  freq1=tapply(k,cut1,length)/length(k)
  xx= as.vector(log10(binned.k))
  lm1=lm(as.numeric(log10(freq1+.000000001))~ xx )
  return(data.frame(Power=power,
                    SFT.R.sq=as.character(round(summary(lm1)$adj.r.squared,2)),
                    slope=round(lm1$coefficients[[2]],2),
                    mean.k=mean(k)))
}

测试一个结果,幂次为10

mypick(10,datExpr)

结果如下,符合预期,其中mean.k,是对所有基因的连通性取均值,代表当前网络的连通性
之后作图的时候需要用到。

Power SFT.R.sq slope   mean.k
1    10     0.81 -1.66 5.193521

现在批量运算

powers = c(c(1:10), seq(from = 12, to=20, by=2))
do.call(rbind,lapply(powers,mypick,datExpr))

得到结果如下

结合这个表格,我会选取6,作为power值,因为从5到6,R平方显著提升。

官方版本的软阈值计算

以上过程只是帮助我们理解无尺度网络和软阈值的概念。
实际预算,WGCNA包中提供了一个函数pickSoftThreshold,可以轻松计算。
这个函数输入的就是不同的power值和表达矩阵

powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers)

因为他使用了分块还有并行化的思想,所以计算速度十分快。

结果返回的是个列表,里面有两个内容
第一个是,他自己确定的软阈值,该函数如果发现了R平法大于0.85的power值,就返回最小的那个。
这里返回的是6
可以用这句命令查看

sft$powerEstimate

第二个返回的就是上面的表格,

sft$fitIndices

我们之前确定6是靠眼睛看,但是很不直观,所以,抽取表格内的数据作图,就会很方便。

sizeGrWindow(9, 5)
par(mfrow = c(1,2))
cex1 = 0.85
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity"type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

这张图是很常见,由两张图组成,都是以不同的软阈值作为横坐标。
第一张图,纵坐标是R平方。横线画在了0.85的地方
从图上看,软阈值为6的时候,R平方第一次有了突破,达到了0.9.此时网络已经符合无尺度分布。
第二张图,纵坐标是连通性的平均值,我们已经计算过,他会越来越小。
这是必然的规律,联合第一张看就行了。

一些差别

如果你仔细看,会发现,我们自己算的R平方和官方的R平法还有差别。

这是为什么呢?
是因为每一次线性模型的计算,都会返回两个R平方。

summary(lm1)

我们提取的是矫正后的0.81,而他提取的是没有矫正的0.83,如果我们想改,也十分容易。
下面的语句就可以提取。

summary(lm1)$r.squared

一些疑惑

把我们的数据变成无尺度是个人为的事情??
是的。
作者也这么说了,真实基因之间的关系应该是符合无尺度分布的,但是相关性计算出来的不具有代表性。
因为相关不是因果嘛。

那么,我们就认为给他一个power值,把他们的之间的关系变成符合无尺度分布就行了。这是作者的原话。

last but not least

基因的相关性有个选项,默认是unsigned,还可以改成sined
有什么区别?
如果是unsigned,我们计算的时候,无论相关性是正和负的,都取的绝对值。
最终得到的基因其实有正相关和负相关,我们一视同仁了。

但是真实情况又不是这样的。为了权衡,就有了signed方法来削减负相关的影响。
就是取0.5之后,再加上0.5,这样如果是负的,最终还是会变成正的值,只是相关性变小了。
如果是正的,那么相关性会变大。最终都会在0到1的区间内,没有负值。
而作者推荐的是signed,他觉得这更符合真实情况。
后面得到的模块里面的也都是正相关的基因。

当然事情并不绝对,他说,如果是单细胞的数据,那么就是因为一群基因正负调控得到的结果,这时候unsigned的比较合适。
这是一个生物学的问题,不是统计学的问题。
可以在果子学生信公众号回复“果子WGCNA”自助获取作者接近2个小时的演讲视频和PPT,自己感受和体会一下。

现在再看这个图就好理解了吧

假如你的数据最后没有合适的阈值怎么办呢?
作者说了,如果是unsigned的就选6,signed就选12,不要纠结了。

接下来我们还要讲讲模块的获取过程,eigengene值的计算。集结完毕后,作为连同多数据的WGCNA一起更新到答疑群中。

生物信息学

图形化生物软件专题(4):MEGAN

2020-6-21 15:03:20

生物信息学

用Tax4Fun2对16S微生物组数据进行功能预测

2020-6-21 22:13:07

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