Steve Horvath UCLA – Horvath Lab UCLA - How to calculate tan squared

[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.

Literature Lottery

Steve Horvath UCLA – Horvath Lab UCLA

To fulfill the demand for quickly locating and searching documents.

Related download

Related searches