原文地址：https://mp.weixin.qq.com/s/3K7vQVOHZkUvYGI8XYUwiA

? ? ? ? ? ? ? ? ? ? 今天推薦給大家一個(gè)R包WGCNA需五，針對(duì)我們的表達(dá)譜數(shù)據(jù)進(jìn)行分析剂习。

簡(jiǎn)單介紹：WGCNA首先假定基因網(wǎng)絡(luò)服從無(wú)尺度分布氧映，并定義基因共表達(dá)相關(guān)矩陣骨田、基因網(wǎng)絡(luò)形成的鄰接函數(shù)瞳抓，然后計(jì)算不同節(jié)點(diǎn)的相異系數(shù)竣贪，并據(jù)此構(gòu)建系統(tǒng)聚類(lèi)樹(shù)研叫。該聚類(lèi)樹(shù)的不同分支代表不同的基因模塊(module)，模塊內(nèi)基因共表達(dá)程度高关面，而分屬不同模塊的基因共表達(dá)程度低坦袍。

主要應(yīng)用：如果某些基因在一個(gè)生理過(guò)程或不同組織中總是具有相類(lèi)似的表達(dá)變化，那么我們有理由認(rèn)為這些基因在功能上是相關(guān)的等太，可以把它們定義為一個(gè)模塊(module)捂齐。當(dāng)基因module被定義出來(lái)后，我們可以利用這些結(jié)果做很多進(jìn)一步的工作缩抡，如篩選module的核心基因奠宜，關(guān)聯(lián)性狀，代謝通路建模瞻想，建立基因互作網(wǎng)絡(luò)等压真。

好了，言歸正傳蘑险，我們開(kāi)始一步步進(jìn)行演示滴肿！

載入WGCNA包，設(shè)置隨機(jī)種子漠其，默認(rèn)數(shù)據(jù)不進(jìn)行因子轉(zhuǎn)換

先把原始數(shù)據(jù)轉(zhuǎn)列嘴高，轉(zhuǎn)成橫排是探針（基因）竿音，縱排是個(gè)體的順序和屎，先變成數(shù)列，用as.data.fame春瞬，然后改列名rownames(design) <- design[,1]

design <- design[,-1]

##datExpr<-as.data.frame(datExpr)? (有可能需要先把數(shù)值轉(zhuǎn)為數(shù)據(jù)集）

##> datExpr1<-read.table("test.txt",header=T,stringsAsFactors = F)

##> datExpr1<-t(datExpr1)

##> colnames(datExpr1)<-datExpr1[1,]

##> datExpr1<-datExpr1[-1,]

##> datExpr1<-as.data.frame(datExpr1)

library(WGCNA)

set.seed(1)

options(stringsAsFactors = F)

構(gòu)造性狀數(shù)據(jù)（亦或是分組數(shù)據(jù)）

obs<-paste("sam",1:10,sep="")

sample<-as.data.frame(diag(x=1,nrow = length(obs)))

rownames(sample)<-obs

colnames(sample)<-1:10

構(gòu)造表達(dá)量數(shù)據(jù)

datExpr<-as.data.frame(t(matrix(runif(30000)+5,3000,10)))

rownames(datExpr)<-obs

names(datExpr)<-paste("transcript",1:3000,sep="")

明確樣本數(shù)和基因數(shù)

nGenes = ncol(datExpr)

nSamples = nrow(datExpr)

首先針對(duì)樣本做個(gè)系統(tǒng)聚類(lèi)樹(shù)

datExpr_tree<-hclust(dist(datExpr), method = "average")

par(mar = c(0,5,2,0))

plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2,

cex.axis = 1, cex.main = 1,cex.lab=1)

針對(duì)前面構(gòu)造的樣品矩陣添加對(duì)應(yīng)顏色

sample_colors <- numbers2colors(sample, colors = c("white","blue"),signed = FALSE)

構(gòu)造10個(gè)樣品的系統(tǒng)聚類(lèi)樹(shù)及性狀熱圖

par(mar = c(1,4,3,1),cex=0.8)

plotDendroAndColors(datExpr_tree, sample_colors,

groupLabels = colnames(sample),

cex.dendroLabels = 0.8,

marAll = c(1, 4, 3, 1),

cex.rowText = 0.01,

main = "Sample dendrogram and trait heatmap")

這個(gè)有什么意義呢柴信？

你可以將樣本分為正常組和對(duì)照組，或者野生型和突變型等宽气，從而可以查看樣本聚類(lèi)情況随常！

針對(duì)10個(gè)樣品繪制主成分圖（在這里不考慮分組情況）

pca = prcomp(datExpr)

sampletype<-rownames(sample)

par(mar = c(4,4,4,6))

plot(pca$x[,c(1,2)],pch=16,col=rep(rainbow(nSamples),each=1),cex=1.5,main = "PCA map")

text(pca$x[,c(1,2)],row.names(pca$x),col="black",pos=3,cex=1)

legend("right",legend=sampletype,ncol = 1,xpd=T,inset = -0.15,

pch=16,cex=1,col=rainbow(length(sampletype)),bty="n")

library(scatterplot3d)

par(mar = c(4,4,4,4))

scatterplot3d(pca$x[,1:3], highlight.3d=F, col.axis="black",color = rep(rainbow(nSamples),each=1),cex.symbols=1.5,cex.lab=1,cex.axis=1, col.grid="lightblue", main="PCA map", pch=16)

legend("topleft",legend = row.names(pca$x) ,pch=16,cex=1,col=rainbow(nSamples), ncol = 2,bty="n")

選擇合適“軟閥值（soft thresholding power）”beta

powers = c(1:30)

pow<-pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

設(shè)置網(wǎng)絡(luò)構(gòu)建參數(shù)選擇范圍潜沦，計(jì)算無(wú)尺度分布拓?fù)渚仃?/p>

par(mfrow = c(1,2))

plot(pow$fitIndices[,1], -sign(pow$fitIndices[,3])*pow$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))

text(pow$fitIndices[,1], -sign(pow$fitIndices[,3])*pow$fitIndices[,2],labels=powers,cex=0.5,col="red")

plot(pow$fitIndices[,1], pow$fitIndices[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))

text(pow$fitIndices[,1], pow$fitIndices[,5], labels=powers, cex=0.6,col="red")

參數(shù)beta取值默認(rèn)是1：30，上述圖形的橫軸均代表權(quán)重參數(shù)β绪氛，左圖縱軸代表對(duì)應(yīng)的網(wǎng)絡(luò)中l(wèi)og(k)與log(p(k))相關(guān)系數(shù)的平方唆鸡。相關(guān)系數(shù)的平方越高，說(shuō)明該網(wǎng)絡(luò)越逼近無(wú)網(wǎng)路尺度的分布枣察。右圖的縱軸代表對(duì)應(yīng)的基因模塊中所有基因鄰接函數(shù)的均值争占。

在這里，我們選擇β=6構(gòu)建基因網(wǎng)絡(luò)序目。

接下來(lái)是非常重要的一塊內(nèi)容就是構(gòu)架基因網(wǎng)絡(luò)

大體思路：計(jì)算基因間的鄰接性臂痕，根據(jù)鄰接性計(jì)算基因間的相似性，然后推出基因間的相異性系數(shù)猿涨，并據(jù)此得到基因間的系統(tǒng)聚類(lèi)樹(shù)握童。然后按照混合動(dòng)態(tài)剪切樹(shù)的標(biāo)準(zhǔn)，設(shè)置每個(gè)基因模塊最少的基因數(shù)目為30叛赚。根據(jù)動(dòng)態(tài)剪切法確定基因模塊后澡绩，再次分析，依次計(jì)算每個(gè)模塊的特征向量值俺附，然后對(duì)模塊進(jìn)行聚類(lèi)分析英古，將距離較近的模塊合并為新的模塊。

1昙读、計(jì)算樹(shù)之間的鄰接性

adjacency = adjacency(datExpr, power = softPower)?

2召调、計(jì)算樹(shù)之間的相異性

TOM = TOMsimilarity(adjacency)

dissTOM = 1-TOM

3、聚類(lèi)分析

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

4蛮浑、設(shè)置基因模塊中最少基因包含30個(gè)

minModuleSize = 30

dynamicMods

= cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit =

2,pamRespectsDendro = FALSE, minClusterSize = minModuleSize)

5唠叛、基因分組上色

dynamicColors = labels2colors(dynamicMods)

table(dynamicColors)

6、計(jì)算基因模塊的特征值

MEList = moduleEigengenes(datExpr, colors = dynamicColors)

MEs = MEList$eigengenes

MEDiss = 1-cor(MEs)

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

7沮稚、建立系統(tǒng)聚類(lèi)樹(shù)

MEDissThres = 0.4

plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")

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

7艺沼、基因模塊合并

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

mergedColors = merge$colors

mergedMEs = merge$newMEs

8、繪制基因網(wǎng)絡(luò)圖

plotDendroAndColors(geneTree,

cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged

dynamic"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang

= 0.05)

拓展分析

不同模塊基因熱圖及關(guān)鍵基因的表達(dá)

person=cor(datExpr,use = 'p')

corr<-TOM

Colors<-mergedColors

colnames(corr)<-colnames(datExpr)

rownames(corr)<-colnames(datExpr)

names(Colors)<-colnames(datExpr)

colnames(person)<-colnames(datExpr)

rownames(person)<-colnames(datExpr)

umc = unique(mergedColors)

lumc = length(umc)

for (i in c(1:lumc)){

if(umc[i]== "grey"){

next

}

ME=MEs[, paste("ME",umc[i], sep="")]

par(mfrow=c(2,1), mar=c(0.3, 5.5, 3, 2))

plotMat(t(scale(datExpr[,Colors==umc[i]])),nrgcols=30,rlabels=F,rcols=umc[i], main=umc[i], cex.main=2)

par(mar=c(5, 4.2, 0, 0.7))

barplot(ME, col=umc[i], main="", cex.main=2,ylab="eigengene expression",xlab="array sample")

}

一共會(huì)生成36個(gè)基因模塊熱圖蕴掏，由于篇幅有限障般，僅僅展示2個(gè)

基因共表達(dá)網(wǎng)絡(luò)熱圖

kME=signedKME(datExpr, mergedMEs, outputColumnName = "kME", corFnc = "cor", corOptions = "use = 'p'")

if (dim(datExpr)[2]>=1500) nSelect=1500 else nSelect=dim(datExpr)[2]

set.seed(1)

select = sample(nGenes, size = nSelect)

selectTOM = dissTOM[select, select]

selectTree = hclust(as.dist(selectTOM), method = "average")

selectColors = moduleColors[select]

plotDiss = selectTOM^7

TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot")

模塊相關(guān)性熱圖

MEs = moduleEigengenes(datExpr, Colors)$eigengenes

MET = orderMEs(MEs)

plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)

模塊與性狀相關(guān)性熱圖

moduleTraitCor = cor(MET, sample, use = "p")

moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "")

dim(textMatrix) = dim(moduleTraitCor)

par(mar = c(2, 4, 2, 1.5))

labeledHeatmap(Matrix = moduleTraitCor,

xLabelsAngle = 0,

cex.lab = 0.5,

xLabels = colnames(sample),

yLabels = names(MET),

ySymbols = names(MET),

colorLabels = FALSE,

colors = blueWhiteRed(100),

textMatrix = textMatrix,

setStdMargins = FALSE,

cex.text = 0.5,

zlim = c(-1,1),

yColorWidth=0.02,

xColorWidth = 0.05,

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

基因的系統(tǒng)樹(shù)圖及性狀相關(guān)性熱圖

geneTraitSignificance = as.data.frame(cor(datExpr, sample, use = "p"))

geneTraitColor=as.data.frame(numbers2colors(geneTraitSignificance,signed=TRUE,colors = colorRampPalette(c("blue","white","red"))(100)))

names(geneTraitColor)= colnames(sample)

par(mar = c(3.5, 7, 2, 1))

plotDendroAndColors(geneTree, cbind(mergedColors, geneTraitColor),

c("Module",colnames(sample)),dendroLabels = FALSE, hang = 0.03,

addGuide = TRUE, guideHang = 0.05)

不同模塊的基因顯著性圖

geneTraitSignificance = as.data.frame(cor(datExpr, sample, use = "p"))

GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))

names(geneTraitSignificance) = paste("GS.", colnames(sample), sep="")

names(GSPvalue) = paste("GS.", colnames(sample), sep="")

modNames = substring(names(MET), 3)

for (module in modNames){

if(module== "grey"){ next }

column = match(module, modNames); # col number of interesting modules

moduleGenes = Colors==module;

par(mfrow = c(1,1))

verboseScatterplot(abs(kME[moduleGenes, column]),

abs(geneTraitSignificance[moduleGenes, 1]),

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

ylab = "Gene significance",

main = paste("Module membership vs. gene significance

"),cex.main = 1.2, cex.lab = 1.2, pch=19,cex.axis = 1.2, col = module)}

好了，今天就到這里……

相關(guān)函數(shù)? plotMEpairs() （這個(gè)也很重要）

一步法：

TOM = TOMsimilarityFromExpr(data_matrix_mv, power = 6);

# Read in the annotation file

# annot = read.csv(file = "GeneAnnotation.csv");

# Select modules需要修改盛杰，選擇需要導(dǎo)出的模塊顏色

modules = c("turquoise");

# Select module probes選擇模塊探測(cè)

probes = colnames(data_matrix_mv)

inModule = is.finite(match(mergedColors, modules));

modProbes = probes[inModule];

#modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];

# Select the corresponding Topological Overlap

modTOM = TOM[inModule, inModule];

dimnames(modTOM) = list(modProbes, modProbes)

# Export the network into edge and node list files Cytoscape can read

cyt = exportNetworkToCytoscape(modTOM,

???????????????????????????????edgeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),

???????????????????????????????nodeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),

???????????????????????????????weighted = TRUE,

???????????????????????????????threshold = 0.02,

???????????????????????????????nodeNames = modProbes,

???????????????????????????????#altNodeNames = modGenes,

???????????????????????????????nodeAttr = mergedColors[inModule]);

作者：蘇慕晨楓

鏈接：http://www.reibang.com/p/72c5b4e9ac3e

來(lái)源：簡(jiǎn)書(shū)

簡(jiǎn)書(shū)著作權(quán)歸作者所有挽荡，任何形式的轉(zhuǎn)載都請(qǐng)聯(lián)系作者獲得授權(quán)并注明出處。