簡(jiǎn)介
WGCNA網(wǎng)絡(luò)分析以后骡技，我們可以對(duì)一些sample構(gòu)建出共表達(dá)的網(wǎng)絡(luò)鸣个，那么這里有一個(gè)問(wèn)題，那就是不同處理下的sample（比方說(shuō)treat和control）布朦，它們的共表達(dá)網(wǎng)絡(luò)是一樣的嗎囤萤？如果不一樣，那些產(chǎn)生差異的gene pair是否與這些處理有關(guān)系是趴，是否是一種比較重要的核心基因呢阁将？
本篇推文老藥新用，這是一篇2010年發(fā)表的文章：《DiffCoEx: a simple and sensitive method to find differentially coexpressed gene modules》
這篇文章基于傳統(tǒng)的WGCNA構(gòu)建的網(wǎng)絡(luò)右遭，來(lái)尋找兩種處理的差異共表達(dá)網(wǎng)絡(luò)
數(shù)學(xué)原理
構(gòu)建差分網(wǎng)絡(luò)大致分為5個(gè)步驟：
step 1：
建立鄰接矩陣C，鄰接矩陣是由兩個(gè)gene pair的相關(guān)性來(lái)表征的缤削，即兩個(gè)基因之間的權(quán)重
step 2：
計(jì)算gene pair在兩個(gè)組別權(quán)重的變化程度窘哈，并構(gòu)建鄰接變化矩陣 D；處理組為 [1]亭敢，對(duì)照組為 [2]
其中:
- C[1]ij 代表處理組的鄰接矩陣滚婉；C[2]ij 代表對(duì)照組的鄰接矩陣
- β 為WGCNA的計(jì)算鄰接矩陣的指數(shù) β
- sign() 為符號(hào)函數(shù)
如果 dij 越高，代表在 A 處理中 gene i 和 gene j 權(quán)重（相關(guān)性帅刀，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性让腹，調(diào)控關(guān)系）的差值越大，即 gene i 與 gene j 的權(quán)重（相關(guān)性扣溺，調(diào)控關(guān)系）在兩個(gè)處理組之間具有顯著差異
step 3：
由鄰接變化矩陣 D 計(jì)算相異矩陣 T
其中骇窍，k 代表與 gene i 和 gene j 不同的 gene k
那么相異矩陣 T 考慮了第三方gene k所帶來(lái)的影響，那么tij 越低锥余，代表在 A 處理中 gene i 和 gene j 權(quán)重（相關(guān)性腹纳，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性，調(diào)控關(guān)系）的差值越大，即 gene i 與 gene j 的權(quán)重（相關(guān)性嘲恍，調(diào)控關(guān)系）在兩個(gè)處理組之間具有顯著差異
理解：假設(shè)C[1]ij很高足画，C[2]ij很低，則 dij 很高佃牛，代表在 A 處理中 gene i 和 gene j 權(quán)重（相關(guān)性淹辞，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性，調(diào)控關(guān)系）的差值越大俘侠，即 gene i 與 gene j 的權(quán)重（相關(guān)性象缀，調(diào)控關(guān)系）在兩個(gè)處理組之間具有顯著差異，此時(shí)兼贡，tij 值就越低
step 4：
進(jìn)行樹(shù)聚類(lèi)攻冷，由step3中可知，tij 值小代表在 A 處理中 gene i 和 gene j 權(quán)重（相關(guān)性遍希，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性等曼，調(diào)控關(guān)系）有顯著差異；即代表在 A 處理中 gene i 和 gene j 權(quán)重（相關(guān)性凿蒜，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性禁谦，調(diào)控關(guān)系）的差值越大，
step 5：
顯著性統(tǒng)計(jì)
代碼：
作者以一套R(shí)NA探針表達(dá)譜的數(shù)據(jù)為例子：
數(shù)據(jù)在這里下載：http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS2901
#required libraries
#WGCNA package can be found at
#http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA
library(WGCNA) ###used for topological overlap calculation and clustering steps
library(RColorBrewer) ###used to create nicer colour palettes
library(preprocessCore) ###used by the quantile normalization function
library(flashClust)
#Note: the data can be downloaded from the Gene Expression Omnibus
# http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS2901
data<-as.matrix(read.csv(file="GDS2901.soft",skip=166,row.names=1,sep="\t",header=T))
data<-data[-15924,]
rawData<-matrix(as.numeric(data[,-1]),nrow=15923)
dimnames(rawData)<-dimnames(data[,-1])
#we create an annotation matrix containing the matches between probesets and gene names
anno<-as.matrix(data[-2475,1])
normData<-normalize.quantiles(log2(rawData))
dimnames(normData)<-dimnames(rawData)
#we remove the probeset at index 2475 because
#after quantile normalization it has zero variance
#(the probeset has the highest signal of all samples)
normData<-normData[-2475,]
# 定義組別 1 的鄰接矩陣
datC1<-t(normData[,c(1:12,25:36,37:48)]) ### these samples correspond to the Eker mutants.
# Note that since the Eker mutants have two sets of 12 control samples (13:24 and 37:48)
# we discard one to have a symmetric perturbation (carcinogenic vs control) between the two conditions (Eker mutants vs wild-types)
# 定義組別 2 的鄰接矩陣
datC2<-t(normData[,49:84]) ###those samples correspond to the wild-types
# 定義 β 值
beta1=6 #user defined parameter for soft thresholding
# 可參考step 2數(shù)學(xué)原理
AdjMatC1<-sign(cor(datC1,method="spearman"))*(cor(datC1,method="spearman"))^2
AdjMatC2<-sign(cor(datC2,method="spearman"))*(cor(datC2,method="spearman"))^2
diag(AdjMatC1)<-0
diag(AdjMatC2)<-0
collectGarbage()
# 計(jì)算相異矩陣 T
dissTOMC1C2=TOMdist((abs(AdjMatC1-AdjMatC2)/2)^(beta1/2))
collectGarbage()
# 對(duì)相異矩陣 T 進(jìn)行聚類(lèi)
#Hierarchical clustering is performed using the Topological Overlap of the adjacency difference as input distance matrix
geneTreeC1C2 = flashClust(as.dist(dissTOMC1C2), method = "average");
# Plot the resulting clustering tree (dendrogram)
png(file="hierarchicalTree.png",height=1000,width=1000)
plot(geneTreeC1C2, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",labels = FALSE, hang = 0.04);
dev.off()
#We now extract modules from the hierarchical tree. This is done using cutreeDynamic. Please refer to WGCNA package documentation for details
dynamicModsHybridC1C2 = cutreeDynamic(dendro = geneTreeC1C2, distM = dissTOMC1C2,method="hybrid",cutHeight=.996,deepSplit = T, pamRespectsDendro = FALSE,minClusterSize = 20);
#Every module is assigned a color. Note that GREY is reserved for genes which do not belong to any differentially coexpressed module
dynamicColorsHybridC1C2 = labels2colors(dynamicModsHybridC1C2)
#the next step merges clusters which are close (see WGCNA package documentation)
mergedColorC1C2<-mergeCloseModules(rbind(datC1,datC2),dynamicColorsHybridC1C2,cutHeight=.2)$color
colorh1C1C2<-mergedColorC1C2
#reassign better colors
colorh1C1C2[which(colorh1C1C2 =="midnightblue")]<-"red"
colorh1C1C2[which(colorh1C1C2 =="lightgreen")]<-"yellow"
colorh1C1C2[which(colorh1C1C2 =="cyan")]<-"orange"
colorh1C1C2[which(colorh1C1C2 =="lightcyan")]<-"green"
# Plot the dendrogram and colors underneath
png(file="module_assignment.png",width=1000,height=1000)
plotDendroAndColors(geneTreeC1C2, colorh1C1C2, "Hybrid Tree Cut",dendroLabels = FALSE, hang = 0.03,addGuide = TRUE, guideHang = 0.05,main = "Gene dendrogram and module colors cells")
dev.off()
#We write each module to an individual file containing affymetrix probeset IDs
modulesC1C2Merged<-extractModules(colorh1C1C2,datC1,anno,dir="modules",file_prefix=paste("Output","Specific_module",sep=''),write=T)
write.table(colorh1C1C2,file="module_assignment.txt",row.names=F,col.names=F,quote=F)
#We plot to a file the comparative heatmap showing correlation changes in the modules
#The code for the function plotC1C2Heatmap and others can be found below under the Supporting Functions section
plotC1C2Heatmap(colorh1C1C2,AdjMatC1,AdjMatC2, datC1, datC2)
png(file="exprChange.png",height=500,width=500)
plotExprChange(datC1,datC2,colorh1C1C2)
dev.off()
進(jìn)行顯著性分析
#This function computes the dispersion value that
#quantifies the change in correlation between two conditions
#for pair of genes drawn from module c1 and module c2
# in case c1 = c2, the function quantifies the differential coexpression in c1.
#cf Choi and Kendziorski 2009
dispersionModule2Module<-function(c1,c2,datC1,datC2,colorh1C1C2)
{
if (c1==c2)
{
difCor<-(cor(datC1[,which(colorh1C1C2 == c1)],method="spearman")-
cor(datC2[,which(colorh1C1C2 == c1)],method="spearman"))^2
n<-length(which(colorh1C1C2 ==c1))
(1/((n^2 -n)/2)*(sum(difCor)/2))^(.5)
}
else if (c1!=c2)
{
difCor<-(cor(datC1[,which(colorh1C1C2 == c1)],datC1[,which(colorh1C1C2==c2)],method="spearman")-
cor(datC2[,which(colorh1C1C2 == c1)],datC2[,which(colorh1C1C2==c2)],method="spearman"))^2
n1<-length(which(colorh1C1C2 ==c1))
n2<-length(which(colorh1C1C2 ==c2))
(1/((n1*n2))*(sum(difCor)))^(.5)
}
}
# we generate a set of 1000 permuted indexes
permutations<-NULL
for (i in 1:1000)
{
permutations<-rbind(permutations,sample(1:(nrow(datC1)+nrow(datC2)),nrow(datC1)))
}
# we scale the data in both conditions to mean 0 and variance 1.
d<-rbind(scale(datC1),scale(datC2))
# This function calculates the dispersion value of a module to module coexpression change on permuted data
permutationProcedureModule2Module<-function(permutation,d,c1,c2,colorh1C1C2)
{
d1<-d[permutation,]
d2<-d[-permutation,]
dispersionModule2Module(c1,c2,d1,d2,colorh1C1C2)
}
#We compute all pairwise module to module dispersion values, and generate a null distribution from permuted scaled data
dispersionMatrix<-matrix(nrow=length(unique(colorh1C1C2))-1,ncol=length(unique(colorh1C1C2))-1)
nullDistrib<-list()
i<-j<-0
for (c1 in setdiff(unique(colorh1C1C2),"grey"))
{
i<-i+1
j<-0
nullDistrib[[c1]]<-list()
for (c2 in setdiff(unique(colorh1C1C2),"grey"))
{
j<-j+1
dispersionMatrix[i,j]<-dispersionModule2Module(c1,c2,datC1,datC2,colorh1C1C2)
nullDistrib[[c1]][[c2]]<-apply(permutations,1,permutationProcedureModule2Module,d,c2,c1,colorh1C1C2)
}
}
#We create a summary matrix indicating for each module to module
#differential coexpression the number of permuted data yielding
#an equal or higher dispersion.
permutationSummary<-matrix(nrow=8,ncol=8)
colnames(permutationSummary)<-setdiff(unique(colorh1C1C2),"grey")
rownames(permutationSummary)<-setdiff(unique(colorh1C1C2),"grey")
for (i in 1:8) {
for (j in 1:8) {
permutationSummary[i,j]<-length(which(nullDistrib[[i]][[j]] >= dispersionMatrix[i,j]))
}
}
#We plot the result (cf supplementary figure 1)
plotMatrix(permutationSummary)
關(guān)于相異矩陣 dissTOMC1C2 的理解：
gene_1 與 gene_2 的系數(shù)為 0.9731268废封，代表在 A 處理中 gene 1 和 gene 2 權(quán)重（相關(guān)性州泊，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性，調(diào)控關(guān)系）的相異系數(shù) tij 為 0.9731268（相異系數(shù) tij 越大代表：在 A 處理中 gene i 和 gene j 權(quán)重（相關(guān)性漂洋，調(diào)控關(guān)系）與 在 B 處理中 gene i 和 gene j 權(quán)重（相關(guān)性遥皂，調(diào)控關(guān)系）差異越小刽漂；反之越大）
