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

#首先安裝一些依賴包，這里有一些包年代久遠可能會報錯
  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ù)處理

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

load("/path/metaAnalysisData.RData")

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

#您需要將分析限制在兩個數(shù)據(jù)集中都表達的基因/探針上。
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ù)可否進行聯(lián)合分析,首先分別計算每個數(shù)據(jù)的的相容性如何？”目前基于轉(zhuǎn)錄組或芯片數(shù)據(jù)進行meta-analysis的方法很多，我們可以在數(shù)據(jù)庫中尋找大量類似的實驗結(jié)果用于meta-analysis，而WGCNA可以在轉(zhuǎn)錄組/芯片數(shù)據(jù)中構(gòu)建無標(biāo)度網(wǎng)絡(luò)并將基因的表達水平與樣本表型關(guān)聯(lián)起來，這在轉(zhuǎn)錄組分析中發(fā)揮了十分重要的作用，在這里我們可以使用WGCNA將兩組/兩組以上的數(shù)據(jù)聯(lián)合起來，用于闡述不同實驗結(jié)果/不同物種間的基因表達模式的關(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. 以下分析是基于該文章的精簡版

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

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

首先篩選探針，進行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ù)集中都表達的基因/探針上。
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ù)可否進行聯(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)建共表達網(wǎng)絡(luò)

出于計算的原因和簡單性，我們首先選擇數(shù)據(jù)集A1中表達量最高的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ù)集的模塊作為參照）。我們使用此代碼自動顯示四個不同的模塊拆分，可以從中進行選擇。

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)
#詳細過程我們需要參照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分數(shù)摘要。
一般來說，Zsummary.pres的值越高，模塊在數(shù)據(jù)集之間的保存程度越高：5<Z<10表示中度保存，而Z>10表示高度保存?；疑K包含未聚類的基因，而金色模塊包含隨機基因。一般來說，這些模塊的Z分數(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值進行比較。具有高相關(guā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

這兩種類型的散點圖傳達了相似但不完全相同的信息。使用所有基因（左）可以包含所有正相關(guān)和負相關(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）排序的所有基因及其平均表達的列表。所以排在第一位的基因是中樞基因。
2） <module>_visantOverall.csv:此文件包含VisANT所需的所有輸入。注意，只有前五列應(yīng)該復(fù)制到VisANT中。前兩列表示模塊中拓撲重疊度最高的基因（第5列）。第3列必須是“0”，第4列是要繪制的線的顏色（“M1002”是橙色）
通過該軟件我們可以查看不同數(shù)據(jù)之間不同的表達模式，該示例中我們可以很清楚的發(fā)現(xiàn)SCAMP5是A1特異表達的，我們可以選擇這樣的基因進行后續(xù)的分析

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

使用WGCNA進行meta-analysis

使用WGCNA進行meta-analysis

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

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

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

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

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

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

使用WGCNA進行meta-analysis

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

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

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

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

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

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av