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

WGCNA是个强悍而且万金油的分析数据分析技能。
我们之前讲过他的常见强悍用法
视频课程：WGCNA，多样本，多分组，多临床信息的数据挖掘利器。

还讲过他的万金油用法：
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一起更新到答疑群中。