WGCNA基本概念

加權(quán)相關(guān)網(wǎng)絡(luò)分析(WGCNA)是一種用于描述不同樣本間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法，可用于識(shí)別高度協(xié)同的基因集，并根據(jù)基因集的互連性和基因集與表型的相關(guān)性識(shí)別候選生物標(biāo)記基因或治療靶點(diǎn)。

與只關(guān)注差異表達(dá)基因相比，WGCNA利用數(shù)千或近萬(wàn)個(gè)變化最大的基因或全部基因的信息來(lái)識(shí)別感興趣的基因集，并與表型進(jìn)行顯著相關(guān)分析。一是充分利用信息；第二，將數(shù)千個(gè)基因和表型之間的關(guān)聯(lián)轉(zhuǎn)化為幾個(gè)基因集和表型之間的關(guān)聯(lián)，從而避免了多重假設(shè)檢驗(yàn)和校正的問(wèn)題。

要了解WGCNA，我們需要了解WGCNA中的以下術(shù)語(yǔ)及其定義。

共表達(dá)網(wǎng)絡(luò)：定義為加權(quán)基因網(wǎng)絡(luò)。點(diǎn)代表基因，邊代表基因表達(dá)相關(guān)性。加權(quán)是指對(duì)相關(guān)性值進(jìn)行冥次運(yùn)算(冥次的值也就是軟閾值 (power,pickSoftThreshold這個(gè)函數(shù)所做的就是確定合適的power))。無(wú)向網(wǎng)絡(luò)的邊屬性計(jì)算方式為abs(cor(genex, geney)) ^ power；有向網(wǎng)絡(luò)的邊屬性計(jì)算方式為(1+cor(genex, geney)/2) ^ power; signhybrid的邊屬性計(jì)算方式為cor(genex, geney)^power if cor>0 else 0。這種處理方式強(qiáng)化了強(qiáng)相關(guān)，弱化了弱相關(guān)或負(fù)相關(guān)，使得相關(guān)性數(shù)值更符合無(wú)標(biāo)度網(wǎng)絡(luò)特征，更具有生物意義。如果沒(méi)有合適的power，一般是由于部分樣品與其它樣品因?yàn)槟撤N原因差別太大導(dǎo)致的，可根據(jù)具體問(wèn)題移除部分樣品或查看后面的經(jīng)驗(yàn)值。Module(模塊)：高度內(nèi)連的基因集。在無(wú)向網(wǎng)絡(luò)中，模塊內(nèi)是高度相關(guān)的基因。在有向網(wǎng)絡(luò)中，模塊內(nèi)是高度正相關(guān)的基因。把基因聚類成模塊后，可以對(duì)每個(gè)模塊進(jìn)行三個(gè)層次的分析：1. 功能富集分析查看其功能特征是否與研究目的相符；2. 模塊與性狀進(jìn)行關(guān)聯(lián)分析，找出與關(guān)注性狀相關(guān)度最高的模塊；3. 模塊與樣本進(jìn)行關(guān)聯(lián)分析，找到樣品特異高表達(dá)的模塊。 基因富集相關(guān)文章 去東方，最好用的在線GO富集分析工具；GO、GSEA富集分析一網(wǎng)打進(jìn)；GSEA富集分析-界面操作。其它關(guān)聯(lián)后面都會(huì)提及。Connectivity (連接度)：類似于網(wǎng)絡(luò)中 “度”(degree)的概念。每個(gè)基因的連接度是與其相連的基因的邊屬性之和。Module eigengene E:給定模型的第一主成分，代表整個(gè)模型的基因表達(dá)譜。這個(gè)是個(gè)很巧妙的梳理，我們之前講過(guò)PCA分析的降維作用，之前主要是拿來(lái)做可視化，現(xiàn)在用到這個(gè)地方，很好的用一個(gè)向量代替了一個(gè)矩陣，方便后期計(jì)算。(降維除了PCA，還可以看看tSNE)Intramodular connectivity:給定基因與給定模型內(nèi)其他基因的關(guān)聯(lián)度，判斷基因所屬關(guān)系。Module membership: 給定基因表達(dá)譜與給定模型的eigengene的相關(guān)性。Hub gene: 關(guān)鍵基因 (連接度最多或連接多個(gè)模塊的基因)。Adjacency matrix(鄰接矩陣)：基因和基因之間的加權(quán)相關(guān)性值構(gòu)成的矩陣。TOM (Topological overlapmatrix)：把鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣，以降低噪音和假相關(guān)，獲得的新距離矩陣，這個(gè)信息可拿來(lái)構(gòu)建網(wǎng)絡(luò)或繪制TOM圖。

基本分析過(guò)程

構(gòu)建基因共表達(dá)網(wǎng)絡(luò)：使用加權(quán)的表達(dá)相關(guān)性。識(shí)別基因集：基于加權(quán)相關(guān)性，進(jìn)行層級(jí)聚類分析，并根據(jù)設(shè)定標(biāo)準(zhǔn)切分聚類結(jié)果，獲得不同的基因模塊，用聚類樹(shù)的分枝和不同顏色表示。如果有表型信息，計(jì)算基因模塊與表型的相關(guān)性，鑒定性狀相關(guān)的模塊。研究模型之間的關(guān)系，從系統(tǒng)層面查看不同模型的互作網(wǎng)絡(luò)。從關(guān)鍵模型中選擇感興趣的驅(qū)動(dòng)基因，或根據(jù)模型中已知基因的功能推測(cè)未知基因的功能。導(dǎo)出TOM矩陣，繪制相關(guān)性圖。

WGCNA包實(shí)戰(zhàn)

R-package WGCNA是一組計(jì)算各種加權(quán)關(guān)聯(lián)分析的函數(shù)，可用于網(wǎng)絡(luò)構(gòu)建、基因篩選、基因簇識(shí)別、拓?fù)涮卣饔?jì)算、數(shù)據(jù)模擬和可視化等。

輸入數(shù)據(jù)和參數(shù)選擇

WGCNA本質(zhì)是基于相關(guān)系數(shù)的網(wǎng)絡(luò)分析方法，適用于多樣品數(shù)據(jù)模式，一般要求樣本數(shù)多于15個(gè)。樣本數(shù)多于20時(shí)效果更好，樣本越多，結(jié)果越穩(wěn)定?；虮磉_(dá)矩陣:常規(guī)表達(dá)矩陣即可，即基因在行，樣品在列，進(jìn)入分析前做一個(gè)轉(zhuǎn)置。RPKM、FPKM或其它標(biāo)準(zhǔn)化方法影響不大，推薦使用Deseq2的varianceStabilizingTransformation或log2(x+1)對(duì)標(biāo)準(zhǔn)化后的數(shù)據(jù)做個(gè)轉(zhuǎn)換。如果數(shù)據(jù)來(lái)自不同的批次，需要先移除批次效應(yīng)(記得上次轉(zhuǎn)錄組培訓(xùn)課講過(guò)如何操作)。如果數(shù)據(jù)存在系統(tǒng)偏移，需要做下quantile normalization。性狀矩陣：用于關(guān)聯(lián)分析的性狀必須是數(shù)值型特征(如下面示例中的Height, Weight,Diameter)。如果是區(qū)域或分類變量，需要轉(zhuǎn)換為0-1矩陣的形式(1表示屬于此組或有此屬性，0表示不屬于此組或無(wú)此屬性，如樣品分組信息WT,KO, OE)。 ID WT KO OE Height Weight Diameter samp1 1 0 0 1 2 3 samp2 1 0 0 2 4 6 samp3 0 1 0 10 20 50 samp4 0 1 0 15 30 80 samp5 0 0 1 NA 9 8 samp6 0 0 1 4 8 7推薦使用Signed network和Robust correlation (bicor)。(這個(gè)根據(jù)自己的需要，看看上面寫(xiě)的每個(gè)網(wǎng)絡(luò)怎么計(jì)算的，更知道怎么選擇)無(wú)向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)，沒(méi)有一個(gè)power值可以使無(wú)標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8或平均連接度降到100以下，可能是由于部分樣品與其他樣品差別太大造成的。這可能由批次效應(yīng)、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì)表達(dá)影響太大等造成,可以通過(guò)繪制樣品聚類查看分組信息、關(guān)聯(lián)批次信息、處理信息和有無(wú)異常樣品(可以使用之前講過(guò)的熱圖簡(jiǎn)化，增加行或列屬性)。如果這確實(shí)是由有意義的生物變化引起的，也可以使用后面程序中的經(jīng)驗(yàn)power值。

安裝WGCNA

WGCNA依賴很多包裝，生物導(dǎo)體上的包裝需要自己安裝，而起重機(jī)上的包裝可以自動(dòng)安裝。一般WGCNA的安裝可以通過(guò)在r中運(yùn)行以下四個(gè)語(yǔ)句來(lái)完成。

建議在編譯安裝r時(shí)加入-with-blas-with-lapack，以提高矩陣運(yùn)算速度。詳見(jiàn)r和Rstudio安裝。

source(“https://bioconductor.org/biocLite.R“) biocLite(c(“AnnotationDbi”, “impute”,”GO.db”, “preprocessCore”)) site=”https://mirrors.tuna.tsinghua.edu.cn/CRAN“ install.packages(c(“WGCNA”, “stringr”, “reshape2”), repos=site)

WGCNA實(shí)戰(zhàn)

實(shí)戰(zhàn)中采用政府提供的清洗矩陣，原始矩陣信息量太大，容易誤導(dǎo)，后臺(tái)回復(fù)WGCNA獲取數(shù)據(jù)。

數(shù)據(jù)讀入

圖書(shū)館

正在加載所需的包:dynamicTreeCut正在加載所需的包:fastcluster正在附加包:“fastcluster”以下對(duì)象被“pa Package:stats”屏蔽:h CLuster = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =* *重要說(shuō)明:您的系統(tǒng)似乎支持多線程，*但在r中的WGCNA內(nèi)未啟用。*要在WGCNA內(nèi)允許所有可用內(nèi)核的多線程，請(qǐng)?jiān)趓中使用* * allowgcnathreads()* *。如有必要，請(qǐng)使用disableWGCNAThreads()禁用線程。*或者，在您的系統(tǒng)上設(shè)置以下環(huán)境變量:* ALLOW _ WGCNA _ THREADS = * *例如* * ALLOW _ WGCNA _ THREADS = 48 * *要在linux bash shell中設(shè)置環(huán)境變量，請(qǐng)?jiān)谶\(yùn)行r之前鍵入* * export ALLOW _ WGCNA _ THREES = 48 * *。其他操作系統(tǒng)或shell將*具有類似的命令來(lái)實(shí)現(xiàn)相同的目的。* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =附加包:“WGCNA”以下對(duì)象被“包:統(tǒng)計(jì)信息”屏蔽:cor

庫(kù)(reshape2)庫(kù)(stringr)

選項(xiàng)(字符串因子=假)

打開(kāi)多線程

enableWGCNAThreads()

允許并行執(zhí)行多達(dá)47個(gè)工作流程。

常規(guī)表達(dá)矩陣，log2轉(zhuǎn)換后或 Deseq2的varianceStabilizingTransformation轉(zhuǎn)換的數(shù)據(jù) 如果有批次效應(yīng)，需要事先移除，可使用removeBatchEffect 如果有系統(tǒng)偏移(可用boxplot查看基因表達(dá)分布是否一致)， 需要quantile normalization

exprMat <。- "WGCNA/LiverFemaleClean.txt "

官方推薦 “signed” 或 “signed hybrid” 為與原文檔一致，故未修改

type = "unsigned "

相關(guān)性計(jì)算 官方推薦 biweight mid-correlation & bicor corType: pearson or bicor 為與原文檔一致，故未修改

corType = "pearson "

corFnc = if else(Cortype = = " Pearson "，cor，bicor)

對(duì)二元變量，如樣本性狀信息計(jì)算相關(guān)性時(shí)， 或基因表達(dá)嚴(yán)重依賴于疾病狀態(tài)時(shí)，需設(shè)置下面參數(shù)

maxPOutliers = if else(Cortype = = " Pearson "，1，0.05)

關(guān)聯(lián)樣品性狀的二元變量時(shí)，設(shè)置

robustY = if else(Cortype = = " Pearson "，T，F(xiàn))

輸入數(shù)據(jù)

dataExpr <。- read.table(exprMat，sep='t '，row.names=1，header=T，quote= "，comment= "，check.names=F)

dim(dataExpr)

[1] 3600 134

head(dataExpr)[，1:8]

F2 _ 2 F2 _ 3 F2 _ 14 F2 _ 15 F2 _ 19 F2 _ 20 MMT 000000044-0.01810 0.0642 6.44 e-05-0.05800 0.04830-0.15197410 MMT 0000046-0.07730-0.0297 1.12 e-01-0.05890 0 0.04 430

數(shù)據(jù)篩選

篩選中位數(shù)絕對(duì)偏差在75%之前的基因，至少M(fèi)AD大于0.01，會(huì)減少計(jì)算量，丟失部分信息或者不篩選，使MAD大于0

m.mad & lt-應(yīng)用(dataExpr，1，mad)DataExprvar & lt；-m . mad = " " gt。max(分位數(shù)(m.mad，probs=seq(0，1，0.25))[2]，0.01))，]

轉(zhuǎn)換成一個(gè)矩陣，行中有樣本，列中有基因

dataExpr <。- as.data.frame(t(dataExprVar))

缺失值的檢測(cè)

gsg = goodSamplesGenes(dataExpr，verbose = 3)

標(biāo)記缺失值過(guò)多的基因和樣本…..步驟1

if(！gsg$allOK){

#或者，打印被刪除的基因和樣本名稱:if (sum(！gsg$goodGenes)>;0) printFlush(粘貼("去除基因:"，粘貼(人名(dataExpr)[！gsg$goodGenes]，collapse = "，")；if (sum(！gsg$goodSamples)>。0) printFlush(粘貼("移除樣本:"，粘貼(行名(dataExpr)[！gsg$goodSamples]，collapse = "，")；#從數(shù)據(jù)中刪除違規(guī)基因和樣本:數(shù)據(jù)表達(dá)式=數(shù)據(jù)表達(dá)式[gsg$goodSamples，gsg $ goodsgenes]

}

nGenes = ncol(DataExpr)NSAmples = nrow(DataExpr)

dim(dataExpr)

[1] 134 2697

head(dataExpr)[，1:8]

MMT 00000051 MMT 000000080 MMT 00000102 MMT 00000149 MMT 00000159 F2 _ 2-0.02260000-0.04870000 0 0.17600000-0.07680000-0.148000000 F2 _ 3 0.061700000

樣本樹(shù)= h樣本(dist(dataExpr)，method = "average ")圖(樣本樹(shù)，main = "樣本聚類以檢測(cè)異常值"，sub= "，xlab= " ")

軟閾值的篩選原則是使構(gòu)建的網(wǎng)絡(luò)更符合無(wú)標(biāo)度網(wǎng)絡(luò)特性。

powers = c(c(1:10)，seq(from = 12，to=30，by = 2))sft = PickSoftthreshold(DataExpr，powerVector=powers，networkType=type，verbose=5)

pickSoftThreshold:將使用塊大小2697。選擇軟閾值:計(jì)算給定功率的連通性…..正在研究SFT 2697人中的基因1到2697。R.sq斜率被截?cái)?。R.sq均值. k .中位數(shù). k .最大值. k . 1 0.1370 0.825 0.412 587.000 5.95 e+02 922.0 2 2 0.0416-0.332 0.630 206.000 2.02 e+02 443.0 3 3 0.2280-0.747 0.920 91.500 8.43 e+01 247。 0.855 4.780 2.01 e+00 60.1 11 12 0.8020-1.410 0.866 3.160 9.61 e-01 53.2 12 14 0.8470-1.340 0.909 2.280 4.84 e-01 47。 7 13 16 0.8850-1.250 0.932 1.750 2.64 e-01 43.1 14 18 0.8830-1.210 0.922 1.400 1.46 e-01 39.1 15 20 0 0.9110-1.180 0.926 1.150 8.35 e-02 35.6 16 22 0.916 0-1.22

par(mfrow = c(1，2)) cex1 = 0.9

橫軸是Soft threshold (power)，縱軸是無(wú)標(biāo)度網(wǎng)絡(luò)的評(píng)估參數(shù)，數(shù)值越高， 網(wǎng)絡(luò)越符合無(wú)標(biāo)度特征 (non-scale)

標(biāo)繪(SFTTFITINDICes[，3])SFTTFITINDICes[，1]，-符號(hào)(SFT $ FITINDICs[，3])SFT $ FITINDICs[，2]，標(biāo)簽=冪，cex=cex1，col="red ")

篩選標(biāo)準(zhǔn)。R-square=0.85

abline(h=0.85，col= "紅色")

Soft threshold與平均連通性

繪圖(SFtfitidices[，5]，xlab= "軟閾值(功率)"，ylab= "平均連通性"，type="n "，main = paste("平均連通性"))文本(SFtfitidices[，5]，標(biāo)簽=功率，cex=cex1，col="red ")

功率= sft $功率估計(jì)功率

[1] 6經(jīng)驗(yàn)功率(在沒(méi)有合格功率時(shí)選擇)

無(wú)向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)，沒(méi)有一個(gè)power值可以使 無(wú)標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8，平均連接度較高如在100以上，可能是由于 部分樣品與其他樣品差別太大。這可能由批次效應(yīng)、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì) 表達(dá)影響太大等造成。可以通過(guò)繪制樣品聚類查看分組信息和有無(wú)異常樣品。 如果這確實(shí)是由有意義的生物變化引起的，也可以使用下面的經(jīng)驗(yàn)power值。

if(is . na(power)){ power = if else(NSA simples & lt；20，ifelse(type == "unsigned "，9，18)，if else(NSAmples & lt；30，ifelse(type == "unsigned "，8，16)，if else(NSAmples & lt；40，ifelse(type == "unsigned "，7，14)，ifelse(type == "unsigned "，6，12)))}

網(wǎng)絡(luò)構(gòu)建一步網(wǎng)絡(luò)構(gòu)建和模塊檢測(cè)

power: 上一步計(jì)算的軟閾值 maxBlockSize: 計(jì)算機(jī)能處理的最大模塊的基因數(shù)量 (默認(rèn)5000)； 4G內(nèi)存電腦可處理8000-10000個(gè)，16G內(nèi)存電腦可以處理2萬(wàn)個(gè)，32G內(nèi)存電腦可 以處理3萬(wàn)個(gè) 計(jì)算資源允許的情況下最好放在一個(gè)block里面。 corType: pearson or bicor numericLabels: 返回?cái)?shù)字而不是顏色作為模塊的名字，后面可以再轉(zhuǎn)換為顏色 saveTOMs：最耗費(fèi)時(shí)間的計(jì)算，存儲(chǔ)起來(lái)，供后續(xù)使用 mergeCutHeight: 合并模塊的閾值，越大模塊越少

net = blockwiseModules(dataExpr，power = power，maxBlockSize = nGenes，TOMType = type，minModuleSize = 30，reassignThreshold = 0，mergeCutHeight = 0.25，numericLabels = TRUE，pamRespectsDendro = FALSE，saveTOMs=TRUE，corType = corType，maxPOutliers=maxPOutliers，loadTOMs=TRUE，saveTOMFileBase = paste0(exprMat，")。tom”)，verbose = 3)

從所有基因分塊計(jì)算模塊特征基因標(biāo)記缺失值過(guò)多的基因和樣本…..步驟1..在第一區(qū)工作。湯姆計(jì)算:鄰接....將使用47個(gè)平行螺紋。慢速計(jì)算分?jǐn)?shù):0.000000..連接....矩陣乘法....正?；?...完成。..正在將塊1的TOM保存到文件WGCNA/LiverFemleclean . txt . TOM-block . 1 . rdata…群聚..….檢測(cè)模塊..….計(jì)算模塊特征基因..….檢查模塊中的kME....從模塊1中刪除3個(gè)基因，因?yàn)樗鼈兊腒ME太低。..從模塊12中移除5個(gè)基因，因?yàn)樗鼈兊腒ME太低。..從模塊14中刪除1個(gè)基因，因?yàn)樗鼈兊腒ME太低。..合并過(guò)于接近的模塊..合并關(guān)閉模塊:合并距離小于0.25的模塊計(jì)算新的MEs…

根據(jù)模塊中基因數(shù)目的多少，降序排列，依次編號(hào)為 1-最大模塊數(shù)。 0 (grey)表示未分入任何模塊的基因。

表格(凈價(jià)$colors)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 135 472 356 333 307 303 177 158 102 94 69 66 63 62分級(jí)聚類樹(shù)顯示每個(gè)模塊的灰色基因，這些基因沒(méi)有被分類到模塊中。

Convert labels to colors for plotting

module labels = net $ colors module colors = labels 2 colors(module labels)

Plot the dendrogram and the module colors underneath 如果對(duì)結(jié)果不滿意，還可以recutBlockwiseTrees，節(jié)省計(jì)算時(shí)間

繪圖樹(shù)和顏色(網(wǎng)絡(luò)塊基因[[1]]]，“模塊顏色”，樹(shù)標(biāo)簽=假，懸掛= 0.03，添加引導(dǎo)=真，引導(dǎo)懸掛= 0.05)

畫(huà)出模塊之間的關(guān)聯(lián)圖

module eigengene, 可以繪制線圖，作為每個(gè)模塊的基因表達(dá)趨勢(shì)的展示

MEs =凈$MEs

不需要重新計(jì)算，只需更改以下名稱即可。官方教程重新計(jì)算，不用一開(kāi)始就費(fèi)心了

MeS _ col = MeS col names(MeS _ col)= paste 0(" ME "，labels 2 colors(as . numeric(str _ replace _ all(col names(MeS)，" ME "，" "))))MEs_col = orderMEs(MEs_col)

根據(jù)基因間表達(dá)量進(jìn)行聚類所得到的各模塊間的相關(guān)性圖 marDendro/marHeatmap 設(shè)置下、左、上、右的邊距

繪圖特征基因網(wǎng)絡(luò)(MEs_col，“特征基因鄰接熱圖”，marDendro = c(3，3，2，4)，marHeatmap = c(3，4，2，2)，繪圖樹(shù)圖= T，xLabelsAngle = 90)

如果有表型數(shù)據(jù)，也可以和ME數(shù)據(jù)放在一起，一起作圖

MEs_colpheno = orderMEs(cbind(MEs_col, traitData)) plotEigengeneNetworks(MEs_colpheno, “Eigengene adjacency heatmap”, marDendro = c(3,3,2,4), marHeatmap = c(3,4,2,2), plotDendrograms = T, xLabelsAngle = 90)

視覺(jué)基因網(wǎng)絡(luò)

如果采用分步計(jì)算，或設(shè)置的blocksize>=總基因數(shù)，直接load計(jì)算好的TOM結(jié)果 否則需要再計(jì)算一遍，比較耗費(fèi)時(shí)間 TOM = TOMsimilarityFromExpr(dataExpr, power=power, corType=corType, networkType=type)

加載(net$TOMFiles[1]，verbose=T)

加載對(duì)象:TOM

湯姆& lt矩陣

dissTOM = 1-TOM

Transform dissTOM with a power to make moderately strong connections more visible in the heatmap

普羅托姆= dissTOM^7

Set diagonal to NA for a nicer plot

diag(繪圖)= NA

Call the plot function 這一部分特別耗時(shí)，行列同時(shí)做層級(jí)聚類

基因圖(網(wǎng)絡(luò)熱圖，所有基因)

細(xì)胞景觀出口網(wǎng)絡(luò)

細(xì)胞角繪制的網(wǎng)絡(luò)圖見(jiàn)我們更新的視頻教程或https://bioinfo.ke.qq.com/。

探針=列名(DataExpr)dim names(TOM)& lt；-列表(探針、探針)

Export the network into edge and node list files Cytoscape can read threshold 默認(rèn)為0.5, 可以根據(jù)自己的需要調(diào)整，也可以都導(dǎo)出后在 cytoscape中再調(diào)整

cyt = exportNetworkToCytoscape(TOM，edgeFile = paste(exprMat，" . edges.txt "，sep= ")，nodeFile = paste(exprMat，" . nodes.txt "，sep= ")，weighted = TRUE，threshold = 0，nodeNames =探測(cè)器，nodeAttr = moduleColors)

相關(guān)表型數(shù)據(jù)

特質(zhì)<。- "WGCNA/TraitsClean.txt "

讀入表型數(shù)據(jù)，不是必須的

if(trait！= " "){ TraitData & lt；- read.table(file=trait，sep='t '，header=T，row.names=1，check.names=FALSE，comment= ' '，quote = ")SampleName = row names(DataExpr)TraitData = TraitData[match(SampleName，rownames(traitData))，] }

模塊與表型數(shù)據(jù)相關(guān)聯(lián)

if(CortType = = " pearsoon "){ ModTraitCor = cor(MeS _ col，traitData，use = " p ")ModTraitTP = corp value student(ModTraitCor，NSAmples)} else { ModTraitCorp = BiCRAndpvalue(MeS _ col，traitData，robustY = robustY)ModTraitCorp = ModTraitCorp }

bicor(x，y，use = use，…)中的警告:bicor:變量' y '中的零MAD。皮爾遜相關(guān)用于具有零(或缺失)MAD的單個(gè)列。

signif表示保留幾位小數(shù)

textMatrix = paste(有效(modTraitCor，2)，" n("，有效(modTraitP，1)，")，sep = ")，dim(TextMatriX)= dim(ModTraitCor)labeledHeatMap(MatriX = ModTraitCor，XLabels = col name(TraitData)，yLabels = col name(MeS _ col)，cex.lab = 0.5，ySymbols = col name(MeS _ col)，colorLabels = FALSE，colors = BlueThread ReD(50)，textMatrix = textMatrix，setStdMargins

模塊中的基因與表型數(shù)據(jù)相關(guān)聯(lián)。從上圖可以看出，MEmagenta與胰島素_ug_l關(guān)聯(lián)，選擇這兩部分進(jìn)行分析。

從上圖我們可以看出，MEmagenta和Insulin_ug_l相關(guān)模塊中的基因與表型數(shù)據(jù)相關(guān)聯(lián)

性狀跟模塊雖然求出了相關(guān)性，可以挑選最相關(guān)的那些模塊來(lái)分析， 但是模塊本身仍然包含非常多的基因，還需進(jìn)一步的尋找最重要的基因。 所有的模塊都可以跟基因算出相關(guān)系數(shù)，所有的連續(xù)型性狀也可以跟基因的表達(dá) 值算出相關(guān)系數(shù)。 如果跟性狀顯著相關(guān)基因也跟某個(gè)模塊顯著相關(guān)，那么這些基因可能就非常重要 。

計(jì)算模塊和基因之間的相關(guān)矩陣

if(CortType = = " Pearson "){ GeneModuleMembership = as . data . frame(cor(DataExpr，MEs_col，use = " p "))MMP VaLue = as . data . frame(CorpValueStudent(as . matrix(GeneModuleMembership)，NSCompleS))} else { GeneModuleMembershipa = BiCRAndVaLue(DataExpr，MEs_col，robustY = robustY)GeneModuleMembership = GeneModuleMembershipap }

計(jì)算性狀與基因的相關(guān)性矩陣

只能計(jì)算連續(xù)字符。如果是離散變量，在構(gòu)造樣本表時(shí)會(huì)轉(zhuǎn)化為0-1矩陣。

if(CortType = = " pearsoon "){ GeneraitCor = as . data . frame(cor(DataExpr，traitData，use = " p "))GeneraItp = as . data . frame(CorpValueStudent(as . matrix(GeneraitCor)，NSaCompleS))} else { GeneraitCora = BiCRAndpvalue(DataExpr，traitData，robustY = robustY)GeneraitCor = as . data . frame(GeneraitCorap)}

bicor(x，y，use = use，…)中的警告:bicor:變量' y '中的零MAD。皮爾遜相關(guān)用于具有零(或缺失)MAD的單個(gè)列。

最后把兩個(gè)相關(guān)性矩陣聯(lián)合起來(lái),指定感興趣模塊進(jìn)行分析

module = "洋紅色" pheno = "胰島素_ ug _ l " ModName = substring(col name(MeS _ col)，3)

獲取關(guān)注的列

module_column = match(module，ModName)pheno _ column = match(pheno，colnames(traitData))

獲取模塊內(nèi)的基因

模塊類=模塊類類==模塊

sizeGrWindow(7，7) par(mfrow = c(1，1))

與性狀高度相關(guān)的基因，也是與性狀相關(guān)的模型的關(guān)鍵基因

verbose散點(diǎn)圖(ABS(GenemmoduleMembership[moduleGenes，module_column])，ABS(GeneraitCor[moduleGenes，pheno_column])，xlab =粘貼(" Module Membership in "，Module，" Module ")，ylab =粘貼(" Gene significance for "，pheno)，main =粘貼(" Module Membership vs Gene signification n ")，cex.main = 1.2，cex.lab = 1.2，cex.axis = 1.2，col = module)

分步方法顯示了計(jì)算鄰接矩陣的每一步

鄰接=鄰接(dataExpr，power = power)

將鄰接矩陣轉(zhuǎn)化為拓?fù)渲丿B矩陣，降低噪聲和虛假相關(guān)性，得到距離矩陣。

湯姆=湯姆相似性(鄰接)dis湯姆= 1-湯姆

層次聚類計(jì)算基因之間的距離樹(shù)

genere = h crest(as . dist(dis STOM)，method = "average ")

模塊合并

We like large modules, so we set the minimum module size relatively high:

最小模塊大小= 30

Module identification using dynamic tree cut:

dynamic modeds = cutreeDynamic(tred = genere，distM = dissTOM，deepSplit = 2，pamRespectsDendro = FALSE，minClusterSize = minModuleSize)

Convert numeric lables into colors

dynamic colors = labels 2 colors(dynamic MODS)

通過(guò)計(jì)算模塊的代表性模型和評(píng)估模塊之間的定量相似性，組合表達(dá)式映射是相似的

模塊melist = moduleigenenes(datexpr，colors = dynamic colors)mes = melist $八個(gè)基因

Calculate dissimilarity of module eigengenes

MEDiss = 1-cor(MEs)

Cluster module eigengenes

METree = h crest(as . dist(MEDiss)，method = " average ")medis thres = 0.25

Call an automatic merging function

merge = mergeCloseModules(datExpr，dynamicColors，cutHeight = MEDissThres，verbose = 3)

The merged module colors

mergedColors = merge $ colors

Eigengenes of the new merged

逐步完成參考:

官網(wǎng)https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/術(shù)語(yǔ)解釋https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/Simulated-00-Background.pdfFAQhttps://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/faq.html生信博客 http://blog.genesino.com