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

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

簡單介紹：WGCNA首先假定基因網(wǎng)絡(luò)服從無尺度分布话瞧，并定義基因共表達(dá)相關(guān)矩陣柜蜈、基因網(wǎng)絡(luò)形成的鄰接函數(shù)啃憎，然后計算不同節(jié)點的相異系數(shù)，并據(jù)此構(gòu)建系統(tǒng)聚類樹秕脓。該聚類樹的不同分支代表不同的基因模塊(module)戚揭，模塊內(nèi)基因共表達(dá)程度高诱告，而分屬不同模塊的基因共表達(dá)程度低。

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

好了，言歸正傳维雇，我們開始一步步進(jìn)行演示淤刃！

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

先把原始數(shù)據(jù)轉(zhuǎn)列钝凶，轉(zhuǎn)成橫排是探針（基因）仪芒，縱排是個體的順序唁影，先變成數(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)

首先針對樣本做個系統(tǒng)聚類樹

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)

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

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

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

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")

這個有什么意義呢据沈？

你可以將樣本分為正常組和對照組，或者野生型和突變型等饺蔑，從而可以查看樣本聚類情況锌介！

針對10個樣品繪制主成分圖（在這里不考慮分組情況）

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ù)選擇范圍，計算無尺度分布拓?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ù)β孔祸，左圖縱軸代表對應(yīng)的網(wǎng)絡(luò)中l(wèi)og(k)與log(p(k))相關(guān)系數(shù)的平方。相關(guān)系數(shù)的平方越高发皿，說明該網(wǎng)絡(luò)越逼近無網(wǎng)路尺度的分布崔慧。右圖的縱軸代表對應(yīng)的基因模塊中所有基因鄰接函數(shù)的均值。

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

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

大體思路：計算基因間的鄰接性，根據(jù)鄰接性計算基因間的相似性玄货，然后推出基因間的相異性系數(shù)皇钞，并據(jù)此得到基因間的系統(tǒng)聚類樹。然后按照混合動態(tài)剪切樹的標(biāo)準(zhǔn)松捉，設(shè)置每個基因模塊最少的基因數(shù)目為30夹界。根據(jù)動態(tài)剪切法確定基因模塊后，再次分析隘世，依次計算每個模塊的特征向量值掉盅，然后對模塊進(jìn)行聚類分析，將距離較近的模塊合并為新的模塊以舒。

1趾痘、計算樹之間的鄰接性

adjacency = adjacency(datExpr, power = softPower)?

2、計算樹之間的相異性

TOM = TOMsimilarity(adjacency)

dissTOM = 1-TOM

3蔓钟、聚類分析

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

4永票、設(shè)置基因模塊中最少基因包含30個

minModuleSize = 30

dynamicMods

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

2,pamRespectsDendro = FALSE, minClusterSize = minModuleSize)

5、基因分組上色

dynamicColors = labels2colors(dynamicMods)

table(dynamicColors)

6、計算基因模塊的特征值

MEList = moduleEigengenes(datExpr, colors = dynamicColors)

MEs = MEList$eigengenes

MEDiss = 1-cor(MEs)

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

7侣集、建立系統(tǒng)聚類樹

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")

}

一共會生成36個基因模塊熱圖世分，由于篇幅有限编振，僅僅展示2個

基因共表達(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)樹圖及性狀相關(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() （這個也很重要）

一步法：

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選擇模塊探測

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]);