• 鑒于大量的轉(zhuǎn)錄組數(shù)據(jù)壮锻，可能出現(xiàn)的一個問題是“如果A組和B組都進(jìn)行了轉(zhuǎn)錄組研究并報告了一些結(jié)果琐旁，這些結(jié)果的相容性如何？”目前基于轉(zhuǎn)錄組或芯片數(shù)據(jù)進(jìn)行meta-analysis的方法很多躯保，我們可以在數(shù)據(jù)庫中尋找大量類似的實驗結(jié)果用于meta-analysis旋膳，而WGCNA可以在轉(zhuǎn)錄組/芯片數(shù)據(jù)中構(gòu)建無標(biāo)度網(wǎng)絡(luò)并將基因的表達(dá)水平與樣本表型關(guān)聯(lián)起來，這在轉(zhuǎn)錄組分析中發(fā)揮了十分重要的作用途事，在這里我們可以使用WGCNA將兩組/兩組以上的數(shù)據(jù)聯(lián)合起來验懊，用于闡述不同實驗結(jié)果/不同物種間的基因表達(dá)模式的關(guān)聯(lián)

• Miller JA, Horvath S, Geschwind DH. (2010) Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc Natl Acad Sci U S A. 2010 Jul 13;107(28):12698-703. 以下分析是基于該文章的精簡版

#首先安裝一些依賴包，這里有一些包年代久遠(yuǎn)可能會報錯
  install.packages(c("impute","dynamicTreeCut","qvalue","flashClust","Hmisc","WGCNA"))

隨后下載示例數(shù)據(jù)尸变，該數(shù)據(jù)包括4組樣本义图，分別是：

• datExprA1 and datExprA2 – two data sets from the Illumina human ref-12 platform
• datExprB1 and datExprB2 – data sets from the Illumina human ref-12 and Affymetrix HG-U133A platforms, respectively (datExprA1 and datExprB1 are the same).
• probesI/A – probe set IDs for human ref-12 and HG-U133A platforms
• genesI/A – gene symbols corresponding to probesI/A

數(shù)據(jù)預(yù)處理

首先篩選探針，進(jìn)行g(shù)ene_symbol的注釋以及初篩等操作

load("/path/metaAnalysisData.RData")

library(WGCNA) 
datExprB1g = (collapseRows(datExprB1,genesI,probesI))[[1]] 
datExprB2g = (collapseRows(datExprB2,genesA,probesA))[[1]]

#您需要將分析限制在兩個數(shù)據(jù)集中都表達(dá)的基因/探針上召烂。
commonProbesA = intersect (rownames(datExprA1),rownames(datExprA2)) 
datExprA1p = datExprA1[commonProbesA,] 
datExprA2p = datExprA2[commonProbesA,] 
commonGenesB = intersect (rownames(datExprB1g),rownames(datExprB2g)) 
datExprB1g = datExprB1g[commonGenesB,] 
datExprB2g = datExprB2g[commonGenesB,]

接下來我們需要評估數(shù)據(jù)集間的關(guān)聯(lián)碱工，以判斷這些數(shù)據(jù)可否進(jìn)行聯(lián)合分析,首先分別計算每個數(shù)據(jù)的的相容性如何？”目前基于轉(zhuǎn)錄組或芯片數(shù)據(jù)進(jìn)行meta-analysis的方法很多奏夫，我們可以在數(shù)據(jù)庫中尋找大量類似的實驗結(jié)果用于meta-analysis怕篷，而WGCNA可以在轉(zhuǎn)錄組/芯片數(shù)據(jù)中構(gòu)建無標(biāo)度網(wǎng)絡(luò)并將基因的表達(dá)水平與樣本表型關(guān)聯(lián)起來，這在轉(zhuǎn)錄組分析中發(fā)揮了十分重要的作用酗昼，在這里我們可以使用WGCNA將兩組/兩組以上的數(shù)據(jù)聯(lián)合起來廊谓，用于闡述不同實驗結(jié)果/不同物種間的基因表達(dá)模式的關(guān)聯(lián)

• Miller JA, Horvath S, Geschwind DH. (2010) Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc Natl Acad Sci U S A. 2010 Jul 13;107(28):12698-703. 以下分析是基于該文章的精簡版

#首先安裝一些依賴包，這里有一些包年代久遠(yuǎn)可能會報錯
  install.packages(c("impute","dynamicTreeCut","qvalue","flashClust","Hmisc","WGCNA"))

數(shù)據(jù)預(yù)處理

首先篩選探針麻削，進(jìn)行g(shù)ene_symbol的注釋以及初篩等操作

load("/path/metaAnalysisData.RData")

library(WGCNA) 
library(flashClust)
datExprB1g = (collapseRows(datExprB1,genesI,probesI))[[1]] 
datExprB2g = (collapseRows(datExprB2,genesA,probesA))[[1]]

#您需要將分析限制在兩個數(shù)據(jù)集中都表達(dá)的基因/探針上蒸痹。
commonProbesA = intersect (rownames(datExprA1),rownames(datExprA2)) 
datExprA1p = datExprA1[commonProbesA,] 
datExprA2p = datExprA2[commonProbesA,] 
commonGenesB = intersect (rownames(datExprB1g),rownames(datExprB2g)) 
datExprB1g = datExprB1g[commonGenesB,] 
datExprB2g = datExprB2g[commonGenesB,]

接下來我們需要評估數(shù)據(jù)集間的關(guān)聯(lián)，以判斷這些數(shù)據(jù)可否進(jìn)行聯(lián)合分析,首先分別計算每個節(jié)點的加權(quán)網(wǎng)絡(luò)連接度呛哟，隨后繪制相關(guān)性散點圖查看相關(guān)性與P值判斷數(shù)據(jù)是否可以聯(lián)合分析叠荠，B組數(shù)據(jù)由于來自不同的平臺，數(shù)據(jù)間的相關(guān)性比較低扫责，但是P值很顯著(P = 2.9e-84)我們也可以將其用于后續(xù)的分析

softPower = 10 
rankExprA1= rank(rowMeans(datExprA1p)) 
rankExprA2= rank(rowMeans(datExprA2p)) 
random5000= sample(commonProbesA,5000) 
rankConnA1= rank(softConnectivity(t(datExprA1p[random5000,]),type="signed",power=softPower)) 
rankConnA2= rank(softConnectivity(t(datExprA2p[random5000,]),type="signed",power=softPower)) 
rankExprB1= rank(rowMeans(datExprB1g)) 
rankExprB2= rank(rowMeans(datExprB2g)) 
random5000= sample(commonGenesB,5000) 
rankConnB1= rank(softConnectivity(t(datExprB1g[random5000,]),type="signed",power=softPower)) 
rankConnB2= rank(softConnectivity(t(datExprB2g[random5000,]),type="signed",power=softPower))

par(mfrow=c(2,2)) 
verboseScatterplot(rankExprA1,rankExprA2, xlab="Ranked Expression (A1)", 
ylab="Ranked Expression (A2)") 
verboseScatterplot(rankConnA1,rankConnA2, xlab="Ranked Connectivity (A1)", 
ylab="Ranked Connectivity (A2)") 
verboseScatterplot(rankExprB1,rankExprB2, xlab="Ranked Expression (B1)", 
ylab="Ranked Expression (B2)") 
verboseScatterplot(rankConnB1,rankConnB2, xlab="Ranked Connectivity (B1)", 
ylab="Ranked Connectivity (B2)")

構(gòu)建共表達(dá)網(wǎng)絡(luò)

出于計算的原因和簡單性榛鼎，我們首先選擇數(shù)據(jù)集A1中表達(dá)量最高的5000個探針（通常不會這樣做），然后每個基因只保留1個探針（如上所述）鳖孤，總共保留4746個基因借帘。

keepGenesExpr = rank(-rowMeans(datExprA1p))<=5000 
datExprA1g = datExprA1p[keepGenesExpr,] 
datExprA2g = datExprA2p[keepGenesExpr,] 
keepGenesDups = (collapseRows(datExprA1g,genesI,probesI))[[2]] 
datExprA1g = datExprA1g[keepGenesDups[,2],] 
datExprA2g = datExprA2g[keepGenesDups[,2],] 
rownames(datExprA1g)<-rownames(datExprA2g)<-keepGenesDups[,1]

隨后計算連接度并構(gòu)建TOM矩陣

adjacencyA1 = adjacency(t(datExprA1g),power=softPower,type="signed"); 
diag(adjacencyA1)=0 
dissTOMA1 = 1-TOMsimilarity(adjacencyA1, TOMType="signed") 
geneTreeA1 = flashClust(as.dist(dissTOMA1), method="average") 
adjacencyA2 = adjacency(t(datExprA2g),power=softPower,type="signed"); 
diag(adjacencyA2)=0 
dissTOMA2 = 1-TOMsimilarity(adjacencyA2, TOMType="signed") 
geneTreeA2 = flashClust(as.dist(dissTOMA2), method="average")

#繪制基因?qū)哟尉垲悩?par(mfrow=c(1,2)) 
plot(geneTreeA1,xlab="",sub="",main="Gene clustering on TOM-based dissimilarity (A1)", 
labels=FALSE,hang=0.04); 
plot(geneTreeA2,xlab="",sub="",main="Gene clustering on TOM-based dissimilarity (A2)", 
labels=FALSE,hang=0.04);

接下來，我們將根據(jù)數(shù)據(jù)集A1確定模塊（在實際情況下淌铐，我們通常會構(gòu)建對照數(shù)據(jù)集的模塊作為參照）。我們使用此代碼自動顯示四個不同的模塊拆分蔫缸，可以從中進(jìn)行選擇腿准。

mColorh=NULL 
for (ds in 0:3){ 
 tree = cutreeHybrid(dendro = geneTreeA1, pamStage=FALSE, 
 minClusterSize = (30-3*ds), cutHeight = 0.99, 
 deepSplit = ds, distM = dissTOMA1) 
 mColorh=cbind(mColorh,labels2colors(tree$labels)); 
} 
#pdf("Module_choices.pdf", height=10,width=25); 
plotDendroAndColors(geneTreeA1, mColorh, paste("dpSplt =",0:3), main = "",dendroLabels=FALSE); 
#dev.off() 
modulesA1 = mColorh[,1] # (Chosen based on plot below)
#詳細(xì)過程我們需要參照cutreeHybrid這個函數(shù)，我們提供了4個整數(shù)。使用四個等級的群集拆分敏感度的粗略控制吐葱。值越高街望，將產(chǎn)生越來越多的小簇。

• 在這里選擇使用deepsplit=0（頂行）弟跑，這樣就會有少量的大模塊灾前。根據(jù)分析的目的，有時最好選擇更多的小模塊孟辑。在這種情況下哎甲，將選擇1-3的深度分割值。

• 接下來饲嗽，我們計算每個模塊的主成分炭玫。第一個PC被稱為模塊的特征基因（module eigengene [ME]），表示模塊中所有基因的最大方差百分比貌虾。換句話說吞加，如果我們展示模塊X的ME做某事，那么模塊X中的大多數(shù)基因也很有可能做同樣的事情尽狠。

PCs1A = moduleEigengenes(t(datExprA1g), colors=modulesA1)
ME_1A = PCs1A$eigengenes
distPC1A = 1-abs(cor(ME_1A,use="p"))
distPC1A = ifelse(is.na(distPC1A), 0, distPC1A)
pcTree1A = hclust(as.dist(distPC1A),method="a")
MDS_1A = cmdscale(as.dist(distPC1A),2)
colorsA1 = names(table(modulesA1))

#繪圖展示模塊之間關(guān)聯(lián)程度
plot(pcTree1A, xlab="",ylab="",main="",sub="")
plot(MDS_1A, col= colorsA1, main="MDS plot", cex=2, pch=19)
ordergenes = geneTreeA1$order
plot.mat(scale(log(datExprA1g[ordergenes,])) , rlabels= modulesA1[ordergenes], clabels=
           colnames(datExprA1g), rcols=modulesA1[ordergenes]) 
for (which.module in names(table(modulesA1))){
  ME = ME_1A[, paste("ME",which.module, sep="")]
  barplot(ME, col=which.module, main="", cex.main=2,
          ylab="eigengene expression",xlab="array sample")
}; 

#現(xiàn)在我們已經(jīng)有了所有的WGCNA變量和模塊定義衔憨，我們可以開始評估網(wǎng)絡(luò)A1中#的模塊在網(wǎng)絡(luò)A2中的情況。作為定性評估袄膏，我們將A1中的模塊強加給數(shù)據(jù)集A2
#的網(wǎng)絡(luò)践图，然后繪制出結(jié)果網(wǎng)絡(luò)。
plotDendroAndColors(geneTreeA1, modulesA1, "Modules", dendroLabels=FALSE, hang=0.03, addGuide=TRUE, 
guideHang=0.05, main="Gene dendrogram and module colors (A1)") 
plotDendroAndColors(geneTreeA2, modulesA1, "Modules", dendroLabels=FALSE, hang=0.03, addGuide=TRUE, 
guideHang=0.05, main="Gene dendrogram and module colors (A2)")

• 我們可以從這些模塊標(biāo)簽在A2的情況發(fā)現(xiàn)哩陕，它們具有非常好的保存性平项。需要注意的是，在第二個數(shù)據(jù)集中通常不可能看到明顯的模型標(biāo)簽分組悍及，即使存在顯著的數(shù)據(jù)關(guān)聯(lián)以及模塊保留能力闽瓢。

• 為了量化這個結(jié)果，我們利用WGCNA庫中內(nèi)置的另一個函數(shù)：modulePreservation心赶。此函數(shù)使用多種策略評估一項研究中的模塊在另一項研究中的保存情況扣讼，并輸出一個Z分?jǐn)?shù)摘要。

• 一般來說缨叫，Zsummary.pres的值越高椭符，模塊在數(shù)據(jù)集之間的保存程度越高：5<Z<10表示中度保存，而Z>10表示高度保存耻姥∠郏灰色模塊包含未聚類的基因，而金色模塊包含隨機基因琐簇。一般來說蒸健，這些模塊的Z分?jǐn)?shù)應(yīng)該比大多數(shù)其他模塊低座享。在這種情況下，我們發(fā)現(xiàn)所有模塊都保存得很好似忧。

multiExpr = list(A1=list(data=t(datExprA1g)),A2=list(data=t(datExprA2g))) 
multiColor = list(A1 = modulesA1) 
mp=modulePreservation(multiExpr,multiColor,referenceNetworks=1,verbose=3,networkType="signed", 
nPermutations=30,maxGoldModuleSize=100,maxModuleSize=400) 
stats = mp$preservation$Z$ref.A1$inColumnsAlsoPresentIn.A2 
stats[order(-stats[,2]),c(1:2)]

計算模塊中的基因與ME的關(guān)聯(lián)程度

Module membership (kME) 是一個有用的值渣叛，因為它可以用來測量每個基因和每個ME之間的相關(guān)性，因此即使是最初沒有分配到模塊的基因也可以包含在網(wǎng)絡(luò)間的比較中盯捌。

geneModuleMembership1 = signedKME(t(datExprA1g), ME_1A) 
colnames(geneModuleMembership1)=paste("PC",colorsA1,".cor",sep=""); 
MMPvalue1=corPvalueStudent(as.matrix(geneModuleMembership1),dim(datExprA1g)[[2]]); 
colnames(MMPvalue1)=paste("PC",colorsA1,".pval",sep=""); 
Gene = rownames(datExprA1g) 
kMEtable1 = cbind(Gene,Gene,modulesA1) 
for (i in 1:length(colorsA1)) 
 kMEtable1 = cbind(kMEtable1, geneModuleMembership1[,i], MMPvalue1[,i]) 
colnames(kMEtable1)=c("PSID","Gene","Module",sort(c(colnames(geneModuleMembership1), 
colnames(MMPvalue1))))

#現(xiàn)在對A2重復(fù)上述步驟淳衙，使用來自A1的模塊賦值來確定A2中的kME值。
# First calculate MEs for A2, since we haven't done that yet 
PCs2A = moduleEigengenes(t(datExprA2g), colors=modulesA1) 
ME_2A = PCs2A$eigengenes 
geneModuleMembership2 = signedKME(t(datExprA2g), ME_2A) 
colnames(geneModuleMembership1)=paste("PC",colorsA1,".cor",sep=""); 
MMPvalue2=corPvalueStudent(as.matrix(geneModuleMembership2),dim(datExprA2g)[[2]]); 
colnames(MMPvalue2)=paste("PC",colorsA1,".pval",sep=""); 
kMEtable2 = cbind(Gene,Gene,modulesA1) 
for (i in 1:length(colorsA1)) 
 kMEtable2 = cbind(kMEtable2, geneModuleMembership2[,i], MMPvalue2[,i]) 
colnames(kMEtable2)=colnames(kMEtable1)

• 我們首先要做的是將A1中每個基因的kME值與A2中每個基因相應(yīng)的kME值進(jìn)行比較饺著。具有高相關(guān)性的點的模塊被高度保留箫攀。下面是所有基因（左）和被分配進(jìn)該模塊的基因子集（右）創(chuàng)建這些圖的示例圖像。
• 我們可以做的第二件事是通過確定哪些基因在兩個網(wǎng)絡(luò)中都具有極高的kME值來確定哪些基因是兩個網(wǎng)絡(luò)中的樞紐瓶籽。

#pdf("all_kMEtable2_vs_kMEtable1.pdf",height=8,width=8) 
for (c in 1:length(colorsA1)){ 
 verboseScatterplot(geneModuleMembership2[,c],geneModuleMembership1[,c],main=colorsA1[c], 
 xlab="kME in A2",ylab="kME in A1") 
}; 
#pdf("inModule_kMEtable2_vs_kMEtable1.pdf",height=8,width=8) 
for (c in 1:length(colorsA1)){ 
 inMod = modulesA1== colorsA1[c] 
 verboseScatterplot(geneModuleMembership2[inMod,c],geneModuleMembership1[inMod,c],main=colorsA1[c], 
 xlab="kME in A2",ylab="kME in A1") 
}; 

topGenesKME = NULL 
for (c in 1:length(colorsA1)){ 
 kMErank1 = rank(-geneModuleMembership1[,c]) 
 kMErank2 = rank(-geneModuleMembership2[,c]) 
 maxKMErank = rank(apply(cbind(kMErank1,kMErank2+.00001),1,max)) 
 topGenesKME = cbind(topGenesKME,Gene[maxKMErank<=10]) 
}; 
colnames(topGenesKME) = colorsA1

這兩種類型的散點圖傳達(dá)了相似但不完全相同的信息匠童。使用所有基因（左）可以包含所有正相關(guān)和負(fù)相關(guān)的基因，但通常也包含很多噪聲（盡管在本例中沒有）塑顺。只使用模塊內(nèi)基因（右）是評估hub基因在不同數(shù)據(jù)集間的保守程度：如果這些基因顯示出集合間的相關(guān)性汤求，那么位于圖右上角的基因很可能是數(shù)據(jù)集之間的公共hub基因。（Hub基因是與MEs和模塊內(nèi)高連通性有顯著相關(guān)性的基因）

注釋模塊以及后續(xù)分析

為了可視化生成的網(wǎng)絡(luò)模塊严拒，我們可以從網(wǎng)絡(luò)輸出數(shù)據(jù)導(dǎo)入VisANT扬绪，VisANT是一個獨立的可視化java程序，可在http://VisANT.bu.edu/上獲得裤唠。VisANT可以讓你直觀地看到一個模塊中的中心基因挤牛，是大多數(shù)人用來制作模塊圖片的程序。關(guān)于如何使用VisANT的教程可以在網(wǎng)站上找到种蘸。我們可以使用作者封裝的代碼提取模塊中的基因用于VisANT的輸入文件

source("tutorialFunctions.R") 
for (co in colorsA1[colorsA1!="grey"])
 visantPrepOverall(modulesA1, co, t(datExprA1g), rownames(datExprA1g), 500, softPower, TRUE)

for (co in colorsA1[colorsA1!="grey"]) 
 visantPrepOverall(modulesA1, co, t(datExprA2g), rownames(datExprA2g), 500, softPower, TRUE)

• 這個函數(shù)會將很多文件輸出到當(dāng)前目錄中墓赴，每個模塊都有兩個文件：
  1） <module>_connectivityOverall.csv:此文件包含給定模塊中從最高到最低模塊內(nèi)連接（kin）排序的所有基因及其平均表達(dá)的列表。所以排在第一位的基因是中樞基因航瞭。
  2） <module>_visantOverall.csv:此文件包含VisANT所需的所有輸入诫硕。注意，只有前五列應(yīng)該復(fù)制到VisANT中刊侯。前兩列表示模塊中拓?fù)渲丿B度最高的基因（第5列）章办。第3列必須是“0”，第4列是要繪制的線的顏色（“M1002”是橙色）

• 通過該軟件我們可以查看不同數(shù)據(jù)之間不同的表達(dá)模式滨彻，該示例中我們可以很清楚的發(fā)現(xiàn)SCAMP5是A1特異表達(dá)的藕届，我們可以選擇這樣的基因進(jìn)行后續(xù)的分析