【WGCNA学习笔记】无尺度分布和软阈值 |
您所在的位置：网站首页 › 阈值是怎么确定的 › 【WGCNA学习笔记】无尺度分布和软阈值 |
【WGCNA学习笔记】无尺度分布和软阈值
|
使用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_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 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 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的情况
无限的分割，就越接近于单个个体 从上面的图中可以大致看出，呈现一个线性的趋势，那么通用线性函数加线和注释，如下所示（也就是检验是否符合线性关系）。 这个时候就是将问题转化为了查看 
提取这个检验中的内容 > 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 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作者是通过基因相关性，认为使用次幂的形式，将这样一个相关性网络符合幂律形式，进而后续分析。但是涉及到的将相关性的取绝对值等一系列问题，是否取绝对值，这不是数学问题，而是生物学问题，更需要结合具体的生物学问题取舍和证明。
|
【本文地址】