#-------------------------------------------------------------------- #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #Change working directory to the folder with downloaded files #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #-------------------------------------------------------------------- #-----Call required libraries--------- library(stats) library(sjstats) library(Rcpp) library(ggplot2) library(ggpubr) library(dplyr) library(caTools) library(overlapping) library(lattice) library(pwr) library(stringr) #--------Read in data------------------ whole_cc <- read.csv("data_with_corpus_callosum_segmented_as_single_ROI.csv", header = TRUE) #cohort HC whole corpus callosum data regions_cc <- read.csv("data_for_corpus_callosum_segmented_by_subregions.csv", header = TRUE) #cohort HC sub-region data gender_cc <- read.csv("data_of_gender_matched_group_HC_sub.csv", header = TRUE) #data for cohort HC_sub (whole corpus callosum) seg_cc <- read.csv("data_for_segmentation_methods_comparison.csv", header = TRUE) #data for Figure 2: 10 matched subjects from HC_sub dataHC_whole<-data.frame(subset(whole_cc, group %in% c("HC"))) dataHC_region<-data.frame(subset(regions_cc, group %in% c("HC"))) l <- read.csv("longitudinal.csv", header = TRUE) long<-data.frame(subset(l, group %in% c("HC","MCI"))) #============================================================================================================================================================================================ #============================================================================================================================================================================================ #-----------------Analysis for figure 2------------------------------------- #-----------------using 10 subjects from HC_sub-------------- seg_cc$type<-factor(seg_cc$type,levels=c("Circular ROI in axial plane (median)", "ROI in sagittal plane (median)", "ROI in sagittal plane (mean)", "Atlas ROI (mean)", "Atlas ROI (median)")) ggplot(seg_cc,aes(x=type,y=fa))+ geom_line(size=1, linetype="dashed", aes(group=id,color=id))+ geom_point(size=5, aes(color=id))+ geom_smooth(method=lm,se=FALSE)+theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=0.5), legend.position='none')+ labs(y=" ",x=" ")+ scale_x_discrete(labels = function(type) str_wrap(type, width = 12)) ggplot(seg_cc, aes(x=type, y=diff)) + geom_boxplot(width=0.25,fill='#A4A4A4') + theme_classic() + theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=0.5), legend.position='none')+ labs(y=" ",x=" ")+ scale_x_discrete(labels = function(type) str_wrap(type, width = 12)) #============================================================================================================================================================================================ #============================================================================================================================================================================================ #------------------Analysis for figure 3 -------------------------------------------------- #------------------Using data from cohort HC------------------------------------------------------------------ #-----------------Figures 3 a - d ----------------------------------------------- #display linear regression plots faplot<-ggplot(dataHC_whole,aes(y=fa,x=age))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,colour="black",linetype="longdash",size=1.5)+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.8,0.95))+ labs(y=" ", x=" ") dev.new() faplot mdplot<-ggplot(dataHC_whole,aes(y=md,x=age))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,colour="black",linetype="longdash",size=1.5)+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.05,0.95))+labs(y=" ",x=" ") dev.new() mdplot rdplot<-ggplot(dataHC_whole,aes(y=rd,x=age))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,colour="black",linetype="longdash",size=1.5)+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.05,0.95))+labs(y=" ",x=" ") dev.new() rdplot adplot<-ggplot(dataHC_whole,aes(y=ad,x=age))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,colour="black",linetype="longdash",size=1.5)+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.05,0.95))+labs(y=" ",x=" ") dev.new() adplot #--------------------------Figures 3 e - h----------------------------------- faplot<-ggplot(dataHC_region,aes(y=fa,x=age,color=factor(region)))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,linetype="longdash",size=1.5)+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position=c(0.5,0.5))+ labs(y=" ", x=" ",color="Region") dev.new() faplot mdplot<-ggplot(dataHC_region,aes(y=md,x=age,color=factor(region)))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,linetype="longdash",size=1.5)+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+labs(y=" ",x=" ",color=" ") dev.new() mdplot rdplot<-ggplot(dataHC_region,aes(y=rd,x=age,color=factor(region)))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,linetype="longdash",size=1.5)+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+labs(y=" ",x=" ",color=" ") dev.new() rdplot adplot<-ggplot(dataHC_region,aes(y=ad,x=age,color=factor(region)))+geom_point(size=2)+ stat_smooth(method="lm",se=FALSE,linetype="longdash",size=1.5)+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+labs(y=" ",x=" ",color=" ") dev.new() adplot #-----------------------------Figures 3 i - l--------------------------- dataHC_region$region<-as.factor(dataHC_region$region) fage<-cut(dataHC_region$age, breaks=4) #convert numerical 'age' to categories ancovaFA3<-aov(fa~region+fage,data=dataHC_region) ancovaMD3<-aov(md~region+fage,data=dataHC_region) ancovaDR3<-aov(rd~region+fage,data=dataHC_region) ancovaDA3<-aov(ad~region+fage,data=dataHC_region) tukeyfa<-as.data.frame(TukeyHSD(ancovaFA3)$region) colnames(tukeyfa)[4]<-"padj" tukeyfa$pair=rownames(tukeyfa) p<-as.matrix(tukeyfa$padj) dev.new() ggplot(tukeyfa,aes(color=cut(p,c(0, 0.001, 0.05, 1), label=c("p < 0.001","p < 0.05","Non-significant"))))+ geom_errorbar(aes(pair,ymin=lwr,ymax=upr),width=0.5,size=2)+ geom_point(aes(pair,diff),size=2)+ labs(color="")+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.65,0.18))+coord_flip()+labs(y=" ",x=" ")+ scale_color_manual(name = " ", values = c("p < 0.001" = "red", "p < 0.05" = "blueviolet", "Non-significant" = "gray60")) tukeymd<-as.data.frame(TukeyHSD(ancovaMD3)$region) colnames(tukeymd)[4]<-"padj" tukeymd$pair=rownames(tukeymd) p<-as.matrix(tukeymd$padj) dev.new() ggplot(tukeymd,aes(color=cut(p,c(0, 0.001, 0.05, 1), label=c("p < 0.001","p < 0.05","Non-significant"))))+ geom_errorbar(aes(pair,ymin=lwr,ymax=upr),width=0.5,size=2)+ geom_point(aes(pair,diff),size=2)+ labs(color="")+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.65,0.18))+coord_flip()+labs(y=" ",x=" ")+ scale_color_manual(name = " ",values = c("p < 0.001" = "red","p < 0.05" = "blueviolet","Non-significant" = "gray60")) tukeyrd<-as.data.frame(TukeyHSD(ancovaDR3)$region) colnames(tukeyrd)[4]<-"padj" tukeyrd$pair=rownames(tukeyrd) p<-as.matrix(tukeyrd$padj) dev.new() ggplot(tukeyrd,aes(color=cut(p,c(0, 0.001, 0.05, 1), label=c("p < 0.001","p < 0.05","Non-significant"))))+ geom_errorbar(aes(pair,ymin=lwr,ymax=upr),width=0.5,size=2)+ geom_point(aes(pair,diff),size=2)+ labs(color="")+theme_classic()+ theme(text=element_text(size=30,face="bold"), axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.65,0.18))+coord_flip()+labs(y=" ",x=" ")+ scale_color_manual(name = " ", values = c("p < 0.001" = "red","p < 0.05" = "blueviolet","Non-significant" = "gray60")) tukeyad<-as.data.frame(TukeyHSD(ancovaDA3)$region) colnames(tukeyad)[4]<-"padj" tukeyad$pair=rownames(tukeyad) p<-as.matrix(tukeyad$padj) dev.new() ggplot(tukeyad,aes(color=cut(p,c(0, 0.001, 0.05, 1), label=c("p < 0.001","p < 0.05","Non-significant"))))+ geom_errorbar(aes(pair,ymin=lwr,ymax=upr),width=0.5,size=2)+ geom_point(aes(pair,diff),size=2)+ labs(color="")+theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1), legend.justification=c(0,1),legend.position=c(0.65,0.18))+coord_flip()+labs(y=" ",x=" ")+ scale_color_manual(name = " ", values = c("p < 0.001" = "red","p < 0.05" = "blueviolet","Non-significant" = "gray60")) #============================================================================================================================================================================================ #============================================================================================================================================================================================ #-------------------------------Analysis for figure 4 (b)------------------------- #-------------------------------Using data from cohort HC----------------------- dev.new() ggplot(dataHC_region, aes(x=region, y=fa, group=region)) + geom_boxplot(color="gray60",lwd=1.2) + theme_classic()+ theme(text=element_text(size=40,face="bold"), axis.text.x=element_text(hjust=1))+ labs(x=" ")+ stat_summary(fun.y=mean, geom="line", aes(group=1),lwd=2,color="red") + stat_summary(fun.y=mean, geom="point",size=2) dev.new() ggplot(dataHC_region, aes(x=region, y=rd, group=region)) + geom_boxplot(color="gray60",lwd=1.2) + theme_classic()+ theme(text=element_text(size=40,face="bold"), axis.text.x=element_text(hjust=1))+ labs(x=" ")+ stat_summary(fun.y=mean, geom="line", aes(group=1),lwd=2,color="red") + stat_summary(fun.y=mean, geom="point",size=2) #============================================================================================================================================================================================ #============================================================================================================================================================================================ #------------------------------Plots for Figure 5------------------------ #------------------------------Using data from cohort HC and MCI------------------------------ #----------------------------Figure 5 (a)-------------------------------- #display linear regression plots separated by disease group faplot2<-ggplot(whole_cc,aes(y=fa,x=age))+geom_point(aes(shape=factor(group),size=2, stroke=2))+ scale_shape_manual(values=c(1,18))+ stat_smooth(method="lm",colour="black",se=FALSE,size=2,aes(linetype=factor(group)))+theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+ labs(y=" ",x=" ",color=" ") dev.new() faplot2 mdplot2<-ggplot(whole_cc,aes(y=md,x=age))+geom_point(aes(shape=factor(group),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ stat_smooth(method="lm",colour="black",se=FALSE,size=2,aes(linetype=factor(group)))+theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+labs(y=" ",x=" ",color=" ") dev.new() mdplot2 rdplot2<-ggplot(whole_cc,aes(y=rd,x=age))+geom_point(aes(shape=factor(group),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ stat_smooth(method="lm",colour="black",se=FALSE,size=2,aes(linetype=factor(group)))+theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+labs(y=" ",x=" ",color=" ") dev.new() rdplot2 adplot2<-ggplot(whole_cc,aes(y=ad,x=age))+geom_point(aes(shape=factor(group),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ stat_smooth(method="lm",colour="black",se=FALSE,size=2,aes(linetype=factor(group)))+theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+labs(y=" ",x=" ",color=" ") dev.new() adplot2 ydensity <- ggplot(regions_cc, aes(fa, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity ydensity <- ggplot(regions_cc, aes(md, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity ydensity <- ggplot(regions_cc, aes(rd, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity ydensity <- ggplot(regions_cc, aes(ad, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity #---------------------------Figure 5 (b)--------------------------------------------------------------------------------------------------- ydensity <- ggplot(regions_cc, aes(fa, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ facet_grid(cols = vars(region))+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity ydensity <- ggplot(regions_cc, aes(md, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ facet_grid(cols = vars(region))+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity ydensity <- ggplot(regions_cc, aes(rd, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ facet_grid(cols = vars(region))+ theme(text=element_text(size=45,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none",axis.title.y=element_blank()) dev.new() ydensity+scale_x_continuous(breaks=c(0.3,0.6,0.9)) ydensity <- ggplot(regions_cc, aes(ad, fill=group)) + geom_density(alpha=.3) + scale_fill_manual(values = c("HC"="white", "MCI"="grey7"))+ facet_grid(cols = vars(region))+ theme(text=element_text(size=45,face="bold"),axis.text.x=element_text(hjust=1),,legend.justification=c(0,1), legend.position="none", axis.title.y=element_blank()) dev.new() ydensity+scale_x_continuous(breaks=c(1.5,1.8,2.1)) #============================================================================================================================================================================================ #============================================================================================================================================================================================ #----------Supplementary figure S1 ------------------------------------------------------- ggplot(long,aes(y=meanfa,x=time))+geom_line(size=3,aes(colour=factor(group),group=factor(group)))+ geom_point(size=7,aes(colour=factor(group)))+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=0.5),legend.justification=c(0,1), legend.position=c(0.5,0.5))+ labs(x=" ", y=" ")+ scale_x_discrete(limits=c("b","3","6","12")) ggplot(long,aes(y=meanmd,x=time))+geom_line(size=3,aes(colour=factor(group),group=factor(group)))+ geom_point(size=7,aes(colour=factor(group)))+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=0.5),legend.justification=c(0,1), legend.position=c(0.5,0.5))+ labs(x=" ", y=" ")+ scale_x_discrete(limits=c("b","3","6","12")) ggplot(long,aes(y=meanrd,x=time))+geom_line(size=3,aes(colour=factor(group),group=factor(group)))+ geom_point(size=7,aes(colour=factor(group)))+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=0.5),legend.justification=c(0,1), legend.position=c(0.5,0.5))+ labs(x=" ", y=" ")+ scale_x_discrete(limits=c("b","3","6","12")) ggplot(long,aes(y=meanaxd,x=time))+geom_line(size=3,aes(colour=factor(group),group=factor(group)))+ geom_point(size=7,aes(colour=factor(group)))+ theme_classic()+ theme(text=element_text(size=30,face="bold"),axis.text.x=element_text(hjust=0.5),legend.justification=c(0,1), legend.position=c(0.5,0.5))+ labs(x=" ", y=" ")+ scale_x_discrete(limits=c("b","3","6","12")) #----------Supplementary figure S2 ------------------------------------------------------- #----------testing effects of gender on ageing using cohort HC_sub---------------------------- #plots produced are without labels etc. for formatting reasons. Please note the name of the plot to understand which parameter it represents #display linear regression plots faplot<-ggplot(gender_cc,aes(y=fa,x=age))+geom_point(aes(shape=factor(gender),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ geom_smooth(method="lm",colour="black",se=FALSE,size=1.5,aes(linetype=factor(gender)))+ theme_classic()+ theme(text=element_text(size=40,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position=c(0.8,0.8))+ labs(y=" ", x=" ",color=" ") dev.new() faplot mdplot<-ggplot(gender_cc,aes(y=md,x=age))+geom_point(aes(shape=factor(gender),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ geom_smooth(method="lm",colour="black",se=FALSE,size=1.5,aes(linetype=factor(gender)))+ theme_classic()+ theme(text=element_text(size=40,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+ labs(y=" ", x=" ",color=" ")+scale_y_continuous(breaks=seq(0.8, 1, by=0.1)) dev.new() mdplot rdplot<-ggplot(gender_cc,aes(y=rd,x=age))+geom_point(aes(shape=factor(gender),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ geom_smooth(method="lm",colour="black",se=FALSE,size=1.5,aes(linetype=factor(gender)))+ theme_classic()+ theme(text=element_text(size=40,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+ labs(y=" ", x=" ",color=" ")+ scale_y_continuous(breaks=seq(0.4, 0.65, by=0.1)) dev.new() rdplot adplot<-ggplot(gender_cc,aes(y=ad,x=age))+geom_point(aes(shape=factor(gender),size=2,stroke=2))+ scale_shape_manual(values=c(1,18))+ geom_smooth(method="lm",colour="black",se=FALSE,size=1.5,aes(linetype=factor(gender)))+ theme_classic()+ theme(text=element_text(size=40,face="bold"),axis.text.x=element_text(hjust=1),legend.justification=c(0,1), legend.position="none")+ labs(y=" ", x=" ",color=" ") dev.new() adplot #===============================================Supplementary figure S3====================================================================================================================== #============================================================================================================================================================================================ faplot<-ggplot(whole_cc,aes(y=md,x=rd,colour=group))+geom_point(size=3)+ stat_smooth(method="lm",se=FALSE,size=3)+ theme_classic()+ theme(text=element_text(size=60,face="bold"),axis.text.x=element_text(hjust=1),legend.position="none")+ labs(y=" ", x=" ")+ scale_color_manual(name = " ", values = c("HC" = "red","MCI" = "black")) dev.new() faplot faplot<-ggplot(whole_cc,aes(y=md,x=ad,colour=group))+geom_point(size=3)+ stat_smooth(method="lm",se=FALSE,size=3)+ theme_classic()+ theme(text=element_text(size=60,face="bold"),axis.text.x=element_text(hjust=1),legend.position="none")+ labs(y=" ", x=" ")+ scale_color_manual(name = " ", values = c("HC" = "red","MCI" = "black")) dev.new() faplot faplot<-ggplot(whole_cc,aes(y=md,x=sum_rd_ad,colour=group))+geom_point(size=3)+ stat_smooth(method="lm",se=FALSE,size=3)+ theme_classic()+ theme(text=element_text(size=60,face="bold"),axis.text.x=element_text(hjust=1),legend.position="none")+ labs(y=" ", x=" ")+ scale_color_manual(name = " ", values = c("HC" = "red","MCI" = "black")) dev.new() faplot