WGCNA用法詳情

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析 (WGCNA, Weighted correlation network analysis)是用來描述不同樣品之間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法，可以用來鑒定高度協(xié)同變化的基因集, 并根據(jù)基因集的內(nèi)連性和基因集與表型之間的關(guān)聯(lián)鑒定候補(bǔ)生物標(biāo)記基因或治療靶點(diǎn)。

相比于只關(guān)注差異表達(dá)的基因互拾，WGCNA利用數(shù)千或近萬個(gè)變化最大的基因或全部基因的信息識(shí)別感興趣的基因集鸵钝，并與表型進(jìn)行顯著性關(guān)聯(lián)分析柜与。一是充分利用了信息妹卿，二是把數(shù)千個(gè)基因與表型的關(guān)聯(lián)轉(zhuǎn)換為數(shù)個(gè)基因集與表型的關(guān)聯(lián)小染，免去了多重假設(shè)檢驗(yàn)校正的問題。

理解WGCNA稿茉，需要先理解下面幾個(gè)術(shù)語和它們?cè)赪GCNA中的定義。

共表達(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ǎng)絡(luò)的邊屬性計(jì)算方式為abs(cor(genex, geney)) ^ power渺蒿；有向網(wǎng)絡(luò)的邊屬性計(jì)算方式為(1+cor(genex, geney)/2) ^ power; sign hybrid的邊屬性計(jì)算方式為cor(genex, geney)^power if cor>0 else 0。這種處理方式強(qiáng)化了強(qiáng)相關(guān)砰琢，弱化了弱相關(guān)或負(fù)相關(guān)蘸嘶，使得相關(guān)性數(shù)值更符合無標(biāo)度網(wǎng)絡(luò)特征，更具有生物意義陪汽。如果沒有合適的power训唱，一般是由于部分樣品與其它樣品因?yàn)槟撤N原因差別太大導(dǎo)致的，可根據(jù)具體問題移除部分樣品或查看后面的經(jīng)驗(yàn)值挚冤。

Module(模塊)：高度內(nèi)連的基因集况增。在無向網(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è)很巧妙的梳理寒矿，我們之前講過PCA分析的降維作用突琳，之前主要是拿來做可視化，現(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 overlap matrix)：把鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣例衍，以降低噪音和假相關(guān)昔期，獲得的新距離矩陣，這個(gè)信息可拿來構(gòu)建網(wǎng)絡(luò)或繪制TOM圖佛玄。

基本分析流程

構(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é)果般贼，獲得不同的基因模塊，用聚類樹的分枝和不同顏色表示。

如果有表型信息哼蛆，計(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包WGCNA是用于計(jì)算各種加權(quán)關(guān)聯(lián)分析的功能集合惕味，可用于網(wǎng)絡(luò)構(gòu)建，基因篩選玉锌，基因簇鑒定名挥，拓?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ù)來自不同的批次，需要先移除批次效應(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)。(這個(gè)根據(jù)自己的需要，看看上面寫的每個(gè)網(wǎng)絡(luò)怎么計(jì)算的煤杀，更知道怎么選擇)

無向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)眷蜈，沒有一個(gè)power值可以使無標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8或平均連接度降到100以下，可能是由于部分樣品與其他樣品差別太大造成的沈自。這可能由批次效應(yīng)酌儒、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì)表達(dá)影響太大等造成, 可以通過繪制樣品聚類查看分組信息、關(guān)聯(lián)批次信息枯途、處理信息和有無異常樣品 (可以使用之前講過的熱圖簡化忌怎，增加行或列屬性)。如果這確實(shí)是由有意義的生物變化引起的酪夷，也可以使用后面程序中的經(jīng)驗(yàn)power值榴啸。

安裝WGCNA

WGCNA依賴的包比較多，bioconductor上的包需要自己安裝晚岭，cran上依賴的包可以自動(dòng)安裝鸥印。通常在R中運(yùn)行下面4條語句就可以完成WGCNA的安裝。

建議在編譯安裝R時(shí)增加--with-blas --with-lapack提高矩陣運(yù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ù)讀入

library(WGCNA)

##?Loading?required?package:?dynamicTreeCut

##?Loading?required?package:?fastcluster

##?

##?Attaching?package:?'fastcluster'

##?The?following?object?is?masked?from?'package:stats':

##?

##?????hclust

##?==========================================================================

##?*

##?*??Package?WGCNA?1.63?loaded.

##?*

##?*????Important?note:?It?appears?that?your?system?supports?multi-threading,

##?*????but?it?is?not?enabled?within?WGCNA?in?R.?

##?*????To?allow?multi-threading?within?WGCNA?with?all?available?cores,?use?

##?*

##?*??????????allowWGCNAThreads()

##?*

##?*????within?R.?Use?disableWGCNAThreads()?to?disable?threading?if?necessary.

##?*????Alternatively,?set?the?following?environment?variable?on?your?system:

##?*

##?*??????????ALLOW_WGCNA_THREADS=<number_of_processors>

##?*

##?*????for?example?

##?*

##?*??????????ALLOW_WGCNA_THREADS=48

##?*

##?*????To?set?the?environment?variable?in?linux?bash?shell,?type?

##?*

##?*???????????export?ALLOW_WGCNA_THREADS=48

##?*

##?*?????before?running?R.?Other?operating?systems?or?shells?will

##?*?????have?a?similar?command?to?achieve?the?same?aim.

##?*

##?==========================================================================

##?

##?Attaching?package:?'WGCNA'

##?The?following?object?is?masked?from?'package:stats':

##?

##?????cor

library(reshape2)

library(stringr)

options(stringsAsFactors?=?FALSE)

#?打開多線程

enableWGCNAThreads()

##?Allowing?parallel?execution?with?up?to?47?working?processes.

#?常規(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?=?ifelse(corType=="pearson",?cor,?bicor)

#?對(duì)二元變量借跪，如樣本性狀信息計(jì)算相關(guān)性時(shí)政己，

#?或基因表達(dá)嚴(yán)重依賴于疾病狀態(tài)時(shí)，需設(shè)置下面參數(shù)

maxPOutliers?=?ifelse(corType=="pearson",1,0.05)

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

robustY?=?ifelse(corType=="pearson",T,F)

##導(dǎo)入數(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

##?MMT00000044?-0.01810??0.0642??6.44e-05?-0.05800??0.04830?-0.15197410

##?MMT00000046?-0.07730?-0.0297??1.12e-01?-0.05890??0.04430?-0.09380000

##?MMT00000051?-0.02260??0.0617?-1.29e-01??0.08710?-0.11500?-0.06502607

##?MMT00000076?-0.00924?-0.1450??2.87e-02?-0.04390??0.00425?-0.23610000

##?MMT00000080?-0.04870??0.0582?-4.83e-02?-0.03710??0.02510??0.08504274

##?MMT00000102??0.17600?-0.1890?-6.50e-02?-0.00846?-0.00574?-0.01807182

##????????????????F2_23????F2_24

##?MMT00000044?-0.00129?-0.23600

##?MMT00000046??0.09340??0.02690

##?MMT00000051??0.00249?-0.10200

##?MMT00000076?-0.06900??0.01440

##?MMT00000080??0.04450??0.00167

##?MMT00000102?-0.12500?-0.06820

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

##?篩選中位絕對(duì)偏差前75%的基因歇由，至少M(fèi)AD大于0.01

##?篩選后會(huì)降低運(yùn)算量，也會(huì)失去部分信息

##?也可不做篩選果港，使MAD大于0即可

m.mad?<-?apply(dataExpr,1,mad)

dataExprVar?<-?dataExpr[which(m.mad?>?

?????????????????max(quantile(m.mad,?probs=seq(0,?1,?0.25))[2],0.01)),]

##?轉(zhuǎn)換為樣品在行沦泌，基因在列的矩陣

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

##?檢測(cè)缺失值

gsg?=?goodSamplesGenes(dataExpr,?verbose?=?3)

##??Flagging?genes?and?samples?with?too?many?missing?values...

##???..step?1

if?(!gsg$allOK){

??#?Optionally,?print?the?gene?and?sample?names?that?were?removed:

??if?(sum(!gsg$goodGenes)>0)?

????printFlush(paste("Removing?genes:",?

?????????????????????paste(names(dataExpr)[!gsg$goodGenes],?collapse?=?",")));

??if?(sum(!gsg$goodSamples)>0)?

????printFlush(paste("Removing?samples:",?

?????????????????????paste(rownames(dataExpr)[!gsg$goodSamples],?collapse?=?",")));

??#?Remove?the?offending?genes?and?samples?from?the?data:

??dataExpr?=?dataExpr[gsg$goodSamples,?gsg$goodGenes]

}

nGenes?=?ncol(dataExpr)

nSamples?=?nrow(dataExpr)

dim(dataExpr)

##?[1]??134?2697

head(dataExpr)[,1:8]

##???????MMT00000051?MMT00000080?MMT00000102?MMT00000149?MMT00000159

##?F2_2??-0.02260000?-0.04870000??0.17600000??0.07680000?-0.14800000

##?F2_3???0.06170000??0.05820000?-0.18900000??0.18600000??0.17700000

##?F2_14?-0.12900000?-0.04830000?-0.06500000??0.21400000?-0.13200000

##?F2_15??0.08710000?-0.03710000?-0.00846000??0.12000000??0.10700000

##?F2_19?-0.11500000??0.02510000?-0.00574000??0.02100000?-0.11900000

##?F2_20?-0.06502607??0.08504274?-0.01807182??0.06222751?-0.05497686

##???????MMT00000207?MMT00000212?MMT00000241

##?F2_2???0.06870000??0.06090000?-0.01770000

##?F2_3???0.10100000??0.05570000?-0.03690000

##?F2_14??0.10900000??0.19100000?-0.15700000

##?F2_15?-0.00858000?-0.12100000??0.06290000

##?F2_19??0.10500000??0.05410000?-0.17300000

##?F2_20?-0.02441415??0.06343181??0.06627665

軟閾值篩選

##?查看是否有離群樣品

sampleTree?=?hclust(dist(dataExpr),?method?=?"average")

plot(sampleTree,?main?=?"Sample?clustering?to?detect?outliers",?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:?will?use?block?size?2697.

##??pickSoftThreshold:?calculating?connectivity?for?given?powers...

##????..working?on?genes?1?through?2697?of?2697

##????Power?SFT.R.sq??slope?truncated.R.sq?mean.k.?median.k.?max.k.

##?1??????1???0.1370??0.825??????????0.412?587.000??5.95e+02??922.0

##?2??????2???0.0416?-0.332??????????0.630?206.000??2.02e+02??443.0

##?3??????3???0.2280?-0.747??????????0.920??91.500??8.43e+01??247.0

##?4??????4???0.3910?-1.120??????????0.908??47.400??4.02e+01??154.0

##?5??????5???0.7320?-1.230??????????0.958??27.400??2.14e+01??102.0

##?6??????6???0.8810?-1.490??????????0.916??17.200??1.22e+01???83.7

##?7??????7???0.8940?-1.640??????????0.869??11.600??7.29e+00???75.4

##?8??????8???0.8620?-1.660??????????0.827???8.250??4.56e+00???69.2

##?9??????9???0.8200?-1.600??????????0.810???6.160??2.97e+00???64.2

##?10????10???0.8390?-1.560??????????0.855???4.780??2.01e+00???60.1

##?11????12???0.8020?-1.410??????????0.866???3.160??9.61e-01???53.2

##?12????14???0.8470?-1.340??????????0.909???2.280??4.84e-01???47.7

##?13????16???0.8850?-1.250??????????0.932???1.750??2.64e-01???43.1

##?14????18???0.8830?-1.210??????????0.922???1.400??1.46e-01???39.1

##?15????20???0.9110?-1.180??????????0.926???1.150??8.35e-02???35.6

##?16????22???0.9160?-1.140??????????0.927???0.968??5.02e-02???32.6

##?17????24???0.9520?-1.120??????????0.961???0.828??2.89e-02???29.9

##?18????26???0.9520?-1.120??????????0.944???0.716??1.77e-02???27.5

##?19????28???0.9380?-1.120??????????0.922???0.626??1.08e-02???25.4

##?20????30???0.9620?-1.110??????????0.951???0.551??6.49e-03???23.5

par(mfrow?=?c(1,2))

cex1?=?0.9

#?橫軸是Soft?threshold?(power)辛掠，縱軸是無標(biāo)度網(wǎng)絡(luò)的評(píng)估參數(shù)谢谦，數(shù)值越高释牺，

#?網(wǎng)絡(luò)越符合無標(biāo)度特征?(non-scale)

plot(sft$fitIndices[,1],?-sign(sft$fitIndices[,3])*sft$fitIndices[,2],

?????xlab="Soft?Threshold?(power)",

?????ylab="Scale?Free?Topology?Model?Fit,signed?R^2",type="n",

?????main?=?paste("Scale?independence"))

text(sft$fitIndices[,1],?-sign(sft$fitIndices[,3])*sft$fitIndices[,2],

?????labels=powers,cex=cex1,col="red")

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

abline(h=0.85,col="red")

#?Soft?threshold與平均連通性

plot(sft$fitIndices[,1],?sft$fitIndices[,5],

?????xlab="Soft?Threshold?(power)",ylab="Mean?Connectivity",?type="n",

?????main?=?paste("Mean?connectivity"))

text(sft$fitIndices[,1],?sft$fitIndices[,5],?labels=powers,?

?????cex=cex1,?col="red")

power?=?sft$powerEstimate

power

##?[1]?6

經(jīng)驗(yàn)power(無滿足條件的power時(shí)選用)

#?無向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)回挽，沒有一個(gè)power值可以使

#?無標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達(dá)到0.8没咙，平均連接度較高如在100以上，可能是由于

#?部分樣品與其他樣品差別太大千劈。這可能由批次效應(yīng)祭刚、樣品異質(zhì)性或?qū)嶒?yàn)條件對(duì)

#?表達(dá)影響太大等造成∏脚疲可以通過繪制樣品聚類查看分組信息和有無異常樣品涡驮。

#?如果這確實(shí)是由有意義的生物變化引起的，也可以使用下面的經(jīng)驗(yàn)power值喜滨。

if?(is.na(power)){

??power?=?ifelse(nSamples<20,?ifelse(type?==?"unsigned",?9,?18),

??????????ifelse(nSamples<30,?ifelse(type?==?"unsigned",?8,?16),

??????????ifelse(nSamples<40,?ifelse(type?==?"unsigned",?7,?14),

??????????ifelse(type?==?"unsigned",?6,?12))???????

??????????)

}

網(wǎng)絡(luò)構(gòu)建

##一步法網(wǎng)絡(luò)構(gòu)建：One-step?network?construction?and?module?detection##

#?power:?上一步計(jì)算的軟閾值

#?maxBlockSize:?計(jì)算機(jī)能處理的最大模塊的基因數(shù)量?(默認(rèn)5000)捉捅；

#??4G內(nèi)存電腦可處理8000-10000個(gè)，16G內(nèi)存電腦可以處理2萬個(gè)虽风，32G內(nèi)存電腦可

#??以處理3萬個(gè)

#??計(jì)算資源允許的情況下最好放在一個(gè)block里面锯梁。

#?corType:?pearson?or?bicor

#?numericLabels:?返回?cái)?shù)字而不是顏色作為模塊的名字，后面可以再轉(zhuǎn)換為顏色

#?saveTOMs：最耗費(fèi)時(shí)間的計(jì)算焰情，存儲(chǔ)起來，供后續(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)

##??Calculating?module?eigengenes?block-wise?from?all?genes

##????Flagging?genes?and?samples?with?too?many?missing?values...

##?????..step?1

##??..Working?on?block?1?.

##?????TOM?calculation:?adjacency..

##?????..will?use?47?parallel?threads.

##??????Fraction?of?slow?calculations:?0.000000

##?????..connectivity..

##?????..matrix?multiplication?(system?BLAS)..

##?????..normalization..

##?????..done.

##????..saving?TOM?for?block?1?into?file?WGCNA/LiverFemaleClean.txt.tom-block.1.RData

##??....clustering..

##??....detecting?modules..

##??....calculating?module?eigengenes..

##??....checking?kME?in?modules..

##??????..removing?3?genes?from?module?1?because?their?KME?is?too?low.

##??????..removing?5?genes?from?module?12?because?their?KME?is?too?low.

##??????..removing?1?genes?from?module?14?because?their?KME?is?too?low.

##??..merging?modules?that?are?too?close..

##??????mergeCloseModules:?Merging?modules?whose?distance?is?less?than?0.25

##????????Calculating?new?MEs...

#?根據(jù)模塊中基因數(shù)目的多少内舟，降序排列，依次編號(hào)為?`1-最大模塊數(shù)`初橘。

#?**0?(grey)**表示**未**分入任何模塊的基因验游。?

table(net$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í)聚類樹展示各個(gè)模塊

##?灰色的為**未分類**到模塊的基因。

#?Convert?labels?to?colors?for?plotting

moduleLabels?=?net$colors

moduleColors?=?labels2colors(moduleLabels)

#?Plot?the?dendrogram?and?the?module?colors?underneath

#?如果對(duì)結(jié)果不滿意保檐，還可以recutBlockwiseTrees耕蝉，節(jié)省計(jì)算時(shí)間

plotDendroAndColors(net$dendrograms[[1]],?moduleColors[net$blockGenes[[1]]],

????????????????????"Module?colors",

????????????????????dendroLabels?=?FALSE,?hang?=?0.03,

????????????????????addGuide?=?TRUE,?guideHang?=?0.05)

繪制模塊之間相關(guān)性圖

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

MEs?=?net$MEs

###?不需要重新計(jì)算夜只，改下列名字就好

###?官方教程是重新計(jì)算的垒在，起始可以不用這么麻煩

MEs_col?=?MEs

colnames(MEs_col)?=?paste0("ME",?labels2colors(

??as.numeric(str_replace_all(colnames(MEs),"ME",""))))

MEs_col?=?orderMEs(MEs_col)

#?根據(jù)基因間表達(dá)量進(jìn)行聚類所得到的各模塊間的相關(guān)性圖

#?marDendro/marHeatmap?設(shè)置下、左扔亥、上场躯、右的邊距

plotEigengeneNetworks(MEs_col,?"Eigengene?adjacency?heatmap",?

??????????????????????marDendro?=?c(3,3,2,4),

??????????????????????marHeatmap?=?c(3,4,2,2),?plotDendrograms?=?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ò) (TOM plot)

#?如果采用分步計(jì)算踢关，或設(shè)置的blocksize>=總基因數(shù)，直接load計(jì)算好的TOM結(jié)果

#?否則需要再計(jì)算一遍粘茄，比較耗費(fèi)時(shí)間

#?TOM?=?TOMsimilarityFromExpr(dataExpr,?power=power,?corType=corType,?networkType=type)

load(net$TOMFiles[1],?verbose=T)

##?Loading?objects:

##???TOM

TOM?<-?as.matrix(TOM)

dissTOM?=?1-TOM

#?Transform?dissTOM?with?a?power?to?make?moderately?strong?

#?connections?more?visible?in?the?heatmap

plotTOM?=?dissTOM^7

#?Set?diagonal?to?NA?for?a?nicer?plot

diag(plotTOM)?=?NA

#?Call?the?plot?function

#?這一部分特別耗時(shí)签舞，行列同時(shí)做層級(jí)聚類

TOMplot(plotTOM,?net$dendrograms,?moduleColors,?

????????main?=?"Network?heatmap?plot,?all?genes")

導(dǎo)出網(wǎng)絡(luò)用于Cytoscape

Cytoscape繪制網(wǎng)絡(luò)圖見我們更新版的視頻教程或https://bioinfo.ke.qq.com/秕脓。

probes?=?colnames(dataExpr)

dimnames(TOM)?<-?list(probes,?probes)

#?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?=?probes,?nodeAttr?=?moduleColors)

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

trait?<-?"WGCNA/TraitsClean.txt"

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

if(trait?!=?"")?{

??traitData?<-?read.table(file=trait,?sep='\t',?header=T,?row.names=1,

??????????????????????????check.names=FALSE,?comment='',quote="")

??sampleName?=?rownames(dataExpr)

??traitData?=?traitData[match(sampleName,?rownames(traitData)),?]

}

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

if?(corType=="pearsoon")?{

??modTraitCor?=?cor(MEs_col,?traitData,?use?=?"p")

??modTraitP?=?corPvalueStudent(modTraitCor,?nSamples)

}?else?{

??modTraitCorP?=?bicorAndPvalue(MEs_col,?traitData,?robustY=robustY)

??modTraitCor?=?modTraitCorP$bicor

??modTraitP???=?modTraitCorP$p

}

##?Warning?in?bicor(x,?y,?use?=?use,?...):?bicor:?zero?MAD?in?variable?'y'.

##?Pearson?correlation?was?used?for?individual?columns?with?zero?(or?missing)

##?MAD.

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

textMatrix?=?paste(signif(modTraitCor,?2),?"\n(",?signif(modTraitP,?1),?")",?sep?=?"")

dim(textMatrix)?=?dim(modTraitCor)

labeledHeatmap(Matrix?=?modTraitCor,?xLabels?=?colnames(traitData),?

???????????????yLabels?=?colnames(MEs_col),?

???????????????cex.lab?=?0.5,?

???????????????ySymbols?=?colnames(MEs_col),?colorLabels?=?FALSE,?

???????????????colors?=?blueWhiteRed(50),?

???????????????textMatrix?=?textMatrix,?setStdMargins?=?FALSE,?

???????????????cex.text?=?0.5,?zlim?=?c(-1,1),

???????????????main?=?paste("Module-trait?relationships"))

模塊內(nèi)基因與表型數(shù)據(jù)關(guān)聯(lián), 從上圖可以看到MEmagenta與Insulin_ug_l相關(guān)吠架，選取這兩部分進(jìn)行分析。

##?從上圖可以看到MEmagenta與Insulin_ug_l相關(guān)

##?模塊內(nèi)基因與表型數(shù)據(jù)關(guān)聯(lián)

#?性狀跟模塊雖然求出了相關(guān)性师妙，可以挑選最相關(guān)的那些模塊來分析诵肛，

#?但是模塊本身仍然包含非常多的基因，還需進(jìn)一步的尋找最重要的基因默穴。

#?所有的模塊都可以跟基因算出相關(guān)系數(shù)怔檩，所有的連續(xù)型性狀也可以跟基因的表達(dá)

#?值算出相關(guān)系數(shù)。

#?如果跟性狀顯著相關(guān)基因也跟某個(gè)模塊顯著相關(guān)蓄诽，那么這些基因可能就非常重要

#?薛训。

###?計(jì)算模塊與基因的相關(guān)性矩陣

if?(corType=="pearsoon")?{

??geneModuleMembership?=?as.data.frame(cor(dataExpr,?MEs_col,?use?=?"p"))

??MMPvalue?=?as.data.frame(corPvalueStudent(

?????????????as.matrix(geneModuleMembership),?nSamples))

}?else?{

??geneModuleMembershipA?=?bicorAndPvalue(dataExpr,?MEs_col,?robustY=robustY)

??geneModuleMembership?=?geneModuleMembershipA$bicor

??MMPvalue???=?geneModuleMembershipA$p

}

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

##?只有連續(xù)型性狀才能進(jìn)行計(jì)算，如果是離散變量仑氛，在構(gòu)建樣品表時(shí)就轉(zhuǎn)為0-1矩陣乙埃。

if?(corType=="pearsoon")?{

??geneTraitCor?=?as.data.frame(cor(dataExpr,?traitData,?use?=?"p"))

??geneTraitP?=?as.data.frame(corPvalueStudent(

?????????????as.matrix(geneTraitCor),?nSamples))

}?else?{

??geneTraitCorA?=?bicorAndPvalue(dataExpr,?traitData,?robustY=robustY)

??geneTraitCor?=?as.data.frame(geneTraitCorA$bicor)

??geneTraitP???=?as.data.frame(geneTraitCorA$p)

}

##?Warning?in?bicor(x,?y,?use?=?use,?...):?bicor:?zero?MAD?in?variable?'y'.

##?Pearson?correlation?was?used?for?individual?columns?with?zero?(or?missing)

##?MAD.

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

module?=?"magenta"

pheno?=?"Insulin_ug_l"

modNames?=?substring(colnames(MEs_col),?3)

#?獲取關(guān)注的列

module_column?=?match(module,?modNames)

pheno_column?=?match(pheno,colnames(traitData))

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

moduleGenes?=?moduleColors?==?module

sizeGrWindow(7,?7)

par(mfrow?=?c(1,1))

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

verboseScatterplot(abs(geneModuleMembership[moduleGenes,?module_column]),

???????????????????abs(geneTraitCor[moduleGenes,?pheno_column]),

???????????????????xlab?=?paste("Module?Membership?in",?module,?"module"),

???????????????????ylab?=?paste("Gene?significance?for",?pheno),

???????????????????main?=?paste("Module?membership?vs.?gene?significance\n"),

???????????????????cex.main?=?1.2,?cex.lab?=?1.2,?cex.axis?=?1.2,?col?=?module)

分步法展示每一步都做了什么

###?計(jì)算鄰接矩陣

adjacency?=?adjacency(dataExpr,?power?=?power)

###?把鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣锯岖，以降低噪音和假相關(guān)介袜，獲得距離矩陣。

TOM?=?TOMsimilarity(adjacency)

dissTOM?=?1-TOM

###?層級(jí)聚類計(jì)算基因之間的距離樹?

geneTree?=?hclust(as.dist(dissTOM),?method?=?"average")

###?模塊合并

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

minModuleSize?=?30

#?Module?identification?using?dynamic?tree?cut:

dynamicMods?=?cutreeDynamic(dendro?=?geneTree,?distM?=?dissTOM,

????????????????????????????deepSplit?=?2,?pamRespectsDendro?=?FALSE,

????????????????????????????minClusterSize?=?minModuleSize)

#?Convert?numeric?lables?into?colors

dynamicColors?=?labels2colors(dynamicMods)

###?通過計(jì)算模塊的代表性模式和模塊之間的定量相似性評(píng)估出吹，合并表達(dá)圖譜相似

的模塊

MEList?=?moduleEigengenes(datExpr,?colors?=?dynamicColors)

MEs?=?MEList$eigengenes

#?Calculate?dissimilarity?of?module?eigengenes

MEDiss?=?1-cor(MEs)

#?Cluster?module?eigengenes

METree?=?hclust(as.dist(MEDiss),?method?=?"average")

MEDissThres?=?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

##?分步法完結(jié)

Reference:

官網(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.pdf

FAQhttps://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/faq.html

生信博客http://blog.genesino.com

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末遇伞，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子捶牢，更是在濱河造成了極大的恐慌鸠珠，老刑警劉巖，帶你破解...
沈念sama閱讀 206,126評(píng)論 6贊 481
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件秋麸，死亡現(xiàn)場(chǎng)離奇詭異渐排，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)灸蟆，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,254評(píng)論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門驯耻，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人炒考，你說我怎么就攤上這事吓歇。” “怎么了票腰？”我有些...
開封第一講書人閱讀 152,445評(píng)論 0贊 341
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵城看，是天一觀的道長。經(jīng)常有香客問我杏慰，道長测柠，這世上最難降的妖魔是什么炼鞠？我笑而不...
開封第一講書人閱讀 55,185評(píng)論 1贊 278
?港島之戀（遺憾婚禮）
正文為了忘掉前任，我火速辦了婚禮轰胁，結(jié)果婚禮上谒主，老公的妹妹穿的比我還像新娘。我一直安慰自己赃阀，他們只是感情好霎肯，可當(dāng)我...
茶點(diǎn)故事閱讀 64,178評(píng)論 5贊 371
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布。她就那樣靜靜地躺著榛斯，像睡著了一般观游。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上驮俗，一...
開封第一講書人閱讀 48,970評(píng)論 1贊 284
城市分裂傳說
那天懂缕，我揣著相機(jī)與錄音，去河邊找鬼王凑。笑死搪柑，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的索烹。我是一名探鬼主播工碾，決...
沈念sama閱讀 38,276評(píng)論 3贊 399
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場(chǎng)噩夢(mèng)啊……” “哼百姓！你這毒婦竟也來了倚喂？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 36,927評(píng)論 0贊 259
萬榮殺人案實(shí)錄
序言：老撾萬榮一對(duì)情侶失蹤瓣戚，失蹤者是張志新（化名）和其女友劉穎，沒想到半個(gè)月后焦读，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體子库，經(jīng)...
沈念sama閱讀 43,400評(píng)論 1贊 300
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 35,883評(píng)論 2贊 323
?白月光啟示錄
正文我和宋清朗相戀三年矗晃，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了仑嗅。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片拷恨。...
茶點(diǎn)故事閱讀 37,997評(píng)論 1贊 333
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡痢掠，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出祷嘶，到底是詐尸還是另有隱情俗他，我是刑警寧澤脖捻，帶...
沈念sama閱讀 33,646評(píng)論 4贊 322
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站兆衅，受9級(jí)特大地震影響地沮，放射性物質(zhì)發(fā)生泄漏嗜浮。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 39,213評(píng)論 3贊 307
男人毒藥：我在死后第九天來索命
文/蒙蒙一摩疑、第九天我趴在偏房一處隱蔽的房頂上張望危融。院中可真熱鬧，春花似錦雷袋、人聲如沸吉殃。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,204評(píng)論 0贊 19
一樁弒父案楷怒，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽蛋勺。三九已至，卻和暖如春率寡，著一層夾襖步出監(jiān)牢的瞬間迫卢，已是汗流浹背。一陣腳步聲響...
開封第一講書人閱讀 31,423評(píng)論 1贊 260
情欲美人皮
我被黑心中介騙來泰國打工冶共，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留乾蛤，地道東北人。一個(gè)月前我還...
沈念sama閱讀 45,423評(píng)論 2贊 352
代替公主和親
正文我出身青樓捅僵，卻偏偏與公主長得像家卖，于是被迫代替她去往敵國和親。傳聞我的和親對(duì)象是個(gè)殘疾皇子庙楚，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 42,722評(píng)論 2贊 345

WGCNA用法詳情

推薦閱讀更多精彩內(nèi)容