WGCNA基本概念

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

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

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

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

基本分析過程

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

WGCNA包實戰(zhàn)

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

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

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

安裝WGCNA

WGCNA依賴很多包裝，生物導(dǎo)體上的包裝需要自己安裝，而起重機上的包裝可以自動安裝。一般WGCNA的安裝可以通過在r中運行以下四個語句來完成。

建議在編譯安裝r時加入-with-blas-with-lapack，以提高矩陣運算速度。詳見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實戰(zhàn)

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

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

圖書館

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

庫(reshape2)庫(stringr)

選項(字符串因子=假)

打開多線程

enableWGCNAThreads()

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

常規(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)性計算 官方推薦 biweight mid-correlation & bicor corType: pearson or bicor 為與原文檔一致，故未修改

corType = "pearson "

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

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

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

關(guān)聯(lián)樣品性狀的二元變量時，設(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ù)絕對偏差在75%之前的基因，至少MAD大于0.01，會減少計算量，丟失部分信息或者不篩選，使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)換成一個矩陣，行中有樣本，列中有基因

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

缺失值的檢測

gsg = goodSamplesGenes(dataExpr，verbose = 3)

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

if(！gsg$allOK){

#或者，打印被刪除的基因和樣本名稱:if (sum(！gsg$goodGenes)>;0)print flush(paste(" remaining genes:"，paste(names(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

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

軟閾值的篩選原則是使構(gòu)建的網(wǎng)絡(luò)更符合無標(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。選擇軟閾值:計算給定功率的連通性…..正在研究SFT 2697人中的基因1到2697。R.sq斜率被截斷。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 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)，縱軸是無標(biāo)度網(wǎng)絡(luò)的評估參數(shù)，數(shù)值越高， 網(wǎng)絡(luò)越符合無標(biāo)度特征 (non-scale)

標(biāo)繪(SFTTFITINDICes[，3])SFTTFITINDICes[，1]，-符號(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 $功率估計功率

[1] 6經(jīng)驗功率(在沒有合格功率時選擇)

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

power: 上一步計算的軟閾值 maxBlockSize: 計算機能處理的最大模塊的基因數(shù)量 (默認(rèn)5000)； 4G內(nèi)存電腦可處理8000-10000個，16G內(nèi)存電腦可以處理2萬個，32G內(nèi)存電腦可 以處理3萬個 計算資源允許的情況下最好放在一個block里面。 corType: pearson or bicor numericLabels: 返回數(shù)字而不是顏色作為模塊的名字，后面可以再轉(zhuǎn)換為顏色 saveTOMs：最耗費時間的計算，存儲起來，供后續(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)

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

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

表格(凈價$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分級聚類樹顯示每個模塊的灰色基因，這些基因沒有被分類到模塊中。

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 如果對結(jié)果不滿意，還可以recutBlockwiseTrees，節(jié)省計算時間

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

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

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

MEs =凈$MEs

不需要重新計算，只需更改以下名稱即可。官方教程重新計算，不用一開始就費心了

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)，繪圖樹圖= 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)

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

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

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

加載對象: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 這一部分特別耗時，行列同時做層級聚類

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

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

細(xì)胞角繪制的網(wǎng)絡(luò)圖見我們更新的視頻教程或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 =探測器，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的單個列。

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

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

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

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

只能計算連續(xù)字符。如果是離散變量，在構(gòu)造樣本表時會轉(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的單個列。

最后把兩個相關(guān)性矩陣聯(lián)合起來,指定感興趣模塊進(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散點圖(ABS(GenemoduleMembership[ModuleGenes，module_column])，ABS(GeneraitCor[ModuleGenes，pheno_column])，xlab =粘貼(" Module Membership in "，Module，" Module ")，ylab =粘貼(" Gene measurement for "，pheno)，main =粘貼(" Module Membership vs Gene measurement n ")，cex.main = 1.2，cex.lab = 1.2，cex.axis = 1.2，col = module)

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

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

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

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

層次聚類計算基因之間的距離樹

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)

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

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

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ù)語解釋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