Steve Horvath UCLA – Horvath Lab UCLA
[Pages:1]
Is human blood a good surrogate for brain tissue in transcriptional studies?
Chaochao Cai, Peter Langfelder, Tova F Fuller, Michael C Oldham, Rui Luo, Leonard H van den Berg, Roel A Ophoff, 4, Steve Horvath
Correspondence: shorvath@mednet.ucla.edu
This R tutorial describes how to carry out a preservation analysis between brain and blood with the R software. The analysis requires WGCNA package, which is available at
#Preprocess the data
#Since we try to find which brain module is preserved in blood, genes expressed in both #tissues are selected out for future analyisis. We do it region by region.
#For CTX and blood
#load the WGCNA package
library(WGCNA)
#load the data
load("CTX_datExpr_module.RData")
load("SAFHS_datExpr.RData")
load("Dutch_datExpr.RData")
#find overlap genes by symbol name
CTX.genes= rownames(datExpr.CTX)
SAFHS.genes= rownames(datExpr.SAFHS)
Dutch.genes= rownames(datExpr.Dutch)
shared.genes =merge(merge(CTX.genes, SAFHS.genes,by=1,sort=F)[,1], rownames(datExpr.Dutch),by=1,sort=F)[,1]
#since some of them are not expressed in blood, we try to move unexpressed genes
shared.genes = as.character(merge(shared.genes, Expressed_blood,by=1,sort=F)[,1])
#finally, we get 2640 expressed genes shared between cortex and blood.
length(shared.genes)
[1] 2640
#Get the datExpr (only contains the shared genes) and module definition for next analysis.
data.CTX=datExpr.CTX[match(shared.genes, CTX.genes),]
data.SAFHS=datExpr.SAFHS[match(shared.genes, SAFHS.genes),]
data.Dutch=datExpr.Dutch[match(shared.genes, Dutch.genes),]
color.CTX=module.CTX[match(shared.genes, CTX.genes)]
#save the files
save(data.CTX,data.SAFHS,data.Dutch,color.CTX,file="CTX_SAFHS_Dutch_for_MP.RData")
#For CN and blood
#load the data
load("CN_datExpr_module.RData")
load("SAFHS_datExpr.RData")
load("Dutch_datExpr.RData")
#find overlap genes by symbol name
CN.genes= rownames()
SAFHS.genes= rownames(datExpr.SAFHS)
Dutch.genes= rownames(datExpr.Dutch)
shared.genes =merge(merge(CN.genes, SAFHS.genes,by=1,sort=F)[,1], rownames(datExpr.Dutch),by=1,sort=F)[,1]
#since some of them are not expressed in blood, we try to move unexpressed genes
shared.genes = as.character(merge(shared.genes, Expressed_blood,by=1,sort=F)[,1])
#finally, we get 2063 expressed genes shared between caudate nucleus and blood.
length(shared.genes)
[1] 2063
#Get the datExpr (only contains the shared genes) and module definition for next analysis.
=[match(shared.genes, CN.genes),]
data.SAFHS=datExpr.SAFHS[match(shared.genes, SAFHS.genes),]
data.Dutch=datExpr.Dutch[match(shared.genes, Dutch.genes),]
= [match(shared.genes, CN.genes)]
#save the files
save(,data.SAFHS,data.Dutch,,file="CN_SAFHS_Dutch_for_MP.RData")
#For CB and blood
#load the data
load("CB_datExpr_module.RData")
load("SAFHS_datExpr.RData")
load("Dutch_datExpr.RData")
#find overlap genes by symbol name
CB.genes= rownames(datExpr.CB)
SAFHS.genes= rownames(datExpr.SAFHS)
Dutch.genes= rownames(datExpr.Dutch)
shared.genes =merge(merge(CB.genes, SAFHS.genes,by=1,sort=F)[,1], rownames(datExpr.Dutch),by=1,sort=F)[,1]
#since some of them are not expressed in blood, we try to move unexpressed genes
shared.genes = as.character(merge(shared.genes, Expressed_blood,by=1,sort=F)[,1])
#finally, we get 2001 expressed genes shared between cerebellum and blood.
length(shared.genes)
[1] 2001
#Get the datExpr (only contains the shared genes) and module definition for next analysis.
data.CB=datExpr.CB[match(shared.genes, CB.genes),]
data.SAFHS=datExpr.SAFHS[match(shared.genes, SAFHS.genes),]
data.Dutch=datExpr.Dutch[match(shared.genes, Dutch.genes),]
color.CB=module.CB [match(shared.genes, CB.genes)]
#save the files
save(data.CB,data.SAFHS,data.Dutch,color.CB,file="CB_SAFHS_Dutch_for_MP.RData")
# Comparison of mean expression level and connectivity between brain and blood
#first, explores the preservation of mean expression between brain and blood, and plots the results. We compare each brain region with both blood data, so we have 6 plots finally.
#load WGCNA package
library(WGCNA)
pdf(8,12,file="Mean_brain_blood_Consistance.pdf")
par(mfrow=c(3,2))
par(mar=c(7,8,4,1))
par(mgp = c(4,1,0))
#for CTX data
load("CTX_SAFHS_Dutch_for_MP.RData")
mean.CTX=apply(data.CTX,1,mean)
mean.SAFHS=apply(data.SAFHS,1,mean)
mean.Dutch=apply(data.Dutch,1,mean)
#plot the result between CTX and Dutch
verboseScatterplot(mean.CTX,mean.Dutch,xlab="Mean expression\nCTX", ylab=" Mean expression\nDutch",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#plot the result between CTX and SAFHS
verboseScatterplot(mean.CTX,mean.SAFHS,xlab="Mean expression\nCTX", ylab=" Mean expression\nSAFHS",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#for CN data
rm(list=ls())
load("CN_SAFHS_Dutch_for_MP.RData")
mean.SAFHS=apply(data.SAFHS,1,mean)
mean.Dutch=apply(data.Dutch,1,mean)
=apply(,1,mean)
#plot the result between CN and Dutch
verboseScatterplot(,mean.Dutch,xlab="Mean expression\nCN", ylab=" Mean expression\nDutch",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#plot the result between CN and SAFHS
verboseScatterplot(,mean.SAFHS,xlab="Mean expression\nCN", ylab=" Mean expression\nSAFHS",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "d", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#for CB data
rm(list=ls())
load("CB_SAFHS_Dutch_for_MP.RData")
mean.SAFHS=apply(data.SAFHS,1,mean)
mean.Dutch=apply(data.Dutch,1,mean)
mean.CB=apply(data.CB,1,mean)
#plot the result between CB and Dutch
verboseScatterplot(mean.CB,mean.Dutch,xlab="Mean expression\nCB", ylab=" Mean expression\nDutch",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "e", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#plot the result between CB and SAFHS
verboseScatterplot(mean.CB,mean.SAFHS,xlab="Mean expression\nCB", ylab=" Mean expression\nSAFHS",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "f", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
dev.off()
# second, explores the preservation of connectivity between brain and blood, and plots the results
pdf(8,12,file=" Connectivity _brain_blood_Consistance.pdf")
par(mfrow=c(3,2))
par(mar=c(7,8,4,1))
par(mgp = c(4,1,0))
#for CTX data
load("CTX_SAFHS_Dutch_for_MP.RData")
sftCTX=softConnectivity(t(data.CTX))
sftSAFHS=softConnectivity(t(data.SAFHS))
sftDutch=softConnectivity(t(data.Dutch))
#plot the result between CTX and Dutch
verboseScatterplot(sftCTX,sftDutch,xlab="Connectivity\nCTX", ylab=" Connectivity\nDutch",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#plot the result between CTX and SAFHS
verboseScatterplot(sftCTX,sftSAFHS,xlab="Connectivity\nCTX", ylab=" Connectivity\nSAFHS",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#for CN data
rm(list=ls())
load("CN_SAFHS_Dutch_for_MP.RData")
sftCN=softConnectivity(t())
sftSAFHS=softConnectivity(t(data.SAFHS))
sftDutch=softConnectivity(t(data.Dutch))
#plot the result between CN and Dutch
verboseScatterplot(sftCN,sftDutch,xlab="Connectivity\nCN", ylab=" Connectivity\nDutch",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#plot the result between CN and SAFHS
verboseScatterplot(sftCN,sftSAFHS,xlab="Connectivity\nCN", ylab=" Connectivity\nSAFHS",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "d", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#for CB data
rm(list=ls())
load("CB_SAFHS_Dutch_for_MP.RData")
sftCB=softConnectivity(t(data.CB))
sftSAFHS=softConnectivity(t(data.SAFHS))
sftDutch=softConnectivity(t(data.Dutch))
#plot the result between CB and Dutch
verboseScatterplot(sftCB,sftDutch,xlab="Connectivity\nCB", ylab=" Connectivity\nDutch",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "e", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#plot the result between CB and SAFHS
verboseScatterplot(sftCB,sftSAFHS,xlab="Connectivity\nCB", ylab=" Connectivity\nSAFHS",corFnc = "bicor")
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.09*(box[4]-box[3]), "f", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
dev.off()
#Preservation module detection
#explores the preservation between each brain region in blood with Modulepreservation #funtion, which is available in WGCNA package.
#MP function for CTX and SAFHS
#load WGCNA
library(WGCNA)
load("CTX_SAFHS_Dutch_for_MP.RData")
colorList=list(color.CTX)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data.CTX))
multiExpr[[2]]=list(data=t(data.SAFHS))
names(multiExpr)=names(colorList)=c("net1","net2")
mp.CTX.SAFHS=modulePreservation(multiExpr,colorList)
save(mp.CTX.SAFHS,file="CTX_SAFHS_MP_Result.RData")
#MP function for CTX and Dutch
colorList=list(color.CTX)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data.CTX))
multiExpr[[2]]=list(data=t(data.Dutch))
names(multiExpr)=names(colorList)=c("net1","net2")
mp.CTX.Dutch=modulePreservation(multiExpr,colorList)
save(mp.CTX.Dutch,file="CTX_Dutch_MP_Result.RData")
#MP function for CN and SAFHS
rm(list=ls())
load("CN_SAFHS_Dutch_for_MP.RData")
colorList=list()
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t())
multiExpr[[2]]=list(data=t(data.SAFHS))
names(multiExpr)=names(colorList)=c("net1","net2")
.SAFHS=modulePreservation(multiExpr,colorList)
save(.SAFHS,file="CN_SAFHS_MP_Result.RData")
#MP function for CN and Dutch
colorList=list()
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t())
multiExpr[[2]]=list(data=t(data.Dutch))
names(multiExpr)=names(colorList)=c("net1","net2")
.Dutch=modulePreservation(multiExpr,colorList)
save(.Dutch,file="CN_Dutch_MP_Result.RData")
#MP function for CB and SAFHS
rm(list=ls())
load("CB_SAFHS_Dutch_for_MP.RData")
colorList=list(color.CB)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data.CB))
multiExpr[[2]]=list(data=t(data.SAFHS))
names(multiExpr)=names(colorList)=c("net1","net2")
mp.CB.SAFHS=modulePreservation(multiExpr,colorList)
save(mp.CB.SAFHS,file="CB_SAFHS_MP_Result.RData")
#MP function for CN and Dutch
colorList=list(color.CB)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data.CB))
multiExpr[[2]]=list(data=t(data.Dutch))
names(multiExpr)=names(colorList)=c("net1","net2")
mp.CB.Dutch=modulePreservation(multiExpr,colorList)
save(mp.CB.Dutch,file="CB_Dutch_MP_Result.RData")
#although module preservation function returns many statistics, Z summary is a good represent ##for them, and we will use Z summary in our case.
#extracts the Z summary for CTX modules
load("CTX_Dutch_MP_Result.RData")
load("CTX_SAFHS_MP_Result.RData")
CTX.SAFHS.preservation=mp.CTX.SAFHS$preservation$Z$1$2 [,1:2]
CTX.Dutch.preservation=mp.CTX.Dutch$preservation$Z$1$ 2 [,1:2]
##remove grey and gold module
CTX.SAFHS.preservation=CTX.SAFHS.preservation[-c(5,8),]
CTX.Dutch.preservation=CTX.Dutch.preservation[-c(5,8),]
##reoder the Z-summary
k.order=order((CTX.Dutch.preservation[,2]+CTX.SAFHS.preservation[,2]),decreasing=T)
CTX.Dutch.preservation=CTX.Dutch.preservation[k.order,]
CTX.SAFHS.preservation=CTX.SAFHS.preservation[k.order,]
CTX.SAFHS.preservation[,1]=paste(rownames(CTX.SAFHS.preservation), " SAFHS",sep="/")
CTX.Dutch.preservation[,1]=paste(rownames(CTX.Dutch.preservation), "Dutch",sep="/")
colnames(CTX.Dutch.preservation)[1]="Module_name"
colnames(CTX.SAFHS.preservation)[1]="Module_name"
bine=data.frame(matrix(NA,38,3))
colnames(bine)=c("Module_Colour","Name","Preservation_Z_Score")
for(i in 1:19){ bine[2*i,1]=rownames(CTX.SAFHS.preservation)[i];
bine[2*i,2:3]=CTX.SAFHS.preservation[i,1:2]}
for(i in 1:19){ bine[2*i-1,1]=rownames(CTX.Dutch.preservation)[i];
bine[2*i-1,2:3]=CTX.Dutch.preservation[i,1:2]}
write.csv(bine,file="CTX_preservation.csv")
| |Module_Colour |Name |Preservation_Z_Score |
|1 |yellow |yellow/Dutch |17.6230047 |
|2 |yellow |yellow/ SAFHS |14.51248412 |
|3 |green |green/Dutch |15.28251201 |
|4 |green |green/ SAFHS |16.51670502 |
|5 |blue |blue/Dutch |10.93244079 |
|6 |blue |blue/ SAFHS |12.02138823 |
|7 |turquoise |turquoise/Dutch |5.462635061 |
|8 |turquoise |turquoise/ SAFHS |5.844942121 |
|9 |orange |orange/Dutch |4.625261909 |
|10 |orange |orange/ SAFHS |6.074912023 |
|11 |palegreen |palegreen/Dutch |1.867224282 |
|12 |palegreen |palegreen/ SAFHS |7.172361084 |
|13 |midnightblue |midnightblue/Dutch |4.021562421 |
|14 |midnightblue |midnightblue/ SAFHS |3.802471484 |
|15 |tan |tan/Dutch |0.758099653 |
|16 |tan |tan/ SAFHS |6.659525167 |
|17 |honeydew |honeydew/Dutch |5.301987295 |
|18 |honeydew |honeydew/ SAFHS |1.244696498 |
|19 |pink |pink/Dutch |2.295085792 |
|20 |pink |pink/ SAFHS |3.449963031 |
|21 |purple |purple/Dutch |0.619838871 |
|22 |purple |purple/ SAFHS |4.84870282 |
|23 |brown |brown/Dutch |4.072739681 |
|24 |brown |brown/ SAFHS |1.301297022 |
|25 |tomato |tomato/Dutch |0.595865764 |
|26 |tomato |tomato/ SAFHS |1.639600186 |
|27 |red |red/Dutch |0.860534257 |
|28 |red |red/ SAFHS |0.568528373 |
|29 |powderblue |powderblue/Dutch |0.253492922 |
|30 |powderblue |powderblue/ SAFHS |1.160586199 |
|31 |salmon |salmon/Dutch |0.222788429 |
|32 |salmon |salmon/ SAFHS |0.1740373 |
|33 |greenyellow |greenyellow/Dutch |0.197154624 |
|34 |greenyellow |greenyellow/ SAFHS |-0.344485201 |
|35 |darkolivegreen |darkolivegreen/Dutch |-0.567849659 |
|36 |darkolivegreen |darkolivegreen/ SAFHS |-0.747719241 |
|37 |black |black/Dutch |-0.448040811 |
|38 |black |black/ SAFHS |-1.290898012 |
#extracts the preservation result for CN modules
#load the data
rm(list=ls())
load("CN_Dutch_MP_Result.RData")
load("CN_SAFHS_MP_Result.RData")
CN.SAFHS.preservation=.SAFHS$preservation$Z$1$2 [,1:2]
CN.Dutch.preservation=.Dutch$preservation$Z$1$ 2 [,1:2]
##remove grey and gold module
CN.Dutch.preservation=CN.Dutch.preservation[-c(8,9),]
CN.SAFHS.preservation=CN.SAFHS.preservation[-c(8,9),]
##reoder the Z-summary
k.order=order((CN.Dutch.preservation[,2]+CN.SAFHS.preservation[,2]),decreasing=T)
CN.Dutch.preservation=CN.Dutch.preservation[k.order,]
CN.SAFHS.preservation=CN.SAFHS.preservation[k.order,]
CN.SAFHS.preservation[,1]=paste(rownames(CN.SAFHS.preservation), " SAFHS",sep="/")
CN.Dutch.preservation[,1]=paste(rownames(CN.Dutch.preservation), "Dutch",sep="/")
colnames(CN.Dutch.preservation)[1]="Module_name"
colnames(CN.SAFHS.preservation)[1]="Module_name"
bine=data.frame(matrix(NA,46,3))
colnames(bine)=c("Module_Colour","Name","Preservation_Z_Score")
for(i in 1:23){ bine[2*i,1]=rownames(CN.SAFHS.preservation)[i];
bine[2*i,2:3]=CN.SAFHS.preservation[i,1:2]}
for(i in 1:23){ bine[2*i-1,1]=rownames(CN.Dutch.preservation)[i];
bine[2*i-1,2:3]=CN.Dutch.preservation[i,1:2]}
write.csv(bine,file="CN_preservation.csv")
| |Module_Colour |Name |Preservation_Z_Score |
|1 |yellow |yellow/Dutch |12.97454994 |
|2 |yellow |yellow/ SAFHS |13.70931579 |
|3 |purple |purple/Dutch |3.145796414 |
|4 |purple |purple/ SAFHS |9.076266574 |
|5 |aliceblue |aliceblue/Dutch |4.212104428 |
|6 |aliceblue |aliceblue/ SAFHS |6.542663275 |
|7 |blue |blue/Dutch |1.773074091 |
|8 |blue |blue/ SAFHS |3.981097565 |
|9 |royalblue |royalblue/Dutch |2.639728838 |
|10 |royalblue |royalblue/ SAFHS |2.604201244 |
|11 |orange |orange/Dutch |1.640152404 |
|12 |orange |orange/ SAFHS |3.570679813 |
|13 |palegreen |palegreen/Dutch |1.932474515 |
|14 |palegreen |palegreen/ SAFHS |3.044461559 |
|15 |powderblue |powderblue/Dutch |1.14310802 |
|16 |powderblue |powderblue/ SAFHS |1.623969169 |
|17 |black |black/Dutch |1.958266918 |
|18 |black |black/ SAFHS |0.374511844 |
|19 |seashell |seashell/Dutch |1.272672518 |
|20 |seashell |seashell/ SAFHS |0.55379216 |
|21 |lavenderblush |lavenderblush/Dutch |0.806347298 |
|22 |lavenderblush |lavenderblush/ SAFHS |0.995589803 |
|23 |aquamarine |aquamarine/Dutch |0.245782268 |
|24 |aquamarine |aquamarine/ SAFHS |1.057314703 |
|25 |mistyrose |mistyrose/Dutch |0.207075493 |
|26 |mistyrose |mistyrose/ SAFHS |0.607545362 |
|27 |chartreuse |chartreuse/Dutch |1.037967391 |
|28 |chartreuse |chartreuse/ SAFHS |-0.30829781 |
|29 |turquoise |turquoise/Dutch |1.164081625 |
|30 |turquoise |turquoise/ SAFHS |-0.650261332 |
|31 |darkgreen |darkgreen/Dutch |-0.48576922 |
|32 |darkgreen |darkgreen/ SAFHS |0.716000313 |
|33 |palevioletred |palevioletred/Dutch |-0.673229763 |
|34 |palevioletred |palevioletred/ SAFHS |0.732141535 |
|35 |brown |brown/Dutch |-0.065184902 |
|36 |brown |brown/ SAFHS |0.047907077 |
|37 |mintcream |mintcream/Dutch |-0.48280317 |
|38 |mintcream |mintcream/ SAFHS |-0.128178728 |
|39 |maroon |maroon/Dutch |-0.545291746 |
|40 |maroon |maroon/ SAFHS |-0.346729284 |
|41 |sandybrown |sandybrown/Dutch |-0.468715399 |
|42 |sandybrown |sandybrown/ SAFHS |-0.545966968 |
|43 |salmon |salmon/Dutch |-1.136439923 |
|44 |salmon |salmon/ SAFHS |-1.056436118 |
|45 |red |red/Dutch |-2.44000899 |
|46 |red |red/ SAFHS |-2.985193505 |
#extracts the preservation result for CB modules
#load the data
rm(list=ls())
load("CB_Dutch_MP_Result.RData")
load("CB_SAFHS_MP_Result.RData")
CB.SAFHS.preservation=mp.CB.SAFHS$preservation$Z$1$2 [,1:2]
CB.Dutch.preservation=mp.CB.Dutch$preservation$Z$1$ 2 [,1:2]
##remove grey and gold module which is useless
CB.Dutch.preservation=CB.Dutch.preservation[-c(10,12),]
CB.SAFHS.preservation=CB.SAFHS.preservation[-c(10,12),]
##reoder the Z-summary
k.order=order((CB.Dutch.preservation[,2]+CB.SAFHS.preservation[,2]),decreasing=T)
CB.Dutch.preservation=CB.Dutch.preservation[k.order,]
CB.SAFHS.preservation=CB.SAFHS.preservation[k.order,]
CB.SAFHS.preservation[,1]=paste(rownames(CB.SAFHS.preservation), "SAFHS",sep="/")
CB.Dutch.preservation[,1]=paste(rownames(CB.Dutch.preservation), "Dutch",sep="/")
colnames(CB.Dutch.preservation)[1]="Module_name"
colnames(CB.SAFHS.preservation)[1]="Module_name"
bine=data.frame(matrix(NA,44,3))
colnames(bine)=c("Module_Colour","Name","Preservation_Z_Score")
for(i in 1:22){ bine[2*i-1,1]=rownames(CB.SAFHS.preservation)[i];
bine[2*i-1,2:3]=CB.SAFHS.preservation[i,1:2]}
for(i in 1:22){ bine[2*i,1]=rownames(CB.Dutch.preservation)[i];
bine[2*i,2:3]=CB.Dutch.preservation[i,1:2]}
write.csv(bine,file=”CB_preservation.csv")
| |Module_Colour |Name |Preservation_Z_Score |
|1 |blue |blue/SAFHS |10.4069478 |
|2 |blue |blue/Dutch |9.06166026 |
|3 |yellow |yellow/SAFHS |3.841849398 |
|4 |yellow |yellow/Dutch |4.346001669 |
|5 |green |green/SAFHS |2.430903649 |
|6 |green |green/Dutch |2.71025349 |
|7 |palegreen |palegreen/SAFHS |3.576070651 |
|8 |palegreen |palegreen/Dutch |0.910745184 |
|9 |purple |purple/SAFHS |2.900397986 |
|10 |purple |purple/Dutch |0.954071119 |
|11 |black |black/SAFHS |2.152411481 |
|12 |black |black/Dutch |1.496312466 |
|13 |pink |pink/SAFHS |1.480209711 |
|14 |pink |pink/Dutch |1.68689298 |
|15 |darkgoldenrod |darkgoldenrod/SAFHS |2.004893963 |
|16 |darkgoldenrod |darkgoldenrod/Dutch |1.101433275 |
|17 |darksalmon |darksalmon/SAFHS |1.156483961 |
|18 |darksalmon |darksalmon/Dutch |1.776318273 |
|19 |turquoise |turquoise/SAFHS |1.320163345 |
|20 |turquoise |turquoise/Dutch |0.825920875 |
|21 |tan |tan/SAFHS |2.94299485 |
|22 |tan |tan/Dutch |-0.900185887 |
|23 |chocolate |chocolate/SAFHS |1.074977345 |
|24 |chocolate |chocolate/Dutch |0.827483111 |
|25 |honeydew |honeydew/SAFHS |0.531819287 |
|26 |honeydew |honeydew/Dutch |0.686770023 |
|27 |cornflowerblue |cornflowerblue/SAFHS |1.569768607 |
|28 |cornflowerblue |cornflowerblue/Dutch |-0.608223848 |
|29 |violetred |violetred/SAFHS |-0.567417797 |
|30 |violetred |violetred/Dutch |1.4341575 |
|31 |darkblue |darkblue/SAFHS |0.856480417 |
|32 |darkblue |darkblue/Dutch |-0.837777913 |
|33 |lemonchiffon |lemonchiffon/SAFHS |-0.248973259 |
|34 |lemonchiffon |lemonchiffon/Dutch |0.120858542 |
|35 |khaki |khaki/SAFHS |-0.131740627 |
|36 |khaki |khaki/Dutch |-0.088070589 |
|37 |wheat |wheat/SAFHS |-0.659923353 |
|38 |wheat |wheat/Dutch |-0.0683524 |
|39 |darkolivegreen |darkolivegreen/SAFHS |-0.240828513 |
|40 |darkolivegreen |darkolivegreen/Dutch |-1.037343996 |
|41 |brown |brown/SAFHS |-0.880146023 |
|42 |brown |brown/Dutch |-1.586842766 |
|43 |red |red/SAFHS |-2.362060973 |
|44 |red |red/Dutch |-2.381356472 |
#plot the result
source("barplot_code.txt")
pdf(10,7.5,file="Figure1_CTX_CN_CB.pdf")
layout(matrix(c(0,1,2,3),1,4,byrow=TRUE), c(0.5,1,1,1), c(2.5), TRUE)
## 1 Plot CTX data
data.CTX=read.csv("CTX_preservation.csv")
data.CTX=data.CTX[,-1]
data.CTX [,1]=as.character(data.CTX [,1])
data.CTX [,2]=as.character(data.CTX [,2])
k=1:38
k=39-k
data.CTX = data.CTX [k,]
## read the color
colors.CTX=c(rep(NA,8),as.character(data.CTX [,1]))
## read the Z-summary
x.CTX= c(rep(NA,8),data.CTX [,3])
##read the name
xlimits=c(0,20)
title=" Preservation of CTX modules"
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x.CTX, colors.CTX, axes = F, main=title, names=c(rep(NA,8),data.CTX[,2]), xlim=xlimits, cex.names=0.7, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.05 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.25*(box[2]-box[1]),box[4]+0.035*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 2)
## 2 Plot CN data
=read.csv("CN_preservation.csv")
=[,-1]
[,1]=as.character( [,1])
[,2]=as.character( [,2])
k=1:46
k=47-k
= [k,]
## read the color
=as.character([,1])
## read the Z-summary
=[,3]
##read the name
xlimits=c(0,20)
title=" Preservation of CN modules"
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(, , axes = F, axisnames = T, main=title, names=[,2], xlim=xlimits, cex.names=0.7, family = "sans", xlab ="red")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.05 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.25*(box[2]-box[1]),box[4]+0.035*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 2)
## 3 Plot CB data
data.CB=read.csv("CB_preservation.csv")
data.CB=data.CB[,-1]
data.CB[,1]=as.character(data.CB [,1])
data.CB [,2]=as.character(data.CB [,2])
k=1:44
k=45-k
data.CB = data.CB [k,]
## read the color
colors.CB=c(rep(NA,2),as.character(data.CB[,1]))
## read the Z-summary
x.CB=c(rep(NA,2),data.CB[,3])
##read the name
xlimits=c(0,20)
title=" Preservation of CB modules"
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x.CB, colors.CB, axes = FALSE, main=title, names=c(rep(NA,2),data.CB[,2]), xlim=xlimits, cex.names=0.7, family = "sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.05 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.25*(box[2]-box[1]),box[4]+0.035*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 2)
dev.off()
#Module relationship
#Module eigengenes are good represents for modules, and we will compare the relationship between modules with eigengenes.
#load WGCNA
library(WGCNA)
#calculate the eigengenes for 3 preserved CTX modules
#load CTX related data
load("CTX_SAFHS_Dutch_for_MP.RData")
#calculate the eigenges for all CTX modules in SAFHS and Dutch data
ME.SAFHS.CTX=moduleEigengenes(t(data.SAFHS),color.CTX)$eigengenes
ME.Dutch.CTX=moduleEigengenes(t(data.Dutch),color.CTX)$eigengenes
#extract the eigengenes for three preserved modules
ME.blue.CTX.Dutch=ME.Dutch.CTX$MEblue
ME.yellow.CTX.Dutch=ME.Dutch.CTX$MEyellow
ME.green.CTX.Dutch=ME.Dutch.CTX$MEgreen
ME.blue.CTX.SAFHS=ME.SAFHS.CTX$MEblue
ME.yellow.CTX.SAFHS=ME.SAFHS.CTX$MEyellow
ME.green.CTX.SAFHS=ME.SAFHS.CTX$MEgreen
#load CN related data
load("CN_SAFHS_Dutch_for_MP.RData")
#calculate the eigenges for all CN modules in SAFHS and Dutch data
ME.=moduleEigengenes(t(data.SAFHS),)$eigengenes
ME.=moduleEigengenes(t(data.Dutch),)$eigengenes
#extract the eigengenes preserved module
ME..Dutch=ME.$MEyellow
ME..SAFHS=ME.$MEyellow
#load CB related data
load("CB_SAFHS_Dutch_for_MP.RData")
#calculate the eigenges for all CN modules in SAFHS and Dutch data
ME.SAFHS.CB=moduleEigengenes(t(data.SAFHS),color.CB)$eigengenes
ME.Dutch.CB=moduleEigengenes(t(data.Dutch),color.CB)$eigengenes
#extract the eigengenes preserved module
ME.blue.CB.Dutch=ME.Dutch.CB$MEblue
ME.blue.CB.SAFHS=ME.SAFHS.CB$MEblue
#put the eigengenes from Dutch data togeather
ME.Dutch=data.frame(
ME.green.CTX.Dutch,
ME.blue.CTX.Dutch,
ME..Dutch,
ME.yellow.CTX.Dutch,
ME.blue.CB.Dutch)
#put the eigengenes from SAFHS data togeather
ME.SAFHS=data.frame(
ME.green.CTX.SAFHS,
ME.yellow.CTX.SAFHS,
ME.blue.CB.SAFHS,
ME.blue.CTX.SAFHS,
ME..SAFHS
)
#plot the relationship of 5 preserved modules
pdf(8,4.2,file="Figure2_Module_relationship.pdf")
layout(matrix(c(0,1,0,2,3,3,4,4),2,4,byrow=TRUE), c(0.5,3.5,0.5,3.5), c(1.2,2.2), TRUE)
par(mar = c(0, 5.2, 1,6))
plot(flashClust(as.dist(1-cor(ME.Dutch,use="p")),method="average"), main="Dutch",sub="",xlab="" ,ylab ="Height", ylim = c(0,1),labels=c("","","","",""))
box = par("usr");
text(box[1]-0.24*(box[2]-box[1]),box[4]+0.0*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 2)
plot(flashClust(as.dist(1-cor(ME.SAFHS,use="p")),method="average"), main="SAFHS",sub="",xlab="" ,ylab ="Height", ylim = c(0,1),labels=c("","","","",""))
box = par("usr");
text(box[1]-0.24*(box[2]-box[1]),box[4]+0.0*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 2)
par(mar = c(5, 8, 0,3))
pME = corPvalueFisher(cor(ME.Dutch),dim(ME.Dutch)[1])
labeledHeatmap(
cor(ME.SAFHS,use="p"),
textMat = matrix(paste(signif(cor(ME.Dutch,use="p"), 2),"\n",signif(pME,1)),5,5),
xLabels=c("ME(green/CTX)",
"ME(blue/CTX)",
"ME(yellow/CN)",
"ME(yellow/CTX)",
"ME(blue/CB)"),
yLabels=c("ME(green/CTX)",
"ME(blue/CTX)",
"ME(yellow/CN)",
"ME(yellow/CTX)",
"ME(blue/CB)"),
colors = greenWhiteRed(50),
setStdMargins = F
)
pME = corPvalueFisher(cor(ME.SAFHS),dim(ME.SAFHS)[1])
labeledHeatmap(
cor(ME.SAFHS,use="p"),
textMat = matrix(paste(signif(cor(ME.SAFHS,use="p"), 2),"\n",signif(pME,1)),5,5),
xLabels=c("ME(green/CTX)",
"ME(blue/CTX)",
"ME(yellow/CN)",
"ME(yellow/CTX)",
"ME(blue/CB)"),
yLabels=c("ME(green/CTX)",
"ME(blue/CTX)",
"ME(yellow/CN)",
"ME(yellow/CTX)",
"ME(blue/CB)"),
colors = greenWhiteRed(50),
setStdMargins =FALSE
)
dev.off()
#Module membership analysis
#calculate the module membership
#load WGCNA package.
library(WGCNA)
#load the data
load("CTX_SAFHS_Dutch_for_MP.RData")
#get the module eigengenes
ME.CTX=moduleEigengenes(t(data.CTX),color.CTX)$eigengenes
ME.Dutch=moduleEigengenes(t(data.Dutch),color.CTX)$eigengenes
ME.SAFHS=moduleEigengenes(t(data.SAFHS),color.CTX)$eigengenes
#calculate the module membership
kME.CTX=signedKME(t(data.CTX),ME.CTX)
kME.Dutch=signedKME(t(data.Dutch),ME.Dutch)
kME.SAFHS=signedKME(t(data.SAFHS),ME.SAFHS)
#extract the blue green and yellow modules’ module membership
KME.ctx=data.frame(color.CTX,kME.CTX$kMEyellow,kME.Dutch$kMEyellow,kME.SAFHS$kMEyellow,kME.CTX$kMEblue,kME.Dutch$kMEblue,kME.SAFHS$kMEblue,kME.CTX$kMEgreen,kME.Dutch$kMEgreen,kME.SAFHS$kMEgreen)
rownames(KME.ctx)=rownames(data.CTX)
write.csv(KME.ctx,file="kME_CTX_yellow_blue_green.csv")
rm(list=ls())
#load the data
load("CN_SAFHS_Dutch_for_MP.RData")
#get the module eigengenes
=moduleEigengenes(t(),)$eigengenes
ME.Dutch=moduleEigengenes(t(data.Dutch),)$eigengenes
ME.SAFHS=moduleEigengenes(t(data.SAFHS),)$eigengenes
#calculate the module membership
=signedKME(t(),)
kME.Dutch=signedKME(t(data.Dutch),ME.Dutch)
kME.SAFHS=signedKME(t(data.SAFHS),ME.SAFHS)
=data.frame(,$kMEyellow,kME.Dutch$kMEyellow,kME.SAFHS$kMEyellow)
rownames()=rownames()
write.csv(,file="kME_CN_yellow.csv")
rm(list=ls())
#load the data
load("CB_SAFHS_Dutch_for_MP.RData")
#get the module eigengenes
ME.CB=moduleEigengenes(t(data.CB),color.CB)$eigengenes
ME.Dutch=moduleEigengenes(t(data.Dutch),color.CB)$eigengenes
ME.SAFHS=moduleEigengenes(t(data.SAFHS),color.CB)$eigengenes
#calculate the module membership
kME.CB=signedKME(t(data.CB),ME.CB)
kME.Dutch=signedKME(t(data.Dutch),ME.Dutch)
kME.SAFHS=signedKME(t(data.SAFHS),ME.SAFHS)
KME.cb=data.frame(color.CB,kME.CB$kMEblue,kME.Dutch$kMEblue,kME.SAFHS$kMEblue)
rownames(KME.cb)=rownames(data.CB)
write.csv(KME.cb,file="kME_CB_blue.csv")
#analysis base on module membership
#load the data from csv files produced by previous steps
rm(list=ls())
data.CTX=read.csv("kME_CTX_yellow_blue_green.csv")
names(data.CTX)
names(data.CTX)=c("Gene.CTX", "color.CTX","CTX.yellow", "Dutch.yellow.ctx", "SAFHS.yellow.ctx", "CTX.blue", "Dutch.blue.ctx", "SAFHS.blue.ctx", "CTX.green", "Dutch.green.ctx", "SAFHS.green.ctx")
attach(data.CTX)
=read.csv("kME_CN_yellow.csv")
names()
names()=c("", "", "CN.yellow", "Dutch.", "SAFHS.")
attach()
data.CB=read.csv("kME_CB_blue.csv")
names(data.CB)
names(data.CB)=c("Gene.CB", "color.CB", "CB.blue", "Dutch.blue.cb", "SAFHS.blue.cb")
attach(data.CB)
#explore the relationship of kME between two blood data sets.
pdf(12,8,file="Figure3_kME_blood_Consistance.pdf")
par(mfrow=c(2,3))
par(mar=c(7,8,4,1))
par(mgp = c(4,1,0))
#1
verboseScatterplot(SAFHS.yellow.ctx, Dutch.yellow.ctx,xlab="Module membership\nSAFHS(yellow/CTX)", ylab="Module membership\nDutch(yellow/CTX)", main = "Yellow/CTX module in Dutch and SAFHS\n",cex.main = 1.4, corFnc = "bicor")
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.12*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#2
verboseScatterplot(SAFHS.blue.ctx,Dutch.blue.ctx,xlab="Module membership\nSAFHS(blue/CTX)", ylab="Module membership\nDutch(blue/CTX)",
main = "Blue/CTX module in Dutch and SAFHS\n",cex.main = 1.4, corFnc = "bicor")
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.12*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#3
verboseScatterplot(SAFHS.green.ctx,Dutch.green.ctx,xlab="Module membership\nSAFHS(green/CTX)", ylab="Module membership\nDutch(green/CTX)",
main = "Green/CTX module in Dutch and SAFHS\n",cex.main = 1.4, corFnc = "bicor")
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.12*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#4
verboseScatterplot(SAFHS., Dutch.,xlab="Module membership\nSAFHS(yellow/CN)", ylab="Module membership\nDutch(yellow/CN)",
main = "Yellow/CN module in Dutch and SAFHS\n",cex.main = 1.4, corFnc = "bicor")
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.12*(box[4]-box[3]), "d", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#5
verboseScatterplot(SAFHS.blue.cb,Dutch.blue.cb,xlab="Module membership\nSAFHS(blue/CB)", ylab="Module membership\nDutch(blue/CB)",
main = "Blue/CB module in Dutch and SAFHS\n",cex.main = 1.4, corFnc = "bicor")
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.12*(box[4]-box[3]), "e", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
#6
yellow_CTX=0.76
blue_CTX=0.74
green_CTX=0.76
yellow_CN=0.73
blue_CB=0.74
barplot(c(yellow_CTX,blue_CTX,green_CTX,yellow_CN,blue_CB),col=c("yellow","blue","green","yellow","blue"),
main = "Summary of MM correlation for preserved\n modules between two blood data sets",
xlab = "Modules",ylab = "Value of Bicorrelation",space =0.9,cex.names = 0.9,names.arg=c("yellow/CTX","blue/CTX","green/CTX","yellow/CN","blue/CB"),
ylim = c(0, 0.9), cex.main=1.3,cex.lab=1.5,las=1)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.12*(box[4]-box[3]), "f", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
dev.off()
#since the kME for each module shows significantly high correlation, we try to combine them in #each module by calibrating the kME values we fit linear models
lmbluectx=lm(SAFHS.blue.ctx~ Dutch.blue.ctx,na.action="na.exclude")
lmgreenctx=lm(SAFHS.green.ctx~ Dutch.green.ctx,na.action="na.exclude")
lmyellowctx=lm(SAFHS.yellow.ctx~ Dutch.yellow.ctx,na.action="na.exclude")
lmbluecb=lm(SAFHS.blue.cb~ Dutch.blue.cb,na.action="na.exclude")
lmyellowcn=lm(SAFHS.~ Dutch.,na.action="na.exclude")
summary(lmbluectx)
Call:
lm(formula = SAFHS.blue.ctx ~ Dutch.blue.ctx, na.action = "na.exclude")
Residuals:
Min 1Q Median 3Q Max
-0.930018 -0.154026 -0.003733 0.153460 1.006253
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.021707 0.004842 4.483 7.69e-06 ***
Dutch.blue.ctx 0.653613 0.011820 55.299 < 2e-16 ***
---
Signif. codes: 0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1
Residual standard error: 0.2406 on 2638 degrees of freedom
Multiple R-squared: 0.5369, Adjusted R-squared: 0.5367
F-statistic: 3058 on 1 and 2638 DF, p-value: < 2.2e-16
summary(lmgreenctx)
Call:
lm(formula = SAFHS.green.ctx ~ Dutch.green.ctx, na.action = "na.exclude")
Residuals:
Min 1Q Median 3Q Max
-0.999302 -0.149303 0.001807 0.151326 0.826849
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.008301 0.004661 -1.781 0.075 .
Dutch.green.ctx 0.667939 0.011403 58.575 |t|)
(Intercept) 0.018606 0.004723 3.939 8.38e-05 ***
Dutch.yellow.ctx 0.673635 0.011511 58.521 < 2e-16 ***
---
Signif. codes: 0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1
Residual standard error: 0.236 on 2638 degrees of freedom
Multiple R-squared: 0.5649, Adjusted R-squared: 0.5647
F-statistic: 3425 on 1 and 2638 DF, p-value: < 2.2e-16
summary(lmbluecb)
Call:
lm(formula = SAFHS.blue.cb ~ Dutch.blue.cb, na.action = "na.exclude")
Residuals:
Min 1Q Median 3Q Max
-0.86727 -0.14717 -0.01038 0.15276 1.00952
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.013953 0.005499 2.537 0.0112 *
Dutch.blue.cb 0.653565 0.013326 49.043 |t|)
(Intercept) 0.026176 0.005466 4.789 1.8e-06 ***
Dutch. 0.641702 0.013613 47.139 < 2e-16 ***
---
Signif. codes: 0 â***â 0.001 â**â 0.01 â*â 0.05 â.â 0.1 â â 1
Residual standard error: 0.2378 on 2061 degrees of freedom
Multiple R-squared: 0.5188, Adjusted R-squared: 0.5186
F-statistic: 2222 on 1 and 2061 DF, p-value: < 2.2e-16
#This suggests to define the following average value:
kMEbloodbluectx= 1/1.66*SAFHS.blue.ctx+.66/1.66*Dutch.blue.ctx
kMEbloodgreenctx=1/1.66*SAFHS.green.ctx+.66/1.66*Dutch.green.ctx
kMEbloodyellowctx= 1/1.66*SAFHS.yellow.ctx+.66/1.66*Dutch.yellow.ctx
kMEbloodbluecb= 1/1.66*SAFHS.blue.cb+.66/1.66*Dutch.blue.cb
kMEbloodyellowcn=1/1.66*SAFHS.+.66/1.66*Dutch.
datBluectx=data.frame(kMEbloodbluectx , SAFHS.blue.ctx, Dutch.blue.ctx)
datGreenctx=data.frame(kMEbloodgreenctx, SAFHS.green.ctx, Dutch.green.ctx)
datYellowctx=data.frame(kMEbloodyellowctx, SAFHS.yellow.ctx, Dutch.yellow.ctx)
datBluecb=data.frame(kMEbloodbluecb, SAFHS.blue.cb, Dutch.blue.cb)
datYellowcn=data.frame(kMEbloodyellowcn, SAFHS., Dutch.)
#how the combined kME represent the old one
signif(cor(datBluectx, use="p"),2)
kMEbloodbluectx SAFHS.blue.ctx Dutch.blue.ctx
kMEbloodbluectx 1.00 0.95 0.91
SAFHS.blue.ctx 0.95 1.00 0.73
Dutch.blue.ctx 0.91 0.73 1.00
signif(cor(datGreenctx, use="p"),2)
kMEbloodgreenctx SAFHS.green.ctx Dutch.green.ctx
kMEbloodgreenctx 1.00 0.95 0.91
SAFHS.green.ctx 0.95 1.00 0.75
Dutch.green.ctx 0.91 0.75 1.00
signif(cor(datYellowctx, use="p"),2)
kMEbloodyellowctx SAFHS.yellow.ctx Dutch.yellow.ctx
kMEbloodyellowctx 1.00 0.95 0.91
SAFHS.yellow.ctx 0.95 1.00 0.75
Dutch.yellow.ctx 0.91 0.75 1.00
signif(cor(datBluecb, use="p"),2)
kMEbloodbluecb SAFHS.blue.cb Dutch.blue.cb
kMEbloodbluecb 1.00 0.95 0.91
SAFHS.blue.cb 0.95 1.00 0.74
Dutch.blue.cb 0.91 0.74 1.00
signif(cor(datYellowcn, use="p"),2)
kMEbloodyellowcn SAFHS. Dutch.
kMEbloodyellowcn 1.00 0.95 0.90
SAFHS. 0.95 1.00 0.72
Dutch. 0.90 0.72 1.00
#explore the relationship of kME between three modules.
pdf(12,8,file="Figure4_kME_Relationship_3_modules.pdf")
par(mfrow=c(2,3))
par(mar=c(7,8,4,1))
par(mgp = c(4,1,0))
verboseScatterplot(kMEbloodyellowctx, kMEbloodbluectx, xlab="Module membership\nBlood(yellow/CTX)", ylab="Module membership\nBlood(blue/CTX)",
main = "Blue/CTX versus Yellow/CTX in Blood\n",cex.main = 1.4)
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.14*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot(kMEbloodyellowctx, kMEbloodgreenctx, xlab="Module membership\nBlood(yellow/CTX)", ylab="Module membership\nBlood(green/CTX)",
main = "Green/CTX versus Yellow/CTX in Blood\n",cex.main = 1.4)
abline(0,-1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.14*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot(kMEbloodbluectx, kMEbloodgreenctx, xlab="Module membership\nBlood(blue/CTX)", ylab="Module membership\nBlood(green/CTX)",
main = "Green/CTX versus Blue/CTX in Blood\n",cex.main = 1.4)
abline(0,-1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.14*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot(CTX.yellow, CTX.blue, xlab="Module membership\nCTX(yellow)", ylab="Module membership\nCTX(blue)",
main = "Blue/CTX versus Yellow/CTX\n",cex.main = 1.4)
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.14*(box[4]-box[3]), "d", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot(CTX.yellow,CTX.green, xlab="Module membership\nCTX(yellow)", ylab="Module membership\nCTX(green)",
main = "Green/CTX versus Yellow/CTX\n",cex.main = 1.4)
abline(0,-1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.15*(box[4]-box[3]), "e", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot(CTX.blue,CTX.green, xlab="Module membership\nCTX(blue)", ylab="Module membership\nCTX(green)", main = "Green/CTX versus Blue/CTX\n",cex.main = 1.4)
abline(0,-1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.14*(box[4]-box[3]), "f", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
dev.off()
#we also combine the kME from 3 preserved CTX modules base on the previous plot
kMEBloodAveragectx= (kMEbloodbluectx+kMEbloodyellowctx-kMEbloodgreenctx)/3
kMEBloodAveragecn= kMEbloodyellowcn
kMEBloodAveragecb= kMEbloodbluecb
kMECTXAverage= (CTX.blue+ CTX.yellow- CTX.green)/3
kMECNAverage= CN.yellow
kMECBAverage= CB.blue
#find how well the new kME can represent the old one
cor(data.frame(kMEBloodAveragectx , kMEbloodbluectx, kMEbloodyellowctx, kMEbloodgreenctx))
kMEBloodAveragectx kMEbloodbluectx kMEbloodyellowctx
kMEBloodAveragectx 1.0000000 0.9974751 0.9993301
kMEbloodbluectx 0.9974751 1.0000000 0.9960215
kMEbloodyellowctx 0.9993301 0.9960215 1.0000000
kMEbloodgreenctx -0.9976269 -0.9909014 -0.9963894
kMEbloodgreenctx
kMEBloodAveragectx -0.9976269
kMEbloodbluectx -0.9909014
kMEbloodyellowctx -0.9963894
kMEbloodgreenctx 1.0000000
cor(data.frame(kMECTXAverage , CTX.blue, CTX.yellow, CTX.green))
kMECTXAverage CTX.blue CTX.yellow CTX.green
kMECTXAverage 1.0000000 0.68593783 0.9697139 -0.73918489
CTX.blue 0.6859378 1.00000000 0.6152591 -0.04479907
CTX.yellow 0.9697139 0.61525910 1.0000000 -0.70790381
CTX.green -0.7391849 -0.04479907 -0.7079038 1.00000000
#try to select the preserved gene base on the threshold |kME|>=0.35
#for CTX
Positive.CTX=data.CTX[kMEBloodAveragectx >.35& kMECTXAverage >.35,]
Negative.CTX=data.CTX[-kMEBloodAveragectx >.35& -kMECTXAverage >.35,]
Genelist.CTX=c(as.character(Positive.CTX[,1]),as.character(Negative.CTX[,1]))
write.table(Genelist.CTX,file="Genelist_CTX.txt", col.names=F,row.names=F)
k.ctx=match(Genelist.CTX,as.character(data.CTX[,1]))
Genelistall.CTX=as.character(data.CTX[,1])
colorctx=as.character(color.CTX)
colorctx[colorctx!="black"]="black"
colorctx[k.ctx]= "red"
#for CN
=[kMEBloodAveragecn >.35& kMECNAverage >.35,]
=[-kMEBloodAveragecn >.35& -kMECNAverage >.35,]
=c(as.character([,1]),as.character([,1]))
write.table(,file="Genelist_CN.txt", col.names=F,row.names=F)
=match(,as.character([,1]))
=as.character([,1])
colorcn=as.character()
colorcn[colorcn!="black"]="black"
colorcn[]= "yellow"
#for CB
Positive.CB=data.CB[kMEBloodAveragecb >.35& kMECBAverage >.35,]
Negative.CB=data.CB[-kMEBloodAveragecb >.35& -kMECBAverage >.35,]
Genelist.CB=c(as.character(Positive.CB[,1]),as.character(Negative.CB[,1]))
write.table(Genelist.CB,file="Genelist_CB.txt", col.names=F,row.names=F)
k.cb=match(Genelist.CB,as.character(data.CB[,1]))
Genelistall.CB=as.character(data.CB[,1])
colorcb=as.character(color.CB)
colorcb[colorcb!="black"]="black"
colorcb[k.cb]= "blue"
#Heritability analysis
#load the heritability data for every gene
Heritability=read.csv("ng2119-S2.csv")[2:3]
Heritability=apply(Heritability[,c(2,2)],2,tapply,Heritability[,1],mean)
Heritability=data.frame(rownames(Heritability),Heritability[,1])
colnames(Heritability)=c("Gene_Symbol","Heritability")
load("SAFHS_datExpr.RData")
all.gene=rownames(Heritability)
k=match(all.gene,rownames(datExpr.SAFHS))
k=k[-(1:length(k))[is.na(k)]]
datExpr.SAFHS=datExpr.SAFHS[k,]
datExpr.SAFHS.expression=data.frame(rownames(datExpr.SAFHS),apply(datExpr.SAFHS,1,mean))
All=Heritability[,2]
CTX=Heritability[match(Genelist.CTX,Heritability[,1]),2]
CTX=CTX[-(1:length(CTX))[is.na(CTX)]]
CTX.all=Heritability[match(Genelistall.CTX,Heritability[,1]),2]
CTX.all=CTX.all[-(1:length(CTX.all))[is.na(CTX.all)]]
CN=Heritability[match(,Heritability[,1]),2]
CN=CN[-(1:length(CN))[is.na(CN)]]
CN.all=Heritability[match(,Heritability[,1]),2]
CN.all=CN.all[-(1:length(CN.all))[is.na(CN.all)]]
CB=Heritability[match(Genelist.CB,Heritability[,1]),2]
CB=CB[-(1:length(CB))[is.na(CB)]]
CB.all=Heritability[match(Genelistall.CB,Heritability[,1]),2]
CB.all=CB.all[-(1:length(CB.all))[is.na(CB.all)]]
datExpr.SAFHS.expression[,1]=as.character(datExpr.SAFHS.expression[,1])
All.1=datExpr.SAFHS.expression[,2]
CTX.1=datExpr.SAFHS.expression[match(Genelist.CTX,datExpr.SAFHS.expression[,1]),2]
CTX.1=CTX.1[-(1:length(CTX.1))[is.na(CTX.1)]]
CTX.all.1=datExpr.SAFHS.expression[match(Genelistall.CTX,datExpr.SAFHS.expression[,1]),2]
CTX.all.1=CTX.all.1[-(1:length(CTX.all.1))[is.na(CTX.all.1)]]
CN.1=datExpr.SAFHS.expression[match(,datExpr.SAFHS.expression[,1]),2]
CN.1=CN.1[-(1:length(CN.1))[is.na(CN.1)]]
CN.all.1=datExpr.SAFHS.expression[match(,datExpr.SAFHS.expression[,1]),2]
CN.all.1=CN.all.1[-(1:length(CN.all.1))[is.na(CN.all.1)]]
CB.1=datExpr.SAFHS.expression[match(Genelist.CB,datExpr.SAFHS.expression[,1]),2]
CB.1=CB.1[-(1:length(CB.1))[is.na(CB.1)]]
CB.all.1=datExpr.SAFHS.expression[match(Genelistall.CB,datExpr.SAFHS.expression[,1]),2]
CB.all.1=CB.all.1[-(1:length(CB.all.1))[is.na(CB.all.1)]]
#We get All, CTX, CN, CB, CTX.all, CN.all, CB.all, and plot the result
pdf(12,12,file="Figure5_Preserved_genes_heritability_Mean_expression.pdf")
par(mfrow=c(3,3))
par(mar=c(7,8,4,1))
par(mgp = c(4,1,0))
verboseScatterplot( kMECTXAverage, kMEBloodAveragectx, ylab="Module membership\nCombined Blood", xlab="Module membership\nCombined CTX",col= colorctx)
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.1*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot( kMECNAverage, kMEBloodAveragecn ,ylab="Module membership\nBlood", xlab="Module membership\nCN",col= colorcn)
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.1*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseScatterplot(kMECBAverage, kMEBloodAveragecb, ylab="Module membership\nBlood", xlab="Module membership\nCB",col= colorcb)
abline(0,1,col="red",lwd=2)
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.1*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
par(mgp = c(4,2,0))
par(mar=c(6,8,4,1))
verboseBarplot(ylim=c(0,0.35),c(All,CTX.all,CTX),c(rep("(n=18525)\nAll genes with\nheritability value",length(All)),rep("(n=2640)\nAll genes in\nCTX network", length(CTX.all)),rep("(n=357)\nHub genes\nin CTX network",length(CTX))), AnovaTest=T, main="CTX",xlab="Group",ylab="Mean Heritability", cex=1, color=c("grey","black","red"))
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.1*(box[4]-box[3]), "d", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseBarplot(ylim=c(0,0.35),c(All,CN.all,CN), c(rep("(n=18525)\nAll genes with\nheritability value",length(All)),rep("(n=2063)\nAll genes in\nCN network", length(CN.all)),rep("(n=305)\nHub genes\nin CN network",length(CN))), AnovaTest=T, main="CN",xlab="Group",ylab="Mean Heritability", cex=1, color=c("grey","black","yellow"))
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.1*(box[4]-box[3]), "e", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseBarplot(ylim=c(0,0.35),c(All,CB.all,CB), c(rep("(n=18525)\nAll genes with\nheritability value",length(All)),rep("(n=2001)\nAll genes in\nCB network", length(CB.all)),rep("(n=277)\nHub genes\nin CB network",length(CB))), AnovaTest=T, main="CB",xlab="Group",ylab="Mean Heritability", cex=1, color=c("grey","black","blue"))
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.1*(box[4]-box[3]), "f", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
ylimits=c(0,400)
verboseBarplot(ylim=c(0,510),c(All.1,CTX.all.1, CTX.1), c(rep("(n=18525)\nAll genes with\nheritability value",length(All.1)),rep("(n=2640)\nAll genes in\nCTX network", length(CTX.all.1)),rep("(n=357)\nHub genes\nin CTX network",length(CTX.1))),
AnovaTest=T, main="CTX",xlab="Group",ylab="Mean Expression", cex=1, color=c("grey","black","red"))
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.11*(box[4]-box[3]), "g", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseBarplot(ylim=c(0,510),c(All.1,CN.all.1,CN.1), c(rep("(n=18525)\nAll genes with\nheritability value",length(All.1)),rep("(n=2063)\nAll genes in\nCN network", length(CN.all.1)),rep("(n=305)\nHub genes\nin CN network",length(CN.1))),
AnovaTest=T, main="CN",xlab="Group",ylab="Mean Expression", cex=1, color=c("grey","black","yellow"))
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.11*(box[4]-box[3]), "h", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
verboseBarplot(ylim=c(0,510),c(All.1,CB.all.1,CB.1), c(rep("(n=18525)\nAll genes with\nheritability value",length(All.1)),rep("(n=2001)\nAll genes in\nCB network", length(CB.all.1)),rep("(n=277)\nHub genes\nin CB network",length(CB.1))),
AnovaTest=T, main="CB",xlab="Group",ylab="Mean Expression", cex=1, color=c("grey","black","blue"))
box = par("usr");
text(box[1]-0.1*(box[2]-box[1]),box[4]+0.11*(box[4]-box[3]), "i", adj = c(0.5, 0.5), xpd = TRUE, cex = 3)
dev.off()
[pic]
#Relate the CD genes with preserved module eigengene
#load WGCNA
library(WGCNA)
#For CTX data
#load the data
load("CTX_SAFHS_Dutch_for_MP.RData")
ME.SAFHS=moduleEigengenes(t(data.SAFHS),color.CTX)$eigengenes
ME.Dutch=moduleEigengenes(t(data.Dutch),color.CTX)$eigengenes
#calculate the combined kME for blue yellow and green
byg.SAFHS=(data.frame(ME.SAFHS$MEblue)+data.frame(ME.SAFHS$MEyellow)+data.frame(ME.SAFHS$MEgreen))/3
byg.Dutch=(data.frame(ME.Dutch$MEblue)+data.frame(ME.Dutch$MEyellow)+data.frame(ME.Dutch$MEgreen))/3
# load the expression of CD genes (all the CD genes in blood expression data) in both Dutch and SAFHS data
load("data.CD.RData")
#calculate the p value for correlation for SAFHS data
pvalue.SAFHS=rep(NA,dim(data.CD.SAFHS)[1])
for(i in 1:dim(data.CD.SAFHS)[1])
{
pvalue.SAFHS[i]= cor.test(byg.SAFHS[,1],data.CD.SAFHS[i,])$p.value
}
SAFHS= data.frame(t(cor((byg.SAFHS),t(data.CD.SAFHS))),pvalue.SAFHS)
SAFHS=data.frame(SAFHS,rep(NA,dim(SAFHS)[1]))
colnames(SAFHS)[1]="Correlation_SAFHS_Yellow"
colnames(SAFHS)[2]="P-value"
colnames(SAFHS)[3]= "Module_CTX"
SAFHS=SAFHS[order(as.numeric(abs(SAFHS[,1])),decreasing=T),]
CTX.module=data.frame(rownames(data.CTX),color.CTX)
k.SAFHS=match(rownames(SAFHS), as.character(CTX.module[,1]))
SAFHS$Module_CTX=CTX.module[k.SAFHS,2]
write.csv(SAFHS,file="SAFHS_CD_CTX.csv")
#calculate the p value for correlation for Dutch data
pvalue.Dutch=rep(NA,dim(data.CD.Dutch)[1])
for(i in 1:dim(data.CD.Dutch)[1])
{
pvalue.Dutch[i]= cor.test(byg.Dutch[,1],data.CD.Dutch[i,])$p.value
}
Dutch= data.frame(t(cor((byg.Dutch),t(data.CD.Dutch))),pvalue.Dutch)
Dutch=data.frame(Dutch,rep(NA,dim(Dutch)[1]))
colnames(Dutch)[1]="Correlation_Dutch_Yellow"
colnames(Dutch)[2]="P-value"
colnames(Dutch)[3]= "Module_CTX"
Dutch=Dutch[order(as.numeric(abs(Dutch[,1])),decreasing=T),]
CTX.module=data.frame(rownames(data.CTX),color.CTX)
k.Dutch=match(rownames(Dutch), as.character(CTX.module[,1]))
Dutch$Module_CTX=CTX.module[k.Dutch,2]
write.csv(Dutch,file="Dutch_CD_CTX.csv")
#For CN data
rm(list=ls())
load("CN_SAFHS_Dutch_for_MP.RData")
ME.SAFHS=moduleEigengenes(t(data.SAFHS),)$eigengenes
ME.Dutch=moduleEigengenes(t(data.Dutch),)$eigengenes
y.SAFHS=data.frame(ME.SAFHS$MEyellow)
y.Dutch=data.frame(ME.Dutch$MEyellow)
#load the expression of CD genes in both Dutch and SAFHS data
load("data.CD.RData")
pvalue.SAFHS=rep(NA,dim(data.CD.SAFHS)[1])
for(i in 1:dim(data.CD.SAFHS)[1])
{
pvalue.SAFHS[i]= cor.test(y.SAFHS[,1],data.CD.SAFHS[i,])$p.value
}
SAFHS= data.frame(t(cor((y.SAFHS),t(data.CD.SAFHS))),pvalue.SAFHS)
SAFHS=data.frame(SAFHS,rep(NA,dim(SAFHS)[1]))
colnames(SAFHS)[1]="Correlation_SAFHS_Yellow"
colnames(SAFHS)[2]="P-value"
colnames(SAFHS)[3]= "Module_CN"
SAFHS=SAFHS[order(as.numeric(abs(SAFHS[,1])),decreasing=T),]
CN.module=data.frame(rownames(),)
k.SAFHS=match(rownames(SAFHS), as.character(CN.module[,1]))
SAFHS$Module_CN=CN.module[k.SAFHS,2]
write.csv(SAFHS,file="SAFHS_CD_CN.csv")
pvalue.Dutch=rep(NA,dim(data.CD.Dutch)[1])
for(i in 1:dim(data.CD.Dutch)[1])
{
pvalue.Dutch[i]= cor.test(y.Dutch[,1],data.CD.Dutch[i,])$p.value
}
Dutch= data.frame(t(cor((y.Dutch),t(data.CD.Dutch))),pvalue.Dutch)
Dutch=data.frame(Dutch,rep(NA,dim(Dutch)[1]))
colnames(Dutch)[1]="Correlation_Dutch_Yellow"
colnames(Dutch)[2]="P-value"
colnames(Dutch)[3]= "Module_CN"
Dutch=Dutch[order(as.numeric(abs(Dutch[,1])),decreasing=T),]
CN.module=data.frame(rownames(),)
k.Dutch=match(rownames(Dutch), as.character(CN.module[,1]))
Dutch$Module_CN=CN.module[k.Dutch,2]
write.csv(Dutch,file="Dutch_CD_CN.csv")
#For CB data
load("CB_SAFHS_Dutch_merge.RData")
ME.SAFHS=moduleEigengenes(t(data.SAFHS),color.CB)$eigengenes
ME.Dutch=moduleEigengenes(t(data.Dutch),color.CB)$eigengenes
b.SAFHS=data.frame(ME.SAFHS$MEblue)
b.Dutch=data.frame(ME.Dutch$MEblue)
load("data.CD.RData")
pvalue.SAFHS=rep(NA,dim(data.CD.SAFHS)[1])
for(i in 1:dim(data.CD.SAFHS)[1])
{
pvalue.SAFHS[i]= cor.test(b.SAFHS[,1],data.CD.SAFHS[i,])$p.value
}
SAFHS= data.frame(t(cor((b.SAFHS),t(data.CD.SAFHS))),pvalue.SAFHS)
SAFHS=data.frame(SAFHS,rep(NA,dim(SAFHS)[1]))
colnames(SAFHS)[1]="Correlation_SAFHS_Blue"
colnames(SAFHS)[2]="P-value"
colnames(SAFHS)[3]= "Module_CB"
SAFHS=SAFHS[order(as.numeric(abs(SAFHS[,1])),decreasing=T),]
CB.module=data.frame(rownames(data.CB),color.CB)
k.SAFHS=match(rownames(SAFHS), as.character(CB.module[,1]))
SAFHS$Module_CB=CB.module[k.SAFHS,2]
write.csv(SAFHS,file="SAFHS_CD_CB.csv")
pvalue.Dutch=rep(NA,dim(data.CD.Dutch)[1])
for(i in 1:dim(data.CD.Dutch)[1])
{
pvalue.Dutch[i]= cor.test(b.Dutch[,1],data.CD.Dutch[i,])$p.value
}
Dutch= data.frame(t(cor((b.Dutch),t(data.CD.Dutch))),pvalue.Dutch)
Dutch=data.frame(Dutch,rep(NA,dim(Dutch)[1]))
colnames(Dutch)[1]="Correlation_Dutch_Blue"
colnames(Dutch)[2]="P-value"
colnames(Dutch)[3]= "Module_CB"
Dutch=Dutch[order(as.numeric(abs(Dutch[,1])),decreasing=T),]
CB.module=data.frame(rownames(data.CB),color.CB)
k.Dutch=match(rownames(Dutch), as.character(CB.module[,1]))
Dutch$Module_CB=CB.module[k.Dutch,2]
write.csv(Dutch,file="Dutch_CD_CB.csv")
#Preservatin study between CTX, CN and CB
#load the data
load("CTX_datExpr_module.RData")
load("CN_datExpr_module.RData")
load("CB_datExpr_module.RData")
#compare CTX and CN with CTX module
#Find out the overlap genes between CTX and CN
name.=as.character(merge(rownames(datExpr.CTX),rownames(),by=1,sort=F)[,1])
datExpr.ctx=datExpr.CTX[match(name.,rownames(datExpr.CTX)),]
=[match(name.,rownames()),]
module.ctx=module.CTX[match(name.,rownames(datExpr.CTX))]
=[match(name.,rownames())]
#Module preservation function
library(WGCNA)
Colors1=module.ctx
data1=datExpr.ctx
data2=
dim(data1)
dim(data2)
colorList=list(Colors1)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data1))
multiExpr[[2]]=list(data=t(data2))
names(multiExpr)=names(colorList)=c("net1","net2")
mp=modulePreservation(multiExpr,colorList,networkType="signed")
save(mp,file=”mp..withctxmodule.RData”)
#compare CTX and CN with CN module
Colors1=
data2=datExpr.ctx
data1=
dim(data1)
dim(data2)
colorList=list(Colors1)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data1))
multiExpr[[2]]=list(data=t(data2))
names(multiExpr)=names(colorList)=c("net1","net2")
mp=modulePreservation(multiExpr,colorList,networkType="signed")
save(mp,file=”mp..withcnmodule.RData”)
rm(list=ls())
#compare CTX and CB with CTX module
rm(list=ls())
load("CTX_datExpr_module.RData")
load("CN_datExpr_module.RData")
load("CB_datExpr_module.RData")
name.ctx.cb=as.character(merge(rownames(datExpr.CTX),rownames(datExpr.CB),by=1,sort=F)[,1])
datExpr.ctx=datExpr.CTX[match(name.ctx.cb,rownames(datExpr.CTX)),]
datExpr.cb=datExpr.CB[match(name.ctx.cb,rownames(datExpr.CB)),]
module.ctx=module.CTX[match(name.ctx.cb,rownames(datExpr.CTX))]
module.cb=module.CB[match(name.ctx.cb,rownames(datExpr.CB))]
library(WGCNA)
Colors1=module.ctx
data1=datExpr.ctx
data2=datExpr.cb
dim(data1)
dim(data2)
colorList=list(Colors1)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data1))
multiExpr[[2]]=list(data=t(data2))
names(multiExpr)=names(colorList)=c("net1","net2")
mp=modulePreservation(multiExpr,colorList,networkType="signed")
save(mp,file=”mp.ctx.cb.withctxmodule.RData”)
#compare CTX and CB with CB module
Colors1=module.cb
data2=datExpr.ctx
data1=datExpr.cb
dim(data1)
dim(data2)
colorList=list(Colors1)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data1))
multiExpr[[2]]=list(data=t(data2))
names(multiExpr)=names(colorList)=c("net1","net2")
mp=modulePreservation(multiExpr,colorList,networkType="signed")
save(mp,file=”mp.ctx.cb.withcbmodule.RData”)
rm(list=ls())
#compare CN and CB with CN module
rm(list=ls())
load("CTX_datExpr_module.RData")
load("CN_datExpr_module.RData")
load("CB_datExpr_module.RData")
.cb=as.character(merge(rownames(),rownames(datExpr.CB),by=1,sort=F)[,1])
=[match(.cb,rownames()),]
datExpr.cb=datExpr.CB[match(.cb,rownames(datExpr.CB)),]
=[match(.cb,rownames())]
module.cb=module.CB[match(.cb,rownames(datExpr.CB))]
library(WGCNA)
Colors1=
data1=
data2=datExpr.cb
dim(data1)
dim(data2)
colorList=list(Colors1)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data1))
multiExpr[[2]]=list(data=t(data2))
names(multiExpr)=names(colorList)=c("net1","net2")
mp=modulePreservation(multiExpr,colorList,networkType="signed")
save(mp,file=”.cb.withcnmodule.RData”)
#compare CN and CB with CB module
library(WGCNA)
Colors1=module.cb
data2=
data1=datExpr.cb
dim(data1)
dim(data2)
colorList=list(Colors1)
colorList[[2]]=NA
multiExpr=list()
multiExpr[[1]]=list(data=t(data1))
multiExpr[[2]]=list(data=t(data2))
names(multiExpr)=names(colorList)=c("net1","net2")
mp=modulePreservation(multiExpr,colorList,networkType="signed")
save(mp,file=”.cb.withcbmodule.RData”)
# Then we extract the Z summary and kept as csv files.
#plot the result
source("barplot_code.txt")
pdf(7,12.5,file=”Supplementary_Figure_CTX_CN_CB_preservation.pdf")
layout(matrix(c(0,0,0,0,1,2,0,3,4,0,5,6),4,3,byrow=TRUE), c(1,3,3), c(0.5,4,4,4), TRUE)
## 1 Plot how well CTX modules preserved in CN data
data1=read.csv ("1..with.ctx.module.csv")
data1 [,1]=as.character(data1[,1])
k=1:23
k=24-k
data1 = data1 [k,]
## read the color
colors1 = as.character(data1[,1])
## read the Z-summary
x1= data1[,3]
##read the name
xlimits=c(0,35)
title="Preservation of CTX modules in CN data "
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x1, colors1, axes = F, main=title,cex.main=1, names= colors1 ,xlim=xlimits, cex.names=1, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.1 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.2*(box[2]-box[1]),box[4]+0.07*(box[4]-box[3]), "a", adj = c(0.5, 0.5), xpd = TRUE, cex = 2,)
## 2 Plot how well CTX modules preserved in CB data
data1=read.csv ("2.ctx.cb.with.ctx.module.csv")
data1 [,1]=as.character(data1[,1])
k=1:23
k=24-k
data1 = data1 [k,]
## read the color
colors1 = as.character(data1[,1])
## read the Z-summary
x1= data1[,3]
##read the name
xlimits=c(0,35)
title="Preservation of CTX modules in CB data "
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x1, colors1, axes = F, main=title,cex.main=1, names= colors1 ,xlim=xlimits, cex.names=1, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.1 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.2*(box[2]-box[1]),box[4]+0.07*(box[4]-box[3]), "b", adj = c(0.5, 0.5), xpd = TRUE, cex = 2,)
## 3 Plot how well CNmodules preserved in CTX data
data1=read.csv ("3...module.csv")
data1 [,1]=as.character(data1[,1])
k=1:23
k=24-k
data1 = data1 [k,]
## read the color
colors1 = as.character(data1[,1])
## read the Z-summary
x1= data1[,3]
##read the name
xlimits=c(0,35)
title="Preservation of CN modules in CTX data "
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x1, colors1, axes = F, main=title,cex.main=1, names= colors1 ,xlim=xlimits, cex.names=1, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.1 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.2*(box[2]-box[1]),box[4]+0.07*(box[4]-box[3]), "c", adj = c(0.5, 0.5), xpd = TRUE, cex = 2,)
## 4 Plot how well CN modules preserved in CB data
data1=read.csv (".cb..module.csv")
data1 [,1]=as.character(data1[,1])
k=1:23
k=24-k
data1 = data1 [k,]
## read the color
colors1 = as.character(data1[,1])
## read the Z-summary
x1= data1[,3]
##read the name
xlimits=c(0,35)
title="Preservation of CN modules in CB data "
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x1, colors1, axes = F, main=title,cex.main=1, names= colors1 ,xlim=xlimits, cex.names=1, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.1 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.2*(box[2]-box[1]),box[4]+0.07*(box[4]-box[3]), "d", adj = c(0.5, 0.5), xpd = TRUE, cex = 2,)
## 5 Plot how well CB modules preserved in CTX data
data1=read.csv ("5.ctx.cb.with.cb.module.csv")
data1 [,1]=as.character(data1[,1])
k=1:23
k=24-k
data1 = data1 [k,]
## read the color
colors1 = as.character(data1[,1])
## read the Z-summary
x1= data1[,3]
##read the name
xlimits=c(0,35)
title="Preservation of CB modules in CTX data "
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x1, colors1, axes = F, main=title,cex.main=1, names= colors1 ,xlim=xlimits, cex.names=1, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.1 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.2*(box[2]-box[1]),box[4]+0.07*(box[4]-box[3]), "e", adj = c(0.5, 0.5), xpd = TRUE, cex = 2,)
## 6 Plot how well CB modules preserved in CN data
data1=read.csv (".cb.with.cb.module.csv")
data1 [,1]=as.character(data1[,1])
k=1:23
k=24-k
data1 = data1 [k,]
## read the color
colors1 = as.character(data1[,1])
## read the Z-summary
x1= data1[,3]
##read the name
xlimits=c(0,35)
title="Preservation of CB modules in CN data "
library(WGCNA)
par(mar = c(3, 3, 3, 8))
par(mgp = c(4,2.8,0))
barplotIt(x1, colors1, axes = F, main=title,cex.main=1, names= colors1 ,xlim=xlimits, cex.names=1, family = "sans", xlab="sans",ylab="sans")
abline(v=5,col="red")
abline(v=10,col="red")
par(mgp = c(2.5, 1.0, 0))
axis(1)
box = par("usr");
text(mean(box[c(1:2)]), box[3] - 0.1 * (box[4]-box[3]), "Preservation Z score", adj = c(0.5, 0.5), xpd = TRUE)
text(box[1]-0.2*(box[2]-box[1]),box[4]+0.07*(box[4]-box[3]), "f", adj = c(0.5, 0.5), xpd = TRUE, cex = 2,)
dev.off()
................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related searches
- steve harvey 21 day detox
- simon sinek s father steve sinek
- steve jobs speech about life
- ucla health new patient forms
- ucla master of financial engineering
- steve spangler science experiment library
- steve harvey detox cleanse
- steve landers toyota
- steve harvey closing remarks today
- steve harvey s devotional
- ucla tuition resident
- ucla credential program