1.邊的理解：WGCNA分析全稱是加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析，而在基因調(diào)控網(wǎng)絡(luò)中兩個重要的概念就是點(diǎn)和邊。點(diǎn)代表基因，邊代表基因表達(dá)相關(guān)性。
2.關(guān)于軟閾值的理解：
①計(jì)算共表達(dá)有很多方法，官網(wǎng)教程中給出的方法是用皮爾森相關(guān)系數(shù)(Pearson) correlation coefficient。加權(quán)是指對相關(guān)性值進(jìn)行冥次運(yùn)算 (冥次的值也就是軟閾值 (power, pickSoftThreshold這個函數(shù)所做的就是確定合適的power))
②基因調(diào)控網(wǎng)絡(luò)更符合無尺度網(wǎng)絡(luò)（scale-free network）：在生物體調(diào)控網(wǎng)絡(luò)中，有少數(shù)的一些基因起到非常重要的一些調(diào)控作用，而其它基因沒有它們的調(diào)控層次高，更符合無尺度網(wǎng)絡(luò)。
③無尺度網(wǎng)絡(luò)的實(shí)現(xiàn)：在無尺度網(wǎng)絡(luò)中，兩個基因之間相關(guān)性的計(jì)算方式為abs(cor(geneX,geneY))^power，冪運(yùn)算方式強(qiáng)化了強(qiáng)相關(guān)，弱化了弱相關(guān)。使得相關(guān)性數(shù)值更符合無尺度網(wǎng)絡(luò)特征，更具有生物意義。解釋：正常情況下，我們可以定義一個閾值比如0.8，大于0.8叫有相關(guān)性，小于0.8叫無相關(guān)性，而這時(shí)0.799的基因可能就會出來干擾，因此可通過給相關(guān)性加冪的方法解決。這樣就會使得大的更大，而小的更小。
3.挑選軟閾值的步驟拆解如下：
以冪為10來計(jì)算一下：
power <- 10
ADJ <- abs(cor(matrix_data))^power
現(xiàn)在就得到了一個基因之間相關(guān)性的矩陣，
現(xiàn)在對每一列的基因進(jìn)行求和，得到的就是這個基因跟其它基因的相關(guān)性之和。減去1就是排除了自身。
k=apply(ADJ, 2, sum)-1
這個K指的就是3600基因組成的網(wǎng)絡(luò)中，每個節(jié)點(diǎn)的連通性。做一個頻次分布圖如下：
橫坐標(biāo)是連通性依次升高，縱坐標(biāo)表示該范圍連通性的頻次，那么究竟符不符合冪律分布還要計(jì)算。
真正的冪律是這樣的，把連通性分隔，分隔內(nèi)連通性的平均值取log10，跟頻率的概率取log10，兩者之間有線性關(guān)系。
先把k從小到大排序，切割成10份
cut1=cut(k,10)
計(jì)算每個區(qū)間的平均值
binned.k=tapply(k,cut1,mean)
然后計(jì)算每個區(qū)間的頻率
freq1=tapply(k,cut1,length)/length(k)
此時(shí)這個均值的對數(shù)和頻率的對數(shù)就是線性的
plot(log10(binned.k),log10(freq1+.000000001),xlab="log10(k)",ylab="log10(p(k))")
如果通用線性函數(shù)加線和注釋就會明顯一點(diǎn)
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平方達(dá)到了0.83
接下來我們可以幾個冪次一起算，然后來選就行了。
把以上結(jié)果寫一個函數(shù)
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)))
測試一個結(jié)果,冪次為10
mypick(10,matrix_data)
結(jié)果如下：
Power SFT.R.sq slope mean.k
10 0.83 -1.98 1.608301
其中mean.k，是對所有基因的連通性取均值，代表當(dāng)前網(wǎng)絡(luò)的連通性
之后作圖的時(shí)候需要用到。
現(xiàn)在批量運(yùn)算
powers = c(c(1:10), seq(from = 12, to=20, by=2))
do.call(rbind,lapply(powers,mypick,matrix_data))
得到結(jié)果如下：
結(jié)合這個表格，我會選取6，作為power值。
4.官方版本的軟閾值計(jì)算
實(shí)際上，WGCNA包中提供了一個函數(shù)pickSoftThreshold,可以幫助我們一鍵計(jì)算。
代碼如下：
powers = c(seq(1,10,by = 1),seq(12,20,by = 2))
sft = pickSoftThreshold(matrix_data,powerVector = powers,verbose = 5)
第一個是，他自己確定的軟閾值，該函數(shù)如果發(fā)現(xiàn)了R平法大于0.85的power值，就返回最小的那個,用下面代碼查看：
power = sft$powerEstimate
power
第二個返回的就是上面的表格：
sft$fitIndices
