網(wǎng)絡圖
在傳統(tǒng)的網(wǎng)絡圖里面懦砂，一般分為非權(quán)重網(wǎng)絡和權(quán)重網(wǎng)絡
對于非權(quán)重網(wǎng)絡：
非權(quán)重網(wǎng)絡的特點是只能表示兩個點之間是否有關(guān)系，單純的“一刀切”组橄，即只有“有關(guān)系”和“無關(guān)系”這兩種類型荞膘，規(guī)定以后閾值，大于該閾值即為有關(guān)系玉工，若小于這個閾值即為無關(guān)系（用 1 表示有關(guān)系羽资；用 0 表示無關(guān)系）
對于權(quán)重網(wǎng)絡：
權(quán)重網(wǎng)絡的特點是不僅能夠表示兩個點之間是否有關(guān)系，還能表示它們關(guān)系的強弱
WGCNA簡介
WGCNA的出發(fā)點是基于系統(tǒng)的基因表達水平來構(gòu)建一個網(wǎng)絡遵班，目的是顯示出基因間的共表達關(guān)系屠升，那么相似表達模式的基因可能存在共調(diào)控潮改、功能相關(guān)或處于同一通路；即如果某些基因的表達趨勢隨著不同處理之間的變化而有相同的變化趨勢（表達模式）腹暖，那么我們認為這些基因很可能在一個通路上汇在，或者在相互調(diào)控的通路上富集。
從而確定在整個網(wǎng)絡上的核心基因
WGCNA原理
1.建立關(guān)系矩陣
由上圖所得：i 和 j 表示任意的兩個基因微服，那么我們定義其相關(guān)系數(shù)矩陣為S趾疚，其中的相關(guān)系數(shù)為第 i 個基因在各處理中的表達量對應第 j 個基因在各處理中的表達量的相關(guān)系數(shù)
計算出S矩陣后，我們可以發(fā)現(xiàn)以蕴，這個矩陣隨關(guān)于對角線對稱的矩陣糙麦，并且對角線的相關(guān)系數(shù)均為1
2.建立關(guān)系矩陣
對于無權(quán)重網(wǎng)絡：
我們根據(jù)上面的S矩陣，設立個閾值丛肮，比方說設立0.6赡磅，那么S矩陣里面的元素大于0.6時，我們認為是有關(guān)系的宝与，小于0.6時我們認為是無關(guān)系的焚廊。其中對于有無關(guān)系：我們利用 1 表示有關(guān)系；0 表示無關(guān)系
那么正如我們之前所介紹的习劫，非權(quán)重網(wǎng)絡只能表示基因間是否有關(guān)系咆瘟，而無法表示基因間這種關(guān)系的強弱
對于權(quán)重網(wǎng)絡：
其中：
gi，gj 代表不同的基因诽里，β代表權(quán)重袒餐，之所以要引入β作為指數(shù)是為了將區(qū)分度不高的幾個基因的相關(guān)系數(shù)給區(qū)分開
我們定義 αij 為第 i 個基因和第 j 個基因的相關(guān)系數(shù)，并取β為指數(shù)谤狡，經(jīng)過這樣的計算以后灸眼，我們定義連通度ki：
其中 n 表示有n 個基因，即連接度ki表示第 i 個基因和其他基因的α值加和
對于權(quán)重網(wǎng)絡圖墓懂，它滿足無尺度分布焰宣，那么什么是無尺度分布呢？
即連通度高的點所占的比例很少捕仔，而連通度低的點所占的比例很高
該圖表示各連接度ki的一個頻率分布直方圖匕积，那么當我們的基因很多的時候，我們的直方圖區(qū)間就可以無限細分榜跌，從而更趨近于概率密度函數(shù)
基于無尺度分布的假設闸天，我們認為p(ki)與ki呈負相關(guān)關(guān)系，為了方便計算斜做，我們定義：
由此，我們利用一元線性回歸取匹配最佳β值湾揽，即用不同的β值去試驗瓤逼，尋找最佳的β值
在線性回歸中笼吟，我們要求 R^2 大于0.8，slope位于 -1 左右霸旗，而平均連接度要盡可能大贷帮，所以該例子中 β=6 是最佳值，如下圖所示：
求出最佳β值以后诱告，帶入即可求得αij撵枢，那么由αij構(gòu)成了鄰接矩陣（軟閾值篩選）：
那么該矩陣里面的數(shù)值介于[0,1]之間，數(shù)值的大小可以表示關(guān)系的強弱
依據(jù)上面S矩陣精居，代入β=6計算αij可得：
（保留4位有效數(shù)字）
注意
這里值得注意的是在實際建模中：
net = blockwiseModules(datExpr, power = 14, maxBlockSize = 6000,
TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs = TRUE,
saveTOMFileBase = "AS-green-FPKM-TOM",
verbose = 3)
TOMType這個參數(shù)表示的意思是：
無向網(wǎng)絡 (unsigned) 的邊屬性計算方式為：abs(cor(genex, geney)) power锄禽；
有向網(wǎng)絡 (signed) 的邊屬性計算方式為： (0.5 * (1+cor(genex, geney)))power;
sign hybrid的邊屬性計算方式為：cor(genex, geney)power : if cor>0 else 0
distance邊屬性計算方式為：(1-(dist/max(dist))2)power
3.TOM矩陣建立
Ω為TOM矩陣
我們基于αij來計算TOM矩陣，從而將基因進行分模塊化靴姿，即哪些基因在一個模塊內(nèi)沃但，哪些基因在另一個模塊內(nèi)
其中 Iij 的目的是引入第三個基因u，即若第 i 個基因和第 j 個基因有連通性佛吓，而第 u 個基因分別于第 i 個基因和第 j 個基因都有連通性宵晚，那么我們認為第 i 個基因和第 j 個基因的連通性通過第 u 個基因的得到了間接加強，所以维雇，當兩個基因連通性很高的時候淤刃，我們通過引入 Iij 來尋找間接使這兩個基因連通性加強的某些基因，如下圖所示：
計算出TOM矩陣后吱型，我們就可以將基因進行聚類分模塊了逸贾，那么我們定義第 i 個基因和第 j 個基因的聚類距離為：
因此根據(jù)這個聚類距離，我們可以得到分模塊的情況：
那么分好模塊以后唁影，我們可以計算模塊與性狀矩陣的相關(guān)性：
其中g(shù)rey模塊是unassigned genes模塊耕陷，其他顏色的模塊是聚類好的模塊
Module-trait relationships這副圖里面模塊與性狀的相關(guān)性（上面的數(shù)值為相關(guān)系數(shù)，下面數(shù)值為顯著性）
其中据沈，性狀矩陣即為每一種處理
而每一模塊的基因相當于整個基因表達譜的子集哟沫，對于整個的基因表達譜來說：
對于某一模塊的基因表達譜：
假設說某一模塊有200個基因，那么怎么計算模塊與性狀之間的相關(guān)性呢锌介？
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0); ##不同顏色的模塊的ME值矩陣(性狀[樣本]vs模塊)
moduleTraitCor = cor(MEs, design , use = "p"); ##這里的design為性狀矩陣（樣本矩陣）
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
這里的性狀1 - 性狀n即為下圖的 D1 - Dn
由代碼可知嗜诀，MEs是每一個模塊中性狀的特征值，該特征值即為在某一模塊中孔祸，對于某一個性狀隆敢，將其基因表達量PCA分解，選取PC1作為該性狀的特征值崔慧，如下圖：
具體為：
橫坐標表示模塊拂蝎，縱坐標表示性狀以及生物學重復，表中的元素表示每個模塊下各個性狀的PC1的值惶室，那么總體上模塊與性狀的相關(guān)性為cor(MEs, design , use = "p")温自，其中design為性狀矩陣：
那么通過這兩個矩陣玄货，我們就可以計算出模塊與性狀的相關(guān)性了：
其中縱坐標為模塊，橫坐標為性狀
這樣一來我們計算出了Module-trait relationships這副圖里面模塊與性狀的相關(guān)性了（上面的數(shù)值為相關(guān)系數(shù)悼泌，下面數(shù)值為顯著性）
4.GS與MM
GS是Gene Signicance松捉，描述的是某一個基因與性狀的相關(guān)性，即某基因（在基因表達矩陣中馆里，該基因隨性狀的變化的向量與模塊特征向量的相關(guān)性（即上一小節(jié)中隘世，MEs其中的某一列即為某個模塊的特征向量，該特征向量不是n階方陣的特征向量）鸠踪；
而MM是Module Membership丙者，描述的是某一個基因與模塊之間的相關(guān)性，即某基因（在基因表達矩陣中慢哈，該基因隨性狀的變化的向量蔓钟，即基因表達譜對應該基因的那一行）與該模塊特征向量的相關(guān)性（即上一小節(jié)中，MEs其中的某一列即為某個模塊的特征向量卵贱，該特征向量不是n階方陣的特征向量）
那么我們怎么照hub基因呢滥沫？
通常需要做出GS和MM的關(guān)系圖
那么如果某個模塊內(nèi)基因向量與某性狀向量的相關(guān)性高，并且這個基因向量與該模塊的PC1向量相關(guān)性高键俱，則這個基因很可能為hub基因兰绣；如上圖中右上角的點代表的基因即為hub基因
基于我們之前說到的MEs矩陣，我們可以利用PC1來看下各個sample之間特征值的變化情況
barplot(as.matrix(t(MEsWW$MEtan)), col='red', main="", cex.main=2,
ylab="eigengene expression",xlab="array sample")
上面的條形圖编振，每一根柱子代表每一個sample（D1_rep1缀辩，D1_rep2，D2_rep1踪央，D2_rep2等等）在對應模塊（比方說MEtan模塊）的PC1值
我們具體看一下各個元素代表什么：
假設MEtan模塊有200個基因臀玄，那么MEtan模塊中的D1_rep1對應的0.47863960代表在D1_rep1中該模塊所有基因的表達特征值，由于劃分模塊以后畅蹂，我們認為在同一模塊里面的基因是共表達的健无，所以某一個模塊內(nèi)對應sample的PC1值（PC1涵蓋整個基因表達譜的最大信息量）可以衡量該模塊內(nèi)這個sample的整體表達特征
5.cytoscope網(wǎng)絡圖
接下來，我們就可以將自己感興趣的液斜，與性狀相關(guān)性比較大的模塊基因?qū)С隼巯停胏ytoscope查看基因間作用關(guān)系圖了。
幾個注意事項
1.input文件
對于WGCNA來說少漆，input數(shù)據(jù)無非就兩種臼膏，一個是基因表達矩陣（一般用FPKM表示），另一個是性狀矩陣
上表中列為性狀（樣品名稱）示损，行為基因名稱
是否有某性狀可以利用0（無）渗磅，1（有）來區(qū)分
2.篩選基因
一般來說，我們盡量保留基因表達矩陣里面的全部基因，但是在WGCNA做基因聚類的時候始鱼，由于計算量過大论巍，對于某些物種來說，可能要進行取舍
但是根據(jù)官網(wǎng)上的提示风响，并不建議利用差異基因來用作篩選條件
這是由于利用差異基因來做可能會破壞之前提到的無尺度分布這個條件。并且對于某些項目丹禀，其核心基因并不一定是差異基因状勤，可能某些核心基因在各性狀組別中表達差異不大，但是很可能參與到多個通路中去双泪，從而影響生物學過程持搜，而差異表達的基因往往都比較偏下游，利用差異表達基因作為篩選條件焙矛，可能會miss掉很多核心基因
所以葫盼，篩選的條件根據(jù)項目的不同而靈活選擇，除了篩選不表達的基因以外村斟。這里還有有一種篩選方法贫导，即在基因表達矩陣中，對每一行基因計算其方差蟆盹，我們盡量選擇方差排序較大的基因孩灯，這樣做出來的效果可能會好一些。
參考：
http://www.reibang.com/p/2d3d079cc338
https://www.omicsclass.com/article/649
《WGCNA: an R package for weighted correlation network analysis》
《Weighted Network analysis》
