使用WGCNA进行分析时，一定会遇到这两幅图，其中有这么几个重要的概念，软阈值，无尺度网络分布。
使用WGCNA分析时，这两幅图肯定会使用到，这两幅图的主要价值就在于寻找软阈值，构建无尺度分布网络。
image.png
无尺度分布
网络中的节点服从幂律分布，叫做无尺度网络分布。为什么是这样呢？
二八法则 就是一种幂律分布，80%的人占据着不到20%的财富，而极少数的人占据着大量的财富。一些少数关键的节点和其他节点有着丰富的联系，而大部分的节点与其他节点联系少，这些关键节点就叫作hub 节点，枢纽节点。
下面这个图，形象的表现了这个关系。钟形曲线是正态分布，幂律分布就是少数节点有极端的节点连接数。
随机网络是可以用一个值来衡量大多数节点之间的距离，或者是6或者是60，也就是你可以用一个尺度去衡量他，可以称为尺度网络。
而有巨人的网络，没办法来衡量两个节点之间的距离，叫作无尺度分布，该分布又叫做幂律分布，我们说的二八定律，长尾定律都是幂律分布的口头化呈现。
生物学相关的连接是符合这种无尺度网络分布的，也就是幂律分布，这样分布的好处在于少了其中的一些节点，这个网络依然能够正常的运转。
真实的生物学网络比如：
human disease network，
gene co-expression networks，
protein-protein interaction networks，
cell-cell interaction networks，
the world wide web and social interaction networks
生物体内，对于一个细胞而言，其中表达了很多蛋白或者RNA，但是不是每一个的作用都一样，其中有一些起到了极端重要的作用，而有一些作用甚微。对于一个微环境中，有一些细胞可能对于微环境的功能起到很重要的作用，某一些作用一般。
这种生物学网络就符合无尺度网络分布
WGCNA就是利用了这样的规律，让基因之间的联系符合无尺度分布。
软阈值
WGCNA全称是Weighted Gene Co-Expression Network Analysis，加权基因共表达网络分析。网络中基因与基因是否连接，取决于是否发生共表达。
真实网络中，计算共表达的方法有很多种，WGCNA官网给出的方法是Pearson相关性分析，我们通过计算得到两个基因的相关性系数后，如何确定两个基因在网络中连边呢。
这个时候，有两种思路，直接按照差异分析的方法，确定一个相关性系数，比如0.7，大于这个相关性系数的，我们就认为共表达，小于这个的，我们就认为这是非共表达。这个时候选用0.7作为分界点一刀切，就是"hard threshold"，即就是不加权。
这种方法的有点就在于，简单直接。
但缺点也很明显：
怎样确定这个阈值，把0.7作为分界点，那0.69的相关性怎么办？
真实的生物学网络，这种非黑即白的定义方法不符合实际。就如同上面所说，大多数生物学网络符合说明无尺度网络分布。
因此，这个时候就引入了软阈值这个概念。
那么怎么确定这个软阈值。
计算软阈值
WGCNA有直接计算软阈值的函数，可以把函数拆解开，自己学习一下。
载入数据
载入WGCNA官方的数据，行是样本名称，列是基因名称
WGCNA首先使用相关性分析，上文中提到Pearson相关性算法，计算任意两个基因之间的相关性。为了避免上文中提到直接使用hard threshold的缺陷，WGCNA给相关性加了一个幂次。
图片来自WGCNA的分析中为什么要挑选软阈值
图中的绿线原来的相关性，进行幂次计算之后(8次)就变成了红色。
可以看到整体的相关性都变小了，但是里面本来就小的变得更加小。
我其实对这个图还是没有特别的理解，相关性是横坐标，纵坐标是邻接性，和得出的结论似乎不太符合。留作以后的思考吧。
为什么要给相关性加上幂次运算呢，这样就使得更符合幂律分布，符合无尺度网络的特点了。
先计算基因两两之间的相关性
cor_data <- cor(datExpr,use = 'p', method = 'pearson')
cor_data[1:5,1:5]
这是采用R自带的函数进行计算，得到的结果是
> cor_data[1:5,1:5]
MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080
MMT00000044 1.00000000 -0.03797242 0.09326526 0.24716099 0.13636580
MMT00000046 -0.03797242 1.00000000 -0.56393957 -0.02934881 -0.06291723
MMT00000051 0.09326526 -0.56393957 1.00000000 0.06120014 -0.05422173
MMT00000076 0.24716099 -0.02934881 0.06120014 1.00000000 -0.02327867
MMT00000080 0.13636580 -0.06291723 -0.05422173 -0.02327867 1.00000000
然后使用WGCNA包的相关性算法计算
cor_data_WGCNA <- WGCNA::cor(datExpr,use = 'p')
cor_data_WGCNA[1:5,1:5]
结果是这样
> cor_data_WGCNA[1:5,1:5]
MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080
MMT00000044 1.00000000 -0.03797242 0.09326526 0.24716099 0.13636580
MMT00000046 -0.03797242 1.00000000 -0.56393957 -0.02934881 -0.06291723
MMT00000051 0.09326526 -0.56393957 1.00000000 0.06120014 -0.05422173
MMT00000076 0.24716099 -0.02934881 0.06120014 1.00000000 -0.02327867
MMT00000080 0.13636580 -0.06291723 -0.05422173 -0.02327867 1.00000000
两种计算的结果是一样的，也就说明WGCNA的算法采用的是Pearson相关性算法。
我在这里试验了一下使用spearman算法，耗时比较长。应该是spearman算法计算相关性时，每次计算需要将每组的顺序进行排列，这个时候就会很耗费算力。
给相关性加一个幂次
为什么要加一个幂次，作者就是通过加上幂次的方式，强行让基因之间的连通性开始符合幂律了。
果子试验的时候是以10次幂进行，我选择8次幂看看效果。
power <- 8
ADJ=abs(cor(datExpr,use = 'p'))^power
这样，将相关性矩阵绝对值之后，都给了一个8次幂，结果如下
> dim(ADJ)
[1] 3600 3600
> ADJ[1:5,1:5]
MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080
MMT00000044 1.000000e+00 4.322615e-12 5.724787e-09 1.392642e-05 1.195759e-07
MMT00000046 4.322615e-12 1.000000e+00 1.022965e-02 5.504572e-13 2.455594e-10
MMT00000051 5.724787e-09 1.022965e-02 1.000000e+00 1.967975e-10 7.471142e-11
MMT00000076 1.392642e-05 5.504572e-13 1.967975e-10 1.000000e+00 8.623133e-14
MMT00000080 1.195759e-07 2.455594e-10 7.471142e-11 8.623133e-14 1.000000e+00
然后对每一列的基因进行求和，得到的就是这个基因跟其他基因相关性之和。减去1是排除了自身。
k=apply(ADJ,2,sum) -1
hist(k)
结果如下
> k[1:10]
MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080 MMT00000102 MMT00000149
0.0814313 15.6397241 10.3721409 0.5134981 11.9155349 3.1397842 16.5059843
MMT00000159 MMT00000207 MMT00000212
15.6716271 22.4635281 1.0466731
横坐标是基因的连接度，纵坐标是拥有该连接点数的频率，大体看上去，是头大尾巴长的幂律分布。
但是，符合与否，还要有相应的检验才可以，具体如何确定一个数据是否符合幂律分布，可看本片文章
这里大致说明，把连通性分隔，分隔内连通性的平均值取log10，跟频率的概率取log10，两者之间有线性关系。为什么是这样，这副图做了一个说明
图片来自无标度网络幂律分布
那么大致画图看一下
先把k从小到大排序，切割成20份
cut1=cut(k,20)
计算每个区间的平均值
binned.k=tapply(k,cut1,mean)
binned.k
> binned.k
(-0.0728,3.82] (3.82,7.63] (7.63,11.5] (11.5,15.3] (15.3,19.1]
2.094294 5.583894 9.466252 13.116735 17.016110
(19.1,22.9] (22.9,26.7] (26.7,30.5] (30.5,34.3] (34.3,38.2]
20.917832 24.597606 28.535248 31.739992 36.147825
(38.2,42] (42,45.8] (45.8,49.6] (49.6,53.4] (53.4,57.2]
40.351527 43.414818 48.002510 51.974849 55.419944
(57.2,61.1] (61.1,64.9] (64.9,68.7] (68.7,72.5] (72.5,76.4]
59.009976 62.642794 66.517636 70.593310 74.597006
然后计算每个区间的频率
freq1=tapply(k,cut1,length)/length(k)
freq1
> freq1
(-0.0728,3.82] (3.82,7.63] (7.63,11.5] (11.5,15.3] (15.3,19.1]
0.392222222 0.224722222 0.137222222 0.083333333 0.051388889
(19.1,22.9] (22.9,26.7] (26.7,30.5] (30.5,34.3] (34.3,38.2]
0.032777778 0.016944444 0.016111111 0.008611111 0.005277778
(38.2,42] (42,45.8] (45.8,49.6] (49.6,53.4] (53.4,57.2]
0.001111111 0.003055556 0.002222222 0.002777778 0.001388889
(57.2,61.1] (61.1,64.9] (64.9,68.7] (68.7,72.5] (72.5,76.4]
0.005555556 0.003611111 0.003888889 0.004722222 0.003055556
画图查看
# 查看均值的对数和频率的对数是否为线性
plot(log10(binned.k),log10(freq1+.000000001),xlab="log10(k)",ylab="log10(p(k))")
将分割范围确定为20的情况
将分割范围确定为10的情况
将分割范围确定为30的情况
无限的分割，就越接近于单个个体
从上面的图中可以大致看出，呈现一个线性的趋势，那么通用线性函数加线和注释，如下所示（也就是检验是否符合线性关系）。
这个时候就是将问题转化为了查看 和 是否符合线性关系。
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=2, lty = 1, lwd = 3)
title(paste( "scale free R^2=",as.character(round(summary(lm1)$adj.r.squared,2)), " slope=", round(lm1$coefficients[[2]],2)))
提取这个检验中的内容
> summary(lm1)
Call:
lm(formula = as.numeric(log10(freq1 + 1e-09)) ~ xx)
Residuals:
Min 1Q Median 3Q Max
-0.8811 -0.1809 0.1189 0.2193 0.3454
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2989 0.2228 1.341 0.191
xx -1.6776 0.1477 -11.360 5.36e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3266 on 28 degrees of freedom
Multiple R-squared: 0.8217, Adjusted R-squared: 0.8154
F-statistic: 129.1 on 1 and 28 DF, p-value: 5.362e-12
这里我们提取的是调整后的R^2
探究不同的次幂的影响
上面是对采用8次幂的结果的观察，那么如果采用其他次幂是什么样子
构建一个计算函数
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)$r.squared,2)),
slope=round(lm1$coefficients[[2]],2),
mean.k=mean(k)))
}
仔细看上面的函数代码，其实是把上面单个的步骤，集合到了一起，写好函数，可以用循环思维，求取其他值。上面的代码，将R^2的从Adjusted R-squared更换为 R-squared，因为WGCNA使用的是 R-squared，而不是Adjusted R-squared。
具体函数结果要返回什么值，主要取决于自己的画图需求
测试一下函数可用否
mypick(6,datExpr)
结果
> mypick(6,datExpr)
Power SFT.R.sq slope mean.k
1 6 0.89 -1.5 19.90579
然后批量计算出其他阈值的情况
powers = c(c(1:10), seq(from = 12, to=20, by=2))
do.call(rbind,lapply(powers,mypick,datExpr))
结果如下
> do.call(rbind,lapply(powers,mypick,datExpr))
Power SFT.R.sq slope mean.k
1 1 -0.09 0.34 746.977852
2 2 0.02 -0.60 254.497776
3 3 0.26 -1.03 111.004761
4 4 0.44 -1.42 56.535969
5 5 0.64 -1.72 32.154311
6 6 0.89 -1.50 19.905790
7 7 0.91 -1.67 13.192338
8 8 0.89 -1.72 9.248524
9 9 0.84 -1.70 6.795315
10 10 0.81 -1.66 5.193521
11 12 0.84 -1.48 3.332680
12 14 0.86 -1.38 2.348905
13 16 0.9 -1.30 1.767929
14 18 0.9 -1.24 1.393515
15 20 0.92 -1.21 1.135087
结合这个表格，一般会选取6，作为power值，因为从5到6，R^2显著提升。也就是所说的软阈值
代码是从公众号摘录的，其中使用到的批量思维do.call + lapply + function，确实值得细细品味，并加以应用。循环的思维在分析中，能够提高效率。
WGCNA官方版本
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers)
根据果子学生信所说，WGCNA采用了分块化和并行方法，所以速度会快很多。
结果如下：
> powers = c(c(1:10), seq(from = 12, to=20, by=2))
> sft = WGCNA::pickSoftThreshold(datExpr, powerVector = powers)
Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
1 1 0.0278 0.345 0.456 747.00 762.0000 1210.0
2 2 0.1260 -0.597 0.843 254.00 251.0000 574.0
3 3 0.3400 -1.030 0.972 111.00 102.0000 324.0
4 4 0.5060 -1.420 0.973 56.50 47.2000 202.0
5 5 0.6810 -1.720 0.940 32.20 25.1000 134.0
6 6 0.9020 -1.500 0.962 19.90 14.5000 94.8
7 7 0.9210 -1.670 0.917 13.20 8.6800 84.1
8 8 0.9040 -1.720 0.876 9.25 5.3900 76.3
9 9 0.8590 -1.700 0.836 6.80 3.5600 70.5
10 10 0.8330 -1.660 0.831 5.19 2.3800 65.8
11 12 0.8530 -1.480 0.911 3.33 1.1500 58.1
12 14 0.8760 -1.380 0.949 2.35 0.5740 51.9
13 16 0.9070 -1.300 0.970 1.77 0.3090 46.8
14 18 0.9120 -1.240 0.973 1.39 0.1670 42.5
15 20 0.9310 -1.210 0.977 1.14 0.0951 38.7
观察上面的结果，依然是5到6的时候，R^2的值会显著变化
sft$powerEstimate
> sft$powerEstimate
[1] 6
该函数如果发现了R平法大于0.85的power值，就返回最小的那个。
这里返回的是6
为了让结果直观，WGCNA文档把这个图也画出来，
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^2。横线画在了0.85的地方
从图上看，软阈值为6的时候，R平方第一次有了突破，达到了0.9，此时网络已经符合无尺度分布。
第二张图，纵坐标是连通性的平均值，我们已经计算过，他会越来越小。（也就是和基因之间的关联性的平均值，又加上了次幂的结果），必然会越来越小。
思考
WGCNA作者是通过基因相关性，认为使用次幂的形式，将这样一个相关性网络符合幂律形式，进而后续分析。但是涉及到的将相关性的取绝对值等一系列问题，是否取绝对值，这不是数学问题，而是生物学问题，更需要结合具体的生物学问题取舍和证明。
WGCNA: an R package for weighted correlation network analysis | BMC Bioinformatics
WGCNA的分析中为什么要挑选软阈值
小鬼的WGCNA图文详解(一)--软阈值