Script_06_aVMCs_validation.Rmd

---
title: "aVMCs validation"
author: "Yunfeng Liu"
date: "1/16/2023"
output: html_document
---

Here,we would like to replicate our aDMCs analysis on X chromosome using two external datasets.

1.Blood DNAm data generated by Johansson et al.: "Continuous Aging of the Human DNA Methylome Throughout the Human Lifespan.PLoS One."
GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE87571

2.Monocytes DNAm data generated by Reynolds et al.: "Age-related variations in the methylome associated with gene expression in human monocytes and T cells. Nat Commun."
GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56046

```{r}
#set library path.
.libPaths("~/researchdrive/yliu/Rlibs")
```

```{r}
# Load data
load(file = "Output_05_DGLM_Johansson_Blood_Xfemales_converge.RData")
load(file = "Output_05_DGLM_Johansson_Blood_Xmales_converge.RData")
load(file = "Output_05_DGLM_Reynolds_Mono_Xfemales_converge.RData")
load(file = "Output_05_DGLM_Reynolds_Mono_Xmales_converge.RData")
load(file = "Output_02_aDMCs_effectsize_BIOS.RData")
```


```{r}
#### aVMCs validation on Johansson blood datasets
## Correct bias and inflation for dispersion model.
# Xfemales
set.seed(1)
bc_DGLM_Johansson_Blood_Xfemales<- bacon(DGLM_Johansson_Blood_Xfemales_converge[,7])

# Calculate bias and inflation
bias(bc_DGLM_Johansson_Blood_Xfemales)         # 0.15  
inflation(bc_DGLM_Johansson_Blood_Xfemales)    # 0.77

# Save histogram plot of test statistics 
tiff("Output_06_Bacon_DGLM_Johansson_Blood_Xfemales.tiff", units="in", width=8, height=6, res=600,compression = 'lzw')
fit(bc_DGLM_Johansson_Blood_Xfemales,xlab="t-statistics",main="aVMCs in Johansson Blood females",n=100)
dev.off()

#Extract the BACON-adjusted t-statistics and p-values
tstats_DGLM_Johansson_Blood_Xfemales<-tstat(bc_DGLM_Johansson_Blood_Xfemales)
pvals_DGLM_Johansson_Blood_Xfemales <-pval(bc_DGLM_Johansson_Blood_Xfemales)
rownames(pvals_DGLM_Johansson_Blood_Xfemales)<-make.names(rownames(DGLM_Johansson_Blood_Xfemales_converge),unique=TRUE) 

# Extract BACON-adjusted effectsize.
set.seed(1)
bc_DGLM_Johansson_Blood_Xfemales_es<-bacon(NULL,DGLM_Johansson_Blood_Xfemales_converge[,5],DGLM_Johansson_Blood_Xfemales_converge[,6])
es_DGLM_Johansson_Blood_Xfemales<-es(bc_DGLM_Johansson_Blood_Xfemales_es)

# Extract BACON-adjusted standard error.
set.seed(1)
bc_DGLM_Johansson_Blood_Xfemales_se<-bacon(NULL,DGLM_Johansson_Blood_Xfemales_converge[,5],DGLM_Johansson_Blood_Xfemales_converge[,6])
se_DGLM_Johansson_Blood_Xfemales<-se(bc_DGLM_Johansson_Blood_Xfemales_se)

# Summary Bacon statistics
dat_DGLM_Johansson_Blood_Xfemales_bacon_disp<-as.data.frame(cbind(es_DGLM_Johansson_Blood_Xfemales,se_DGLM_Johansson_Blood_Xfemales,tstats_DGLM_Johansson_Blood_Xfemales,pvals_DGLM_Johansson_Blood_Xfemales))
colnames(dat_DGLM_Johansson_Blood_Xfemales_bacon_disp)<-c("bacon_effectsize_disp","bacon_std.error_disp","bacon_tstatistics_disp","bacon_pval_disp")
rownames(dat_DGLM_Johansson_Blood_Xfemales_bacon_disp)<-make.names(rownames(DGLM_Johansson_Blood_Xfemales_converge),unique = TRUE)

# check females effecsize of aVMCs between BIOS and Johansson 
dat_Johansson_Blood_Xfemales_aVMCs<-dat_DGLM_Johansson_Blood_Xfemales_bacon_disp[rownames(df_aVMCs_effectsize_BIOS),]
df_Johansson_validation_females<-cbind(dat_Johansson_Blood_Xfemales_aVMCs[,1],df_aVMCs_effectsize_BIOS[,-2])
colnames(df_Johansson_validation_females)<-c("Johansson","BIOS","aVMCs")

# Xmales
set.seed(1)
bc_DGLM_Johansson_Blood_Xmales<- bacon(DGLM_Johansson_Blood_Xmales_converge[,7])

# Calculate bias and inflation
bias(bc_DGLM_Johansson_Blood_Xmales)         # 0.39 
inflation(bc_DGLM_Johansson_Blood_Xmales)    # 0.91

# Save histogram plot of test statistics 
tiff("Output_06_Bacon_DGLM_Johansson_Blood_Xmales.tiff", units="in", width=8, height=6, res=600,compression = 'lzw')
fit(bc_DGLM_Johansson_Blood_Xmales,xlab="t-statistics",main="aVMCs in Johansson Blood males",n=100)
dev.off()

#Extract the BACON-adjusted t-statistics and p-values
tstats_DGLM_Johansson_Blood_Xmales<-tstat(bc_DGLM_Johansson_Blood_Xmales)
pvals_DGLM_Johansson_Blood_Xmales <-pval(bc_DGLM_Johansson_Blood_Xmales)
rownames(pvals_DGLM_Johansson_Blood_Xmales)<-make.names(rownames(DGLM_Johansson_Blood_Xmales_converge),unique=TRUE) 

# Extract BACON-adjusted effectsize.
set.seed(1)
bc_DGLM_Johansson_Blood_Xmales_es<-bacon(NULL,DGLM_Johansson_Blood_Xmales_converge[,5],DGLM_Johansson_Blood_Xmales_converge[,6])
es_DGLM_Johansson_Blood_Xmales<-es(bc_DGLM_Johansson_Blood_Xmales_es)

# Extract BACON-adjusted standard error.
set.seed(1)
bc_DGLM_Johansson_Blood_Xmales_se<-bacon(NULL,DGLM_Johansson_Blood_Xmales_converge[,5],DGLM_Johansson_Blood_Xmales_converge[,6])
se_DGLM_Johansson_Blood_Xmales<-se(bc_DGLM_Johansson_Blood_Xmales_se)

# Summary Bacon statistics
dat_DGLM_Johansson_Blood_Xmales_bacon_disp<-as.data.frame(cbind(es_DGLM_Johansson_Blood_Xmales,se_DGLM_Johansson_Blood_Xmales,tstats_DGLM_Johansson_Blood_Xmales,pvals_DGLM_Johansson_Blood_Xmales))
colnames(dat_DGLM_Johansson_Blood_Xmales_bacon_disp)<-c("bacon_effectsize_disp","bacon_std.error_disp","bacon_tstatistics_disp","bacon_pval_disp")
rownames(dat_DGLM_Johansson_Blood_Xmales_bacon_disp)<-make.names(rownames(DGLM_Johansson_Blood_Xmales_converge),unique = TRUE)

# check males effecsize of aVMCs between BIOS and Johansson 
dat_Johansson_Blood_Xmales_aVMCs<-dat_DGLM_Johansson_Blood_Xmales_bacon_disp[rownames(df_aVMCs_effectsize_BIOS),]
df_Johansson_validation_males<-cbind(dat_Johansson_Blood_Xmales_aVMCs[,1],df_aVMCs_effectsize_BIOS[,-1])
colnames(df_Johansson_validation_males)<-c("Johansson","BIOS","aVMCs")
```

```{r}
# save the validation data of effect size 
save(df_Johansson_validation_females,df_Johansson_validation_males,file="Output_06_Johansson_validation.RData")
```

```{r}
# Extract 1098 female p-values and 39 male p-values.
dat_Johansson_Blood_Xfemales_aVMCs$aVMCs<-df_aVMCs_effectsize_BIOS$aVMCs
df_aVMCs_females_Johansson_Blood<-filter(dat_Johansson_Blood_Xfemales_aVMCs,aVMCs=="female-specific"|aVMCs=="both-sex")
pvals_aVMCs_Xfemales_Johansson_Blood<-as.data.frame(pvals_DGLM_Johansson_Blood_Xfemales[rownames(df_aVMCs_females_Johansson_Blood),])
colnames(pvals_aVMCs_Xfemales_Johansson_Blood)<-"pval_females"

dat_Johansson_Blood_Xmales_aVMCs$aVMCs<-df_aVMCs_effectsize_BIOS$aVMCs
df_aVMCs_males_Johansson_Blood<-filter(dat_Johansson_Blood_Xmales_aVMCs,aVMCs=="male-specific"|aVMCs=="both-sex")
pvals_aVMCs_Xmales_Johansson_Blood<-as.data.frame(pvals_DGLM_Johansson_Blood_Xmales[rownames(df_aVMCs_males_Johansson_Blood),])
colnames(pvals_aVMCs_Xmales_Johansson_Blood)<-"pval_males"

# Perform FDR correction for 1137 pvals(1098+39)
padjs_aVMCs_Johansson_Blood<-as.data.frame(p.adjust(c(pvals_aVMCs_Xfemales_Johansson_Blood$pval_females,pvals_aVMCs_Xmales_Johansson_Blood$pval_males),method="BH"))
colnames(padjs_aVMCs_Johansson_Blood)<-"bacon_p.adj_dispersion"

# check females (how many of 1098) and males (how many of 39) separately how many were FDR significant 
# 1034 female aVMCs are still FDR significant in Hannum blood data.
padjs_aVMCs_females_Johansson_Blood<-as.data.frame(padjs_aVMCs_Johansson_Blood[1:1098,])
colnames(padjs_aVMCs_females_Johansson_Blood)<-"bacon_p.adj_disp"
df_aVMCs_females_Johansson_Blood$bacon_p.adj_disp<-padjs_aVMCs_females_Johansson_Blood$bacon_p.adj_disp

df_aVMCs_females_Johansson_Blood<-filter(df_aVMCs_females_Johansson_Blood,bacon_p.adj_disp<0.05)
df_aVMCs_females_effectsize_up_Johansson<-filter(df_aVMCs_females_Johansson_Blood,bacon_effectsize_disp>0) # 1030
df_aVMCs_females_effectsize_down_Johansson<-filter(df_aVMCs_females_Johansson_Blood,bacon_effectsize_disp<0) # 4

# 37 male aVMCs are still FDR significant in Hannum blood data.
padjs_aVMCs_males_Johansson_Blood<-as.data.frame(padjs_aVMCs_Johansson_Blood[1099:1137,])
colnames(padjs_aVMCs_males_Johansson_Blood)<-"bacon_p.adj_disp"
df_aVMCs_males_Johansson_Blood$bacon_p.adj_disp<-padjs_aVMCs_males_Johansson_Blood$bacon_p.adj_disp

df_aVMCs_males_Johansson_Blood<-filter(df_aVMCs_males_Johansson_Blood,bacon_p.adj_disp<0.05)
df_aVMCs_males_effectsize_up_Johansson<-filter(df_aVMCs_males_Johansson_Blood,bacon_effectsize_disp>0) # 36
df_aVMCs_males_effectsize_down_Johansson<-filter(df_aVMCs_males_Johansson_Blood,bacon_effectsize_disp<0) # 1

save(df_aVMCs_females_effectsize_up_Johansson,df_aVMCs_females_effectsize_down_Johansson,df_aVMCs_females_Johansson_Blood,df_aVMCs_males_effectsize_up_Johansson,df_aVMCs_males_effectsize_down_Johansson,df_aVMCs_males_Johansson_Blood,file = "Output_06_Johansson_Catogue of aVMCs.RData")
```

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

```{r}
#### aVMCs validation on Reynolds Monocytes datasets
## Correct bias and inflation for dispersion model.
# Xfemales
set.seed(1)
bc_DGLM_Reynolds_Mono_Xfemales<- bacon(DGLM_Reynolds_Mono_Xfemales_converge[,7])

# Calculate bias and inflation
bias(bc_DGLM_Reynolds_Mono_Xfemales)         # 0.17  
inflation(bc_DGLM_Reynolds_Mono_Xfemales)    # 0.71

# Save histogram plot of test statistics 
tiff("Output_06_Bacon_DGLM_Reynolds_Mono_Xfemales.tiff", units="in", width=8, height=6, res=600,compression = 'lzw')
fit(bc_DGLM_Reynolds_Mono_Xfemales,xlab="t-statistics",main="aVMCs in Reynolds Monocytes females",n=100)
dev.off()

#Extract the BACON-adjusted t-statistics and p-values
tstats_DGLM_Reynolds_Mono_Xfemales<-tstat(bc_DGLM_Reynolds_Mono_Xfemales)
pvals_DGLM_Reynolds_Mono_Xfemales <-pval(bc_DGLM_Reynolds_Mono_Xfemales)
rownames(pvals_DGLM_Reynolds_Mono_Xfemales)<-make.names(rownames(DGLM_Reynolds_Mono_Xfemales_converge),unique=TRUE) 

# Extract BACON-adjusted effectsize.
set.seed(1)
bc_DGLM_Reynolds_Mono_Xfemales_es<-bacon(NULL,DGLM_Reynolds_Mono_Xfemales_converge[,5],DGLM_Reynolds_Mono_Xfemales_converge[,6])
es_DGLM_Reynolds_Mono_Xfemales<-es(bc_DGLM_Reynolds_Mono_Xfemales_es)

# Extract BACON-adjusted standard error.
set.seed(1)
bc_DGLM_Reynolds_Mono_Xfemales_se<-bacon(NULL,DGLM_Reynolds_Mono_Xfemales_converge[,5],DGLM_Reynolds_Mono_Xfemales_converge[,6])
se_DGLM_Reynolds_Mono_Xfemales<-se(bc_DGLM_Reynolds_Mono_Xfemales_se)

# Summary Bacon statistics
dat_DGLM_Reynolds_Mono_Xfemales_bacon_disp<-as.data.frame(cbind(es_DGLM_Reynolds_Mono_Xfemales,se_DGLM_Reynolds_Mono_Xfemales,tstats_DGLM_Reynolds_Mono_Xfemales,pvals_DGLM_Reynolds_Mono_Xfemales))
colnames(dat_DGLM_Reynolds_Mono_Xfemales_bacon_disp)<-c("bacon_effectsize_disp","bacon_std.error_disp","bacon_tstatistics_disp","bacon_pval_disp")
rownames(dat_DGLM_Reynolds_Mono_Xfemales_bacon_disp)<-make.names(rownames(DGLM_Reynolds_Mono_Xfemales_converge),unique = TRUE)

# check females effecsize of aVMCs between BIOS and Reynolds 
dat_Reynolds_Mono_Xfemales_aVMCs<-dat_DGLM_Reynolds_Mono_Xfemales_bacon_disp[rownames(df_aVMCs_effectsize_BIOS),]
df_Reynolds_validation_females<-cbind(dat_Reynolds_Mono_Xfemales_aVMCs[,1],df_aVMCs_effectsize_BIOS[,-2])
colnames(df_Reynolds_validation_females)<-c("Reynolds","BIOS","aVMCs")

# Xmales
set.seed(1)
bc_DGLM_Reynolds_Mono_Xmales<- bacon(DGLM_Reynolds_Mono_Xmales_converge[,7])

# Calculate bias and inflation
bias(bc_DGLM_Reynolds_Mono_Xmales)         # 0.30  
inflation(bc_DGLM_Reynolds_Mono_Xmales)    # 0.79

# Save histogram plot of test statistics 
tiff("Output_06_Bacon_DGLM_Reynolds_Mono_Xmales.tiff", units="in", width=8, height=6, res=600,compression = 'lzw')
fit(bc_DGLM_Reynolds_Mono_Xmales,xlab="t-statistics",main="aVMCs in Reynolds Monocytes males",n=100)
dev.off()

#Extract the BACON-adjusted t-statistics and p-values
tstats_DGLM_Reynolds_Mono_Xmales<-tstat(bc_DGLM_Reynolds_Mono_Xmales)
pvals_DGLM_Reynolds_Mono_Xmales <-pval(bc_DGLM_Reynolds_Mono_Xmales)
rownames(pvals_DGLM_Reynolds_Mono_Xmales)<-make.names(rownames(DGLM_Reynolds_Mono_Xmales_converge),unique=TRUE) 

# Extract BACON-adjusted effectsize.
set.seed(1)
bc_DGLM_Reynolds_Mono_Xmales_es<-bacon(NULL,DGLM_Reynolds_Mono_Xmales_converge[,5],DGLM_Reynolds_Mono_Xmales_converge[,6])
es_DGLM_Reynolds_Mono_Xmales<-es(bc_DGLM_Reynolds_Mono_Xmales_es)

# Extract BACON-adjusted standard error.
set.seed(1)
bc_DGLM_Reynolds_Mono_Xmales_se<-bacon(NULL,DGLM_Reynolds_Mono_Xmales_converge[,5],DGLM_Reynolds_Mono_Xmales_converge[,6])
se_DGLM_Reynolds_Mono_Xmales<-se(bc_DGLM_Reynolds_Mono_Xmales_se)

# Summary Bacon statistics
dat_DGLM_Reynolds_Mono_Xmales_bacon_disp<-as.data.frame(cbind(es_DGLM_Reynolds_Mono_Xmales,se_DGLM_Reynolds_Mono_Xmales,tstats_DGLM_Reynolds_Mono_Xmales,pvals_DGLM_Reynolds_Mono_Xmales))
colnames(dat_DGLM_Reynolds_Mono_Xmales_bacon_disp)<-c("bacon_effectsize_disp","bacon_std.error_disp","bacon_tstatistics_disp","bacon_pval_disp")
rownames(dat_DGLM_Reynolds_Mono_Xmales_bacon_disp)<-make.names(rownames(DGLM_Reynolds_Mono_Xmales_converge),unique = TRUE)

# check males effecsize of aVMCs between BIOS and Reynolds 
dat_Reynolds_Mono_Xmales_aVMCs<-dat_DGLM_Reynolds_Mono_Xmales_bacon_disp[rownames(df_aVMCs_effectsize_BIOS),]
df_Reynolds_validation_males<-cbind(dat_Reynolds_Mono_Xmales_aVMCs[,1],df_aVMCs_effectsize_BIOS[,-1])
colnames(df_Reynolds_validation_males)<-c("Reynolds","BIOS","aVMCs")
```

```{r}
# save the validation data of effect size 
save(df_Reynolds_validation_females,df_Reynolds_validation_males,file="Output_06_Reynolds_validation.RData")
```

```{r}
# Extract 1038 female p-values and 39 male p-values.
dat_Reynolds_Mono_Xfemales_aVMCs$aVMCs<-df_aVMCs_effectsize_BIOS$aVMCs
df_aVMCs_females_Reynolds_Mono<-filter(dat_Reynolds_Mono_Xfemales_aVMCs,aVMCs=="female-specific"|aVMCs=="both-sex")
pvals_aVMCs_Xfemales_Reynolds_Mono<-as.data.frame(pvals_DGLM_Reynolds_Mono_Xfemales[rownames(df_aVMCs_females_Reynolds_Mono),])
colnames(pvals_aVMCs_Xfemales_Reynolds_Mono)<-"pval_females"

dat_Reynolds_Mono_Xmales_aVMCs$aVMCs<-df_aVMCs_effectsize_BIOS$aVMCs
df_aVMCs_males_Reynolds_Mono<-filter(dat_Reynolds_Mono_Xmales_aVMCs,aVMCs=="male-specific"|aVMCs=="both-sex")
pvals_aVMCs_Xmales_Reynolds_Mono<-as.data.frame(pvals_DGLM_Reynolds_Mono_Xmales[rownames(df_aVMCs_males_Reynolds_Mono),])
colnames(pvals_aVMCs_Xmales_Reynolds_Mono)<-"pval_males"

# Perform FDR correction for 1137 pvals(1038+39)
padjs_aVMCs_Reynolds_Mono<-as.data.frame(p.adjust(c(pvals_aVMCs_Xfemales_Reynolds_Mono$pval_females,pvals_aVMCs_Xmales_Reynolds_Mono$pval_males),method="BH"))
colnames(padjs_aVMCs_Reynolds_Mono)<-"bacon_p.adj_disp"

# check females (how many of 1098) and males (how many of 39) separately how many were FDR significant 
# 1035 female aVMCs are still FDR significant in Reynolds Monocytes data.
padjs_aVMCs_females_Reynolds_Mono<-as.data.frame(padjs_aVMCs_Reynolds_Mono[1:1098,])
colnames(padjs_aVMCs_females_Reynolds_Mono)<-"bacon_p.adj_disp"
df_aVMCs_females_Reynolds_Mono$bacon_p.adj_disp<-padjs_aVMCs_females_Reynolds_Mono$bacon_p.adj_disp

df_aVMCs_females_Reynolds_Mono<-filter(df_aVMCs_females_Reynolds_Mono,bacon_p.adj_disp<0.05)
df_aVMCs_females_effectsize_up_Reynolds<-filter(df_aVMCs_females_Reynolds_Mono,bacon_effectsize_disp>0) # 1035
df_aVMCs_females_effectsize_down_Reynolds<-filter(df_aVMCs_females_Reynolds_Mono,bacon_effectsize_disp<0) # 0

# 39 male aVMCs are still FDR significant in Reynolds blood data.
padjs_aVMCs_males_Reynolds_Mono<-as.data.frame(padjs_aVMCs_Reynolds_Mono[1099:1137,])
colnames(padjs_aVMCs_males_Reynolds_Mono)<-"bacon_p.adj_disp"
df_aVMCs_males_Reynolds_Mono$bacon_p.adj_disp<-padjs_aVMCs_males_Reynolds_Mono$bacon_p.adj_disp

df_aVMCs_males_Reynolds_Mono<-filter(df_aVMCs_males_Reynolds_Mono,bacon_p.adj_disp<0.05)
df_aVMCs_males_effectsize_up_Reynolds<-filter(df_aVMCs_males_Reynolds_Mono,bacon_effectsize_disp>0) # 38
df_aVMCs_males_effectsize_down_Reynolds<-filter(df_aVMCs_males_Reynolds_Mono,bacon_effectsize_disp<0) # 1

save(df_aVMCs_females_effectsize_up_Reynolds,df_aVMCs_females_Reynolds_Mono,df_aVMCs_males_effectsize_up_Reynolds,df_aVMCs_males_effectsize_down_Reynolds,df_aVMCs_males_Reynolds_Mono,file = "Output_06_Reynolds_Catogue of aVMCs.RData")
```

```{r}
# Effect size of aVMCs between external datasets and BIOS.
df_Johansson_validation_females$tag<-"Johansson Blood"
df_Johansson_validation_females$sex<-"Females"
df_Johansson_validation_males$tag<-"Johansson Blood"
df_Johansson_validation_males$sex<-"Males"

df_Reynolds_validation_females$tag<-"Reynolds Monocytes"
df_Reynolds_validation_females$sex<-"Females"
df_Reynolds_validation_males$tag<-"Reynolds Monocytes"
df_Reynolds_validation_males$sex<-"Males"

# Combined 4 datasets.
colnames(df_Johansson_validation_females)[1]<-"External.datasets"
colnames(df_Johansson_validation_males)[1]<-"External.datasets" 

colnames(df_Reynolds_validation_females)[1]<-"External.datasets"
colnames(df_Reynolds_validation_males)[1]<-"External.datasets"

df_aVMCs_validation<-rbind(df_Johansson_validation_females,df_Johansson_validation_males,df_Reynolds_validation_females,df_Reynolds_validation_males)
df_aVMCs_validation$tag<-factor(df_aVMCs_validation$tag,levels = c("Johansson Blood","Reynolds Monocytes"))
df_aVMCs_validation$sex<-factor(df_aVMCs_validation$sex,levels = c("Males","Females"))
df_aVMCs_validation$aVMCs<-factor(df_aVMCs_validation$aVMCs,levels = c("female-specific","male-specific","both-sex"))

tiff("Output_06_Validation of aVMCs.tiff", units="in", width=16, height=13, res=600,compression = 'lzw')
ggplot(df_aVMCs_validation) +
    aes(x = BIOS,y = External.datasets,color=aVMCs) +
      facet_grid(tag~sex,switch ="y",scales = "fixed") +
     geom_point(size = 3, alpha = 0.5) +
     geom_vline(aes(xintercept=0)) +
     geom_hline(aes(yintercept=0)) +
     geom_abline(slope=1, intercept = 0)+
  scale_x_continuous(limits=c(-0.04, 0.04))+
  scale_y_continuous(limits=c(-0.04, 0.04))+
  scale_color_manual(values = c("#F08080", "#56B4E9","#32CD32"))+
  theme_bw()+
  theme(panel.spacing.x = unit(.5, "cm"),legend.position = "bottom",legend.title = element_text(size = 34),legend.text = element_text(size = 34),axis.title.y = element_text(size = 34),axis.title.x = element_text(size = 34),plot.title = element_text(hjust = 0.5,size = 34),axis.text.x = element_text(size = 34,angle = 45,hjust = 1),axis.text.y = element_text(size = 34),strip.text.x = element_text(size = 28),strip.text.y = element_text(size = 28))+
    labs(x="Effect size in BIOS",y="Effect size")
dev.off()
```


```{r}
# The overlap aVMCs between BIOS data and three external datasets.
# Females:925
Overlap_females<-Reduce(intersect, list(rownames(dat_DGLM_Xfemales_bacon_disp),rownames(df_aVMCs_females_Johansson_Blood),rownames(df_aVMCs_females_Reynolds_Mono)))
dat_aVMCs_females_validated<-dat_DGLM_Xfemales_bacon_disp[Overlap_females,]

# Males: 37
Overlap_males<-Reduce(intersect, list(rownames(dat_DGLM_Xmales_bacon_disp),rownames(df_aVMCs_males_Johansson_Blood),rownames(df_aVMCs_males_Reynolds_Mono)))
dat_aVMCs_males_validated<-dat_DGLM_Xmales_bacon_disp[Overlap_males,]

# Statistically significant aVMCs both in females and males.
both_sex<-intersect(rownames(dat_aVMCs_females_validated),rownames(dat_aVMCs_males_validated))
df_aVMCs_both_sex_validated<-dat_DGLM_Xfemales_bacon_disp[both_sex,]

# Statistically significant aVMCs only in females.
df_aVMCs_female_specific_validated<-subset(dat_aVMCs_females_validated,!rownames(dat_aVMCs_females_validated)%in% both_sex)

# Statistically significant aVMCs only in males.
df_aVMCs_male_specific_validated<-subset(dat_aVMCs_males_validated,!rownames(dat_aVMCs_males_validated)%in% both_sex)

save(df_aVMCs_both_sex_validated,df_aVMCs_female_specific_validated,df_aVMCs_male_specific_validated,dat_aVMCs_females_validated,dat_aVMCs_males_validated,file = "Output_06_Catalogue_aVMCs_validated.RData")
```