1.边的理解：WGCNA分析全称是加权基因共表达网络分析，而在基因调控网络中两个重要的概念就是点和边。点代表基因，边代表基因表达相关性。
2.关于软阈值的理解：
①计算共表达有很多方法，官网教程中给出的方法是用皮尔森相关系数(Pearson) correlation coefficient。加权是指对相关性值进行冥次运算 (冥次的值也就是软阈值 (power, pickSoftThreshold这个函数所做的就是确定合适的power))
②基因调控网络更符合无尺度网络（scale-free network）：在生物体调控网络中，有少数的一些基因起到非常重要的一些调控作用，而其它基因没有它们的调控层次高，更符合无尺度网络。
③无尺度网络的实现：在无尺度网络中，两个基因之间相关性的计算方式为abs(cor(geneX,geneY))^power，幂运算方式强化了强相关，弱化了弱相关。使得相关性数值更符合无尺度网络特征，更具有生物意义。解释：正常情况下，我们可以定义一个阈值比如0.8，大于0.8叫有相关性，小于0.8叫无相关性，而这时0.799的基因可能就会出来干扰，因此可通过给相关性加幂的方法解决。这样就会使得大的更大，而小的更小。
3.挑选软阈值的步骤拆解如下：
以幂为10来计算一下：
power <- 10
ADJ <- abs(cor(matrix_data))^power
现在就得到了一个基因之间相关性的矩阵，
现在对每一列的基因进行求和，得到的就是这个基因跟其它基因的相关性之和。减去1就是排除了自身。
k=apply(ADJ, 2, sum)-1
这个K指的就是3600基因组成的网络中，每个节点的连通性。做一个频次分布图如下：
横坐标是连通性依次升高，纵坐标表示该范围连通性的频次，那么究竟符不符合幂律分布还要计算。
真正的幂律是这样的，把连通性分隔，分隔内连通性的平均值取log10，跟频率的概率取log10，两者之间有线性关系。
先把k从小到大排序，切割成10份
cut1=cut(k,10)
计算每个区间的平均值
binned.k=tapply(k,cut1,mean)
然后计算每个区间的频率
freq1=tapply(k,cut1,length)/length(k)
此时这个均值的对数和频率的对数就是线性的
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.83
接下来我们可以几个幂次一起算，然后来选就行了。
把以上结果写一个函数
mypick <- function(powers,matrix_data){
power <- powers
cor<-stats::cor
ADJ=abs(cor(matrix_data))^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,matrix_data)
结果如下：
Power SFT.R.sq slope mean.k
10 0.83 -1.98 1.608301
其中mean.k，是对所有基因的连通性取均值，代表当前网络的连通性
之后作图的时候需要用到。
现在批量运算
powers = c(c(1:10), seq(from = 12, to=20, by=2))
do.call(rbind,lapply(powers,mypick,matrix_data))
得到结果如下：
结合这个表格，我会选取6，作为power值。
4.官方版本的软阈值计算
实际上，WGCNA包中提供了一个函数pickSoftThreshold,可以帮助我们一键计算。
代码如下：
powers = c(seq(1,10,by = 1),seq(12,20,by = 2))
sft = pickSoftThreshold(matrix_data,powerVector = powers,verbose = 5)
第一个是，他自己确定的软阈值，该函数如果发现了R平法大于0.85的power值，就返回最小的那个,用下面代码查看：
power = sft$powerEstimate
power
第二个返回的就是上面的表格：
sft$fitIndices
