Script_05_aDMCs_validation.Rmd

---
title: "aDMCs 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 all necessary libraries.
library(minfi)
library(DNAmArray)
library(tidyverse)
library(sva)
library(bacon)
```


```{r}
# Load data
load(file = "Output_01_chrX_betas.RData")
load(file = "Output_02_aDMCs_effectsize_BIOS.RData")
load(file = "Input_05_Johansson_Blood_450k.RData")
load(file = "Input_05_Reynolds_Monocytes_450k.RData")
```

```{r}
# X-methylation profile of females and males using discovery and replication cohort 
# BIOS males
dat_mean_BIOS_Xmales<-as.data.frame(rowMeans(assay(betas_Xmales),na.rm = TRUE))
colnames(dat_mean_BIOS_Xmales)<-"BIOS.males"

# BIOS females
dat_mean_BIOS_Xfemales<-as.data.frame(rowMeans(assay(betas_Xfemales),na.rm = TRUE))
colnames(dat_mean_BIOS_Xfemales)<-"BIOS.females"

# Johansson females
dat_mean_Johansson_Blood_Xfemales<-as.data.frame(rowMeans(betas_Johansson_Blood_Xfemales,na.rm = TRUE))
colnames(dat_mean_Johansson_Blood_Xfemales)<-"Johansson.females"

# Johansson males
dat_mean_Johansson_Blood_Xmales<-as.data.frame(rowMeans(betas_Johansson_Blood_Xmales,na.rm = TRUE))
colnames(dat_mean_Johansson_Blood_Xmales)<-"Johansson.males"

# Reynolds females
dat_mean_Reynolds_Mono_Xfemales<-as.data.frame(rowMeans(betas_Reynolds_Mono_Xfemales,na.rm = TRUE))
colnames(dat_mean_Reynolds_Mono_Xfemales)<-"Reynolds.females"

# Reynolds males
dat_mean_Reynolds_Mono_Xmales<-as.data.frame(rowMeans(betas_Reynolds_Mono_Xmales,na.rm = TRUE))
colnames(dat_mean_Reynolds_Mono_Xmales)<-"Reynolds.males"

# Prepare mean X-methylation data in females and males.
dat_mean_BIOS<-cbind(dat_mean_BIOS_Xfemales,dat_mean_BIOS_Xmales)
colnames(dat_mean_BIOS)<-c("BIOS females","BIOS males")

dat_mean_Johansson_Blood<-cbind(dat_mean_Johansson_Blood_Xfemales,dat_mean_Johansson_Blood_Xmales)
colnames(dat_mean_Johansson_Blood)<-c("Johansson females","Johansson males")

dat_mean_Reynolds_Mono<-cbind(dat_mean_Reynolds_Mono_Xfemales,dat_mean_Reynolds_Mono_Xmales)
colnames(dat_mean_Reynolds_Mono)<-c("Reynolds females","Reynolds males")

save(dat_mean_BIOS,dat_mean_Johansson_Blood,dat_mean_Reynolds_Mono,file = "mean X-methylation profile.RData")

tiff("X methylation profile.tiff", units="in", width=8, height=6, res=600,compression = 'lzw')

par(lwd=4)
densityPlot(as.matrix(dat_mean_BIOS),add=T,main="X methylation profile",xlab="Beta values",sampGroups=c("BIOS Females","BIOS Males"),xlim=c(0,1),pal=c("#F08080","#87CEEB"),cex.main=1.5,cex.lab=1.5, cex.axis=1.5,legend = F) 

par(lwd=1)
densityPlot(as.matrix(dat_mean_Johansson_Blood),add=F,sampGroups=c("Johansson Females","Johansson Males"),pal=c("#F08080","#87CEEB"),legend = F)

par(lwd=1,lty=2)
densityPlot(as.matrix(dat_mean_Reynolds_Mono),add=F,sampGroups=c("Reynolds Females","Reynolds Males"),pal=c("#F08080","#87CEEB"),legend = F)
dev.off()
```


```{r}
#### aDMCs validation on Johansson blood datasets.
# rank-inverse normal(RIN) transformation for Xfemales and Xmales separately.
RIN <- function(x) {
  y <- rank(x, NA)
  y <- ppoints(y)[y]
  y <- qnorm(y)
  x[!is.na(x)] <- y
  x
 }
RIN.mvalues_Johansson_Blood_Xfemales <- t(apply(mvalues_Johansson_Blood_Xfemales, 1, RIN))
RIN.mvalues_Johansson_Blood_Xmales <- t(apply(mvalues_Johansson_Blood_Xmales, 1, RIN))
```


```{r}
## Prepare SVA
# Johansson Xfemales
## null model:only exclude interest variable
design0_Johansson_Blood_Xfemales = model.matrix(~ CD8T + CD4T + NK + Bcell + Mono + sentrix_position,data=Covariates_Johansson_Blood_Xfemales)
## Full model:all covarites
design_Johansson_Blood_Xfemales <- model.matrix(~Age + CD8T + CD4T + NK + Bcell + Mono + sentrix_position, data=Covariates_Johansson_Blood_Xfemales)

## Estimate Johansson female latent factors by SVA
svobj_Johansson_Blood_Xfemales = sva(RIN.mvalues_Johansson_Blood_Xfemales,design_Johansson_Blood_Xfemales,design0_Johansson_Blood_Xfemales,n.sv=5)
sv_Johansson_Blood_Xfemales <- as.data.frame(svobj_Johansson_Blood_Xfemales$sv)
colnames(sv_Johansson_Blood_Xfemales) <- c('SV1', 'SV2', 'SV3', 'SV4', 'SV5')

# Johansson Xmales
## null model:only exclude interest variable
design0_Johansson_Blood_Xmales = model.matrix(~ CD8T + CD4T + NK + Bcell + Mono + sentrix_position,data=Covariates_Johansson_Blood_Xmales)
## Full model:all covarites
design_Johansson_Blood_Xmales <- model.matrix(~Age + CD8T + CD4T + NK + Bcell + Mono + sentrix_position, data=Covariates_Johansson_Blood_Xmales)

# Estimate Johansson male latent factors by SVA
svobj_Johansson_Blood_Xmales = sva(RIN.mvalues_Johansson_Blood_Xmales,design_Johansson_Blood_Xmales,design0_Johansson_Blood_Xmales,n.sv=5)
sv_Johansson_Blood_Xmales <- as.data.frame(svobj_Johansson_Blood_Xmales$sv)
colnames(sv_Johansson_Blood_Xmales) <- c('SV1', 'SV2', 'SV3', 'SV4', 'SV5')
```

```{r}
# Add estimated latent factors to covariates dataframe. 
#Xfemales
Covariates_sv_Johansson_Blood_Xfemales<-cbind(Covariates_Johansson_Blood_Xfemales,sv_Johansson_Blood_Xfemales)

#Xmales
Covariates_sv_Johansson_Blood_Xmales<-cbind(Covariates_Johansson_Blood_Xmales,sv_Johansson_Blood_Xmales)

# Change sentrixt position into factor
Covariates_sv_Johansson_Blood_Xfemales$sentrix_position<-factor(Covariates_sv_Johansson_Blood_Xfemales$sentrix_position)
Covariates_sv_Johansson_Blood_Xmales$sentrix_position<-factor(Covariates_sv_Johansson_Blood_Xmales$sentrix_position)
```

```{r}
#### Construct a Function called runDGLM for Xfemales and Xmales separately.
## Identify aDMCs in Johansson Blood females.
runDGLM_Johansson_Blood_Xfemales <- function(x)
{
  dat_Xfemales_dglm <- matrix(nrow=nrow(x), ncol=9)
  colnames(dat_Xfemales_dglm)<-c("cgID","effectsize.linear","std.error.linear","t.linear","p.linear","effectsize.disp","std.error.disp","t.disp","p.disp")
  for(CpG in 1:nrow(x))
  {
    print(CpG) #to check progress
    probe <- rownames(x)[CpG]
    tempdata <- data.frame(Age=Covariates_sv_Johansson_Blood_Xfemales[,2],data = x[probe,],sentrix_position=Covariates_sv_Johansson_Blood_Xfemales[,4],CD8T=Covariates_sv_Johansson_Blood_Xfemales[,5], CD4T=Covariates_sv_Johansson_Blood_Xfemales[,6], NK=Covariates_sv_Johansson_Blood_Xfemales[,7], Bcell=Covariates_sv_Johansson_Blood_Xfemales[,8], Mono=Covariates_sv_Johansson_Blood_Xfemales[,9],  SV1=Covariates_sv_Johansson_Blood_Xfemales[,11], SV2=Covariates_sv_Johansson_Blood_Xfemales[,12], SV3=Covariates_sv_Johansson_Blood_Xfemales[,13], SV4=Covariates_sv_Johansson_Blood_Xfemales[,14], SV5=Covariates_sv_Johansson_Blood_Xfemales[,15])
    
    fit_Xfemales <- tryCatch({ dglm(formula = data ~ Age + CD8T + CD4T + NK + Bcell + Mono + sentrix_position + SV1 + SV2 + SV3 + SV4 + SV5, dformula = ~ Age, data=tempdata, family = gaussian(link = identity))},error=identity)
    if(is.null(fit_Xfemales$message)){ #good
      #extract info about the linear model
      effectsize.linear <- summary(fit_Xfemales)$coefficients[2,1]
      std.error.linear<-summary(fit_Xfemales)$coefficients[2,2]
      t.value.linear <-summary(fit_Xfemales)$coefficients[2,3]
      p.value.linear <- summary(fit_Xfemales)$coefficients[2,4]
      
      #extract info about the dispersion model
      effectsize.disp <- summary(fit_Xfemales$dispersion.fit)$coefficients[2]
      std.error.disp<-summary(fit_Xfemales$dispersion.fit)$coefficients[4]
      t.value.disp <- summary(fit_Xfemales$dispersion.fit)$coefficients[6]
      p.value.disp <- summary(fit_Xfemales$dispersion.fit)$coefficients[8]
      
      # Summary statistics Output
      out <- matrix(c(probe,effectsize.linear,std.error.linear,t.value.linear,p.value.linear,effectsize.disp,std.error.disp,t.value.disp,p.value.disp), ncol=9)
      dat_Xfemales_dglm[CpG,] <- out
    }else{                  #bad
      print("Didn't converge")
      out <- matrix(c(probe, fit_Xfemales$message, rep(NA, times = 7)), ncol=9)
      dat_Xfemales_dglm[CpG,] <- out
    }
    
  }
  return(dat_Xfemales_dglm)
}

# Run runDGLM function
DGLM_Johansson_Blood_Xfemales<- runDGLM_Johansson_Blood_Xfemales(RIN.mvalues_Johansson_Blood_Xfemales)
DGLM_Johansson_Blood_Xfemales<-data.frame(DGLM_Johansson_Blood_Xfemales) 
rownames(DGLM_Johansson_Blood_Xfemales)<-make.names(DGLM_Johansson_Blood_Xfemales$cgID,unique=TRUE)
DGLM_Johansson_Blood_Xfemales<-DGLM_Johansson_Blood_Xfemales[,2:9]
DGLM_Johansson_Blood_Xfemales_errors<-subset(DGLM_Johansson_Blood_Xfemales,is.na(DGLM_Johansson_Blood_Xfemales$t.linear)) # No errors(all CpGs converge)
DGLM_Johansson_Blood_Xfemales_converge<-DGLM_Johansson_Blood_Xfemales

# Convert each column to numerical variable.
DGLM_Johansson_Blood_Xfemales_converge<-as.data.frame(lapply(DGLM_Johansson_Blood_Xfemales_converge,as.numeric)) 
rownames(DGLM_Johansson_Blood_Xfemales_converge)<-make.names(rownames(DGLM_Johansson_Blood_Xfemales),unique=TRUE) 
save(DGLM_Johansson_Blood_Xfemales_converge,file = "Output_05_DGLM_Johansson_Blood_Xfemales_converge.RData")
```

```{r}
## Correct bias and inflation for linear model.
# Xfemales
set.seed(1)
bc_DGLM_Johansson_Blood_Xfemales<- bacon(DGLM_Johansson_Blood_Xfemales_converge[,3])

# Calculate bias and inflation
bias(bc_DGLM_Johansson_Blood_Xfemales)         # 0.12  
inflation(bc_DGLM_Johansson_Blood_Xfemales)    # 1.8

# Save histogram plot of test statistics 
tiff("Output_05_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="aDMCs 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[,1],DGLM_Johansson_Blood_Xfemales_converge[,2])
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[,1],DGLM_Johansson_Blood_Xfemales_converge[,2])
se_DGLM_Johansson_Blood_Xfemales<-se(bc_DGLM_Johansson_Blood_Xfemales_se)

# Summary Bacon statistics
dat_DGLM_Johansson_Blood_Xfemales_bacon_linear<-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_linear)<-c("bacon_effectsize_linear","bacon_std.error_linear","bacon_tstatistics_linear","bacon_pval_linear")
rownames(dat_DGLM_Johansson_Blood_Xfemales_bacon_linear)<-make.names(rownames(DGLM_Johansson_Blood_Xfemales_converge),unique = TRUE)  

# check females effecsize of aDMCs between BIOS and Johansson 
dat_Johansson_Blood_Xfemales_aDMCs<-dat_DGLM_Johansson_Blood_Xfemales_bacon_linear[rownames(df_aDMCs_effectsize_BIOS),]
df_Johansson_validation_females<-cbind(dat_Johansson_Blood_Xfemales_aDMCs[,1],df_aDMCs_effectsize_BIOS[,-2])
colnames(df_Johansson_validation_females)<-c("Johansson","BIOS","aDMCs")
```

```{r}
#### Construct a Function called runDGLM for Xfemales and Xmales separately.
## Identify aDMCs in Johansson Blood males.
runDGLM_Johansson_Blood_Xmales <- function(x)
{
  dat_Xmales_dglm <- matrix(nrow=nrow(x), ncol=9)
  colnames(dat_Xmales_dglm)<-c("cgID","effectsize.linear","std.error.linear","t.linear","p.linear","effectsize.disp","std.error.disp","t.disp","p.disp")
  for(CpG in 1:nrow(x))
  {
    print(CpG) #to check progress
    probe <- rownames(x)[CpG]
    tempdata <- data.frame(Age=Covariates_sv_Johansson_Blood_Xmales[,2],data = x[probe,],sentrix_position=Covariates_sv_Johansson_Blood_Xmales[,4],CD8T=Covariates_sv_Johansson_Blood_Xmales[,5], CD4T=Covariates_sv_Johansson_Blood_Xmales[,6], NK=Covariates_sv_Johansson_Blood_Xmales[,7], Bcell=Covariates_sv_Johansson_Blood_Xmales[,8], Mono=Covariates_sv_Johansson_Blood_Xmales[,9],  SV1=Covariates_sv_Johansson_Blood_Xmales[,11], SV2=Covariates_sv_Johansson_Blood_Xmales[,12], SV3=Covariates_sv_Johansson_Blood_Xmales[,13], SV4=Covariates_sv_Johansson_Blood_Xmales[,14], SV5=Covariates_sv_Johansson_Blood_Xmales[,15])
    
    fit_Xmales <- tryCatch({ dglm(formula = data ~ Age + CD8T + CD4T + NK + Bcell + Mono + sentrix_position + SV1 + SV2 + SV3 + SV4 + SV5, dformula = ~ Age, data=tempdata, family = gaussian(link = identity))},error=identity)
    if(is.null(fit_Xmales$message)){ #good
      #extract info about the linear model
      effectsize.linear <- summary(fit_Xmales)$coefficients[2,1]
      std.error.linear<-summary(fit_Xmales)$coefficients[2,2]
      t.value.linear <-summary(fit_Xmales)$coefficients[2,3]
      p.value.linear <- summary(fit_Xmales)$coefficients[2,4]
      
      #extract info about the dispersion model
      effectsize.disp <- summary(fit_Xmales$dispersion.fit)$coefficients[2]
      std.error.disp<-summary(fit_Xmales$dispersion.fit)$coefficients[4]
      t.value.disp <- summary(fit_Xmales$dispersion.fit)$coefficients[6]
      p.value.disp <- summary(fit_Xmales$dispersion.fit)$coefficients[8]
      
      # Summary statistics Output
      out <- matrix(c(probe,effectsize.linear,std.error.linear,t.value.linear,p.value.linear,effectsize.disp,std.error.disp,t.value.disp,p.value.disp), ncol=9)
      dat_Xmales_dglm[CpG,] <- out
    }else{                  #bad
      print("Didn't converge")
      out <- matrix(c(probe, fit_Xmales$message, rep(NA, times = 7)), ncol=9)
      dat_Xmales_dglm[CpG,] <- out
    }
    
  }
  return(dat_Xmales_dglm)
}

# Run runDGLM function
DGLM_Johansson_Blood_Xmales<- runDGLM_Johansson_Blood_Xmales(RIN.mvalues_Johansson_Blood_Xmales)
DGLM_Johansson_Blood_Xmales<-data.frame(DGLM_Johansson_Blood_Xmales) 
rownames(DGLM_Johansson_Blood_Xmales)<-make.names(DGLM_Johansson_Blood_Xmales$cgID,unique=TRUE)
DGLM_Johansson_Blood_Xmales<-DGLM_Johansson_Blood_Xmales[,2:9]
DGLM_Johansson_Blood_Xmales_errors<-subset(DGLM_Johansson_Blood_Xmales,is.na(DGLM_Johansson_Blood_Xmales$t.linear)) # No errors(all CpGs converge)
DGLM_Johansson_Blood_Xmales_converge<-DGLM_Johansson_Blood_Xmales

# Convert each column to numerical variable.
DGLM_Johansson_Blood_Xmales_converge<-as.data.frame(lapply(DGLM_Johansson_Blood_Xmales_converge,as.numeric)) 
rownames(DGLM_Johansson_Blood_Xmales_converge)<-make.names(rownames(DGLM_Johansson_Blood_Xmales),unique=TRUE)
save(DGLM_Johansson_Blood_Xmales_converge,file = "Output_05_DGLM_Johansson_Blood_Xmales_converge.RData")
```

```{r}
## Correct bias and inflation for linear model.
# Xmales
set.seed(1)
bc_DGLM_Johansson_Blood_Xmales<- bacon(DGLM_Johansson_Blood_Xmales_converge[,3])

# Calculate bias and inflation
bias(bc_DGLM_Johansson_Blood_Xmales)         # -0.2  
inflation(bc_DGLM_Johansson_Blood_Xmales)    # 1.6

# Save histogram plot of test statistics 
tiff("Output_05_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="aDMCs 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[,1],DGLM_Johansson_Blood_Xmales_converge[,2])
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[,1],DGLM_Johansson_Blood_Xmales_converge[,2])
se_DGLM_Johansson_Blood_Xmales<-se(bc_DGLM_Johansson_Blood_Xmales_se)

# Summary Bacon statistics
dat_DGLM_Johansson_Blood_Xmales_bacon_linear<-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_linear)<-c("bacon_effectsize_linear","bacon_std.error_linear","bacon_tstatistics_linear","bacon_pval_linear")
rownames(dat_DGLM_Johansson_Blood_Xmales_bacon_linear)<-make.names(rownames(DGLM_Johansson_Blood_Xmales_converge),unique = TRUE)

# check males effecsize of aDMCs between BIOS and Hannum 
dat_Johansson_Blood_Xmales_aDMCs<-dat_DGLM_Johansson_Blood_Xmales_bacon_linear[rownames(df_aDMCs_effectsize_BIOS),]
df_Johansson_validation_males<-cbind(dat_Johansson_Blood_Xmales_aDMCs[,1],df_aDMCs_effectsize_BIOS[,-1])
colnames(df_Johansson_validation_males)<-c("Johansson","BIOS","aDMCs")
```

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

```{r}
# Extract 80(33+47) female p-values and 1837(1790+47) male p-values.
dat_Johansson_Blood_Xfemales_aDMCs$aDMCs<-df_aDMCs_effectsize_BIOS$aDMCs
df_aDMCs_females_Johansson_Blood<-filter(dat_Johansson_Blood_Xfemales_aDMCs,aDMCs=="female-specific"|aDMCs=="both-sex")
pvals_aDMCs_Xfemales_Johansson_Blood<-as.data.frame(pvals_DGLM_Johansson_Blood_Xfemales[rownames(df_aDMCs_females_Johansson_Blood),])
colnames(pvals_aDMCs_Xfemales_Johansson_Blood)<-"pval_females"

dat_Johansson_Blood_Xmales_aDMCs$aDMCs<-df_aDMCs_effectsize_BIOS$aDMCs
df_aDMCs_males_Johansson_Blood<-filter(dat_Johansson_Blood_Xmales_aDMCs,aDMCs=="male-specific"|aDMCs=="both-sex")
pvals_aDMCs_Xmales_Johansson_Blood<-as.data.frame(pvals_DGLM_Johansson_Blood_Xmales[rownames(df_aDMCs_males_Johansson_Blood),])
colnames(pvals_aDMCs_Xmales_Johansson_Blood)<-"pval_males"

# Perform FDR correction for 1588 pvals(80+1837)
padjs_aDMCs_Johansson_Blood<-as.data.frame(p.adjust(c(pvals_aDMCs_Xfemales_Johansson_Blood$pval_females,pvals_aDMCs_Xmales_Johansson_Blood$pval_males),method="BH"))
colnames(padjs_aDMCs_Johansson_Blood)<-"bacon_p.adj_linear"

# check males (how many of 1837) and females (how many of 80) separately how many were FDR significant 
# 79 female aDMCs are still FDR significant in Hannum blood data.
padjs_aDMCs_females_Johansson_Blood<-as.data.frame(padjs_aDMCs_Johansson_Blood[1:80,])
colnames(padjs_aDMCs_females_Johansson_Blood)<-"bacon_p.adj_linear"
df_aDMCs_females_Johansson_Blood$bacon_p.adj_linear<-padjs_aDMCs_females_Johansson_Blood$bacon_p.adj_linear

df_aDMCs_females_Johansson_Blood<-filter(df_aDMCs_females_Johansson_Blood,bacon_p.adj_linear<0.05)
df_aDMCs_females_effectsize_up_Johansson<-filter(df_aDMCs_females_Johansson_Blood,bacon_effectsize_linear>0) # 38
df_aDMCs_females_effectsize_down_Johansson<-filter(df_aDMCs_females_Johansson_Blood,bacon_effectsize_linear<0) # 41

# 1371 male aDMCs are still FDR significant in Johansson blood data.
padjs_aDMCs_males_Johansson_Blood<-as.data.frame(padjs_aDMCs_Johansson_Blood[81:1917,])
colnames(padjs_aDMCs_males_Johansson_Blood)<-"bacon_p.adj_linear"
df_aDMCs_males_Johansson_Blood$bacon_p.adj_linear<-padjs_aDMCs_males_Johansson_Blood$bacon_p.adj_linear

df_aDMCs_males_Johansson_Blood<-filter(df_aDMCs_males_Johansson_Blood,bacon_p.adj_linear<0.05)
df_aDMCs_males_effectsize_up_Johansson<-filter(df_aDMCs_males_Johansson_Blood,bacon_effectsize_linear>0) # 839
df_aDMCs_males_effectsize_down_Johansson<-filter(df_aDMCs_males_Johansson_Blood,bacon_effectsize_linear<0) # 532

save(df_aDMCs_females_effectsize_up_Johansson,df_aDMCs_females_effectsize_down_Johansson,df_aDMCs_females_Johansson_Blood,df_aDMCs_males_effectsize_up_Johansson,df_aDMCs_males_effectsize_down_Johansson,df_aDMCs_males_Johansson_Blood,file = "Output_05_Johansson_Catogue of aDMCs.RData")
```

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

```{r}
## aDMCs validation on Reynolds Mono datasets.
# rank-inverse normal(RIN) transformation for Xfemales and Xmales separately.
RIN.mvalues_Reynolds_Mono_Xfemales <- t(apply(mvalues_Reynolds_Mono_Xfemales, 1, RIN))
RIN.mvalues_Reynolds_Mono_Xmales <- t(apply(mvalues_Reynolds_Mono_Xmales, 1, RIN))
```


```{r}
## Prepare SVA
# Reynolds Mono Xfemales
## null model:only exclude interest variable
design0_Reynolds_Mono_Xfemales = model.matrix(~ Bcell + Tcell + NK + Neutro + sentrix_position,data=Covariates_Reynolds_Mono_Xfemales)
## Full model:all covarites
design_Reynolds_Mono_Xfemales <- model.matrix(~Age + Bcell + Tcell + NK + Neutro + sentrix_position, data=Covariates_Reynolds_Mono_Xfemales)

## Estimate Reynolds female latent factors by SVA
svobj_Reynolds_Mono_Xfemales = sva(RIN.mvalues_Reynolds_Mono_Xfemales,design_Reynolds_Mono_Xfemales,design0_Reynolds_Mono_Xfemales,n.sv=5)
sv_Reynolds_Mono_Xfemales <- as.data.frame(svobj_Reynolds_Mono_Xfemales$sv)
colnames(sv_Reynolds_Mono_Xfemales) <- c('SV1', 'SV2', 'SV3', 'SV4', 'SV5')

# Reynolds Mono Xmales
## null model:only exclude interest variable
design0_Reynolds_Mono_Xmales = model.matrix(~ Bcell + Tcell + NK + Neutro + sentrix_position,data=Covariates_Reynolds_Mono_Xmales)
## Full model:all covarites
design_Reynolds_Mono_Xmales <- model.matrix(~Age + Bcell + Tcell + NK + Neutro + sentrix_position, data=Covariates_Reynolds_Mono_Xmales)

# Estimate Reynolds male latent factors by SVA
svobj_Reynolds_Mono_Xmales = sva(RIN.mvalues_Reynolds_Mono_Xmales,design_Reynolds_Mono_Xmales,design0_Reynolds_Mono_Xmales,n.sv=5)
sv_Reynolds_Mono_Xmales <- as.data.frame(svobj_Reynolds_Mono_Xmales$sv)
colnames(sv_Reynolds_Mono_Xmales) <- c('SV1', 'SV2', 'SV3', 'SV4', 'SV5')
```

```{r}
# Add estimated latent factors to covariates dataframe. 
#Xfemales
Covariates_sv_Reynolds_Mono_Xfemales<-cbind(Covariates_Reynolds_Mono_Xfemales,sv_Reynolds_Mono_Xfemales)

#Xmales
Covariates_sv_Reynolds_Mono_Xmales<-cbind(Covariates_Reynolds_Mono_Xmales,sv_Reynolds_Mono_Xmales)

# Change sentrixt position into factor
Covariates_sv_Reynolds_Mono_Xfemales$sentrix_position<-factor(Covariates_sv_Reynolds_Mono_Xfemales$sentrix_position)
Covariates_sv_Reynolds_Mono_Xmales$sentrix_position<-factor(Covariates_sv_Reynolds_Mono_Xmales$sentrix_position)
```

```{r}
#### Construct a Function called runDGLM for Xfemales and Xmales separately.
## Identify aDMCs in Reynolds Mono females.
runDGLM_Reynolds_Mono_Xfemales <- function(x)
{
  dat_Xfemales_dglm <- matrix(nrow=nrow(x), ncol=9)
  colnames(dat_Xfemales_dglm)<-c("cgID","effectsize.linear","std.error.linear","t.linear","p.linear","effectsize.disp","std.error.disp","t.disp","p.disp")
  for(CpG in 1:nrow(x))
  {
    print(CpG) #to check progress
    probe <- rownames(x)[CpG]
    tempdata <- data.frame(Age=Covariates_sv_Reynolds_Mono_Xfemales[,2],data = x[probe,],sentrix_position=Covariates_sv_Reynolds_Mono_Xfemales[,3],Bcell=Covariates_sv_Reynolds_Mono_Xfemales[,4],Tcell=Covariates_sv_Reynolds_Mono_Xfemales[,5], NK=Covariates_sv_Reynolds_Mono_Xfemales[,6], Neutro=Covariates_sv_Reynolds_Mono_Xfemales[,7],   SV1=Covariates_sv_Reynolds_Mono_Xfemales[,9], SV2=Covariates_sv_Reynolds_Mono_Xfemales[,10], SV3=Covariates_sv_Reynolds_Mono_Xfemales[,11], SV4=Covariates_sv_Reynolds_Mono_Xfemales[,12], SV5=Covariates_sv_Reynolds_Mono_Xfemales[,13])
    
    fit_Xfemales <- tryCatch({ dglm(formula = data ~ Age + Bcell + Tcell + NK + Neutro + sentrix_position + SV1 + SV2 + SV3 + SV4 + SV5, dformula = ~ Age, data=tempdata, family = gaussian(link = identity))},error=identity)
    if(is.null(fit_Xfemales$message)){ #good
      #extract info about the linear model
      effectsize.linear <- summary(fit_Xfemales)$coefficients[2,1]
      std.error.linear<-summary(fit_Xfemales)$coefficients[2,2]
      t.value.linear <-summary(fit_Xfemales)$coefficients[2,3]
      p.value.linear <- summary(fit_Xfemales)$coefficients[2,4]
      
      #extract info about the dispersion model
      effectsize.disp <- summary(fit_Xfemales$dispersion.fit)$coefficients[2]
      std.error.disp<-summary(fit_Xfemales$dispersion.fit)$coefficients[4]
      t.value.disp <- summary(fit_Xfemales$dispersion.fit)$coefficients[6]
      p.value.disp <- summary(fit_Xfemales$dispersion.fit)$coefficients[8]
      
      # Summary statistics Output
      out <- matrix(c(probe,effectsize.linear,std.error.linear,t.value.linear,p.value.linear,effectsize.disp,std.error.disp,t.value.disp,p.value.disp), ncol=9)
      dat_Xfemales_dglm[CpG,] <- out
    }else{                  #bad
      print("Didn't converge")
      out <- matrix(c(probe, fit_Xfemales$message, rep(NA, times = 7)), ncol=9)
      dat_Xfemales_dglm[CpG,] <- out
    }
    
  }
  return(dat_Xfemales_dglm)
}

# Run runDGLM function
DGLM_Reynolds_Mono_Xfemales<- runDGLM_Reynolds_Mono_Xfemales(RIN.mvalues_Reynolds_Mono_Xfemales)
DGLM_Reynolds_Mono_Xfemales<-data.frame(DGLM_Reynolds_Mono_Xfemales) 
rownames(DGLM_Reynolds_Mono_Xfemales)<-make.names(DGLM_Reynolds_Mono_Xfemales$cgID,unique=TRUE)
DGLM_Reynolds_Mono_Xfemales<-DGLM_Reynolds_Mono_Xfemales[,2:9]
DGLM_Reynolds_Mono_Xfemales_errors<-subset(DGLM_Reynolds_Mono_Xfemales,is.na(DGLM_Reynolds_Mono_Xfemales$t.linear)) # No errors(all CpGs converge)
DGLM_Reynolds_Mono_Xfemales_converge<-DGLM_Reynolds_Mono_Xfemales

# Convert each column to numerical variable.
DGLM_Reynolds_Mono_Xfemales_converge<-as.data.frame(lapply(DGLM_Reynolds_Mono_Xfemales_converge,as.numeric)) 
rownames(DGLM_Reynolds_Mono_Xfemales_converge)<-make.names(rownames(DGLM_Reynolds_Mono_Xfemales),unique=TRUE) 
save(DGLM_Reynolds_Mono_Xfemales_converge,file = "Output_05_DGLM_Reynolds_Mono_Xfemales_converge.RData")
```

```{r}
## Correct bias and inflation for linear model.
# Xfemales
set.seed(1)
bc_DGLM_Reynolds_Mono_Xfemales<- bacon(DGLM_Reynolds_Mono_Xfemales_converge[,3])

# Calculate bias and inflation
bias(bc_DGLM_Reynolds_Mono_Xfemales)         # -0.06  
inflation(bc_DGLM_Reynolds_Mono_Xfemales)    # 1.2

# Save histogram plot of test statistics 
tiff("Output_05_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="aDMCs in Reynolds Mono 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[,1],DGLM_Reynolds_Mono_Xfemales_converge[,2])
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[,1],DGLM_Reynolds_Mono_Xfemales_converge[,2])
se_DGLM_Reynolds_Mono_Xfemales<-se(bc_DGLM_Reynolds_Mono_Xfemales_se)

# Summary Bacon statistics
dat_DGLM_Reynolds_Mono_Xfemales_bacon_linear<-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_linear)<-c("bacon_effectsize_linear","bacon_std.error_linear","bacon_tstatistics_linear","bacon_pval_linear")
rownames(dat_DGLM_Reynolds_Mono_Xfemales_bacon_linear)<-make.names(rownames(DGLM_Reynolds_Mono_Xfemales_converge),unique = TRUE)

# check females effecsize of aDMCs between BIOS and Reynolds Mono 
dat_Reynolds_Mono_Xfemales_aDMCs<-dat_DGLM_Reynolds_Mono_Xfemales_bacon_linear[rownames(df_aDMCs_effectsize_BIOS),]
df_Reynolds_validation_females<-cbind(dat_Reynolds_Mono_Xfemales_aDMCs[,1],df_aDMCs_effectsize_BIOS[,-2])
colnames(df_Reynolds_validation_females)<-c("Reynolds","BIOS","aDMCs")
```

```{r}
#### Construct a Function called runDGLM for Xfemales and Xmales separately.
## Identify aDMCs in Liu Mono males.
runDGLM_Reynolds_Mono_Xmales <- function(x)
{
  dat_Xmales_dglm <- matrix(nrow=nrow(x), ncol=9)
  colnames(dat_Xmales_dglm)<-c("cgID","effectsize.linear","std.error.linear","t.linear","p.linear","effectsize.disp","std.error.disp","t.disp","p.disp")
  for(CpG in 1:nrow(x))
  {
    print(CpG) #to check progress
    probe <- rownames(x)[CpG]
    tempdata <- data.frame(Age=Covariates_sv_Reynolds_Mono_Xmales[,2],data = x[probe,],sentrix_position=Covariates_sv_Reynolds_Mono_Xmales[,3],Bcell=Covariates_sv_Reynolds_Mono_Xmales[,4],Tcell=Covariates_sv_Reynolds_Mono_Xmales[,5], NK=Covariates_sv_Reynolds_Mono_Xmales[,6],Neutro=Covariates_sv_Reynolds_Mono_Xmales[,7],SV1=Covariates_sv_Reynolds_Mono_Xmales[,9],SV2=Covariates_sv_Reynolds_Mono_Xmales[,10],SV3=Covariates_sv_Reynolds_Mono_Xmales[,11],SV4=Covariates_sv_Reynolds_Mono_Xmales[,12],SV5=Covariates_sv_Reynolds_Mono_Xmales[,13])
    
    fit_Xmales <- tryCatch({ dglm(formula = data ~ Age + Bcell + Tcell + NK + Neutro + sentrix_position + SV1 + SV2 + SV3 + SV4 + SV5, dformula = ~ Age, data=tempdata, family = gaussian(link = identity))},error=identity)
    if(is.null(fit_Xmales$message)){ #good
      #extract info about the linear model
      effectsize.linear <- summary(fit_Xmales)$coefficients[2,1]
      std.error.linear<-summary(fit_Xmales)$coefficients[2,2]
      t.value.linear <-summary(fit_Xmales)$coefficients[2,3]
      p.value.linear <- summary(fit_Xmales)$coefficients[2,4]
      
      #extract info about the dispersion model
      effectsize.disp <- summary(fit_Xmales$dispersion.fit)$coefficients[2]
      std.error.disp<-summary(fit_Xmales$dispersion.fit)$coefficients[4]
      t.value.disp <- summary(fit_Xmales$dispersion.fit)$coefficients[6]
      p.value.disp <- summary(fit_Xmales$dispersion.fit)$coefficients[8]
      
      # Summary statistics Output
      out <- matrix(c(probe,effectsize.linear,std.error.linear,t.value.linear,p.value.linear,effectsize.disp,std.error.disp,t.value.disp,p.value.disp), ncol=9)
      dat_Xmales_dglm[CpG,] <- out
    }else{                  #bad
      print("Didn't converge")
      out <- matrix(c(probe, fit_Xmales$message, rep(NA, times = 7)), ncol=9)
      dat_Xmales_dglm[CpG,] <- out
    }
    
  }
  return(dat_Xmales_dglm)
}

# Run runDGLM function
DGLM_Reynolds_Mono_Xmales<- runDGLM_Reynolds_Mono_Xmales(RIN.mvalues_Reynolds_Mono_Xmales)
DGLM_Reynolds_Mono_Xmales<-data.frame(DGLM_Reynolds_Mono_Xmales) 
rownames(DGLM_Reynolds_Mono_Xmales)<-make.names(DGLM_Reynolds_Mono_Xmales$cgID,unique=TRUE)
DGLM_Reynolds_Mono_Xmales<-DGLM_Reynolds_Mono_Xmales[,2:9]
DGLM_Reynolds_Mono_Xmales_errors<-subset(DGLM_Reynolds_Mono_Xmales,is.na(DGLM_Reynolds_Mono_Xmales$t.linear)) # No errors(all CpGs converge)
DGLM_Reynolds_Mono_Xmales_converge<-DGLM_Reynolds_Mono_Xmales

# Convert each column to numerical variable.
DGLM_Reynolds_Mono_Xmales_converge<-as.data.frame(lapply(DGLM_Reynolds_Mono_Xmales_converge,as.numeric)) 
rownames(DGLM_Reynolds_Mono_Xmales_converge)<-make.names(rownames(DGLM_Reynolds_Mono_Xmales),unique=TRUE) 
save(DGLM_Reynolds_Mono_Xmales_converge,file = "Output_05_DGLM_Reynolds_Mono_Xmales_converge.RData")
```

```{r}
## Correct bias and inflation for linear model.
# Xmales
set.seed(1)
bc_DGLM_Reynolds_Mono_Xmales<- bacon(DGLM_Reynolds_Mono_Xmales_converge[,3])

# Calculate bias and inflation
bias(bc_DGLM_Reynolds_Mono_Xmales)         # -0.05  
inflation(bc_DGLM_Reynolds_Mono_Xmales)    # 1.15

# Save histogram plot of test statistics 
tiff("Output_05_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="aDMCs in Reynolds Mono 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[,1],DGLM_Reynolds_Mono_Xmales_converge[,2])
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[,1],DGLM_Reynolds_Mono_Xmales_converge[,2])
se_DGLM_Reynolds_Mono_Xmales<-se(bc_DGLM_Reynolds_Mono_Xmales_se)

# Summary Bacon statistics
dat_DGLM_Reynolds_Mono_Xmales_bacon_linear<-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_linear)<-c("bacon_effectsize_linear","bacon_std.error_linear","bacon_tstatistics_linear","bacon_pval_linear")
rownames(dat_DGLM_Reynolds_Mono_Xmales_bacon_linear)<-make.names(rownames(DGLM_Reynolds_Mono_Xmales_converge),unique = TRUE)

# check males effecsize of aDMCs between BIOS and Reynolds Mono 
dat_Reynolds_Mono_Xmales_aDMCs<-dat_DGLM_Reynolds_Mono_Xmales_bacon_linear[rownames(df_aDMCs_effectsize_BIOS),]
df_Reynolds_validation_males<-cbind(dat_Reynolds_Mono_Xmales_aDMCs[,1],df_aDMCs_effectsize_BIOS[,-1])
colnames(df_Reynolds_validation_males)<-c("Reynolds","BIOS","aDMCs")
```

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

```{r}
# Extract 80(33+47) female p-values and 1837(1790+47) male p-values.
dat_Reynolds_Mono_Xfemales_aDMCs$aDMCs<-df_aDMCs_effectsize_BIOS$aDMCs
df_aDMCs_females_Reynolds_Mono<-filter(dat_Reynolds_Mono_Xfemales_aDMCs,aDMCs=="female-specific"|aDMCs=="both-sex")
pvals_aDMCs_Xfemales_Reynolds_Mono<-as.data.frame(pvals_DGLM_Reynolds_Mono_Xfemales[rownames(df_aDMCs_females_Reynolds_Mono),])
colnames(pvals_aDMCs_Xfemales_Reynolds_Mono)<-"pval_females"

dat_Reynolds_Mono_Xmales_aDMCs$aDMCs<-df_aDMCs_effectsize_BIOS$aDMCs
df_aDMCs_males_Reynolds_Mono<-filter(dat_Reynolds_Mono_Xmales_aDMCs,aDMCs=="male-specific"|aDMCs=="both-sex")
pvals_aDMCs_Xmales_Reynolds_Mono<-as.data.frame(pvals_DGLM_Reynolds_Mono_Xmales[rownames(df_aDMCs_males_Reynolds_Mono),])
colnames(pvals_aDMCs_Xmales_Reynolds_Mono)<-"pval_males"

# Perform FDR correction for 1917 pvals(80+1837)
padjs_aDMCs_Reynolds_Mono<-as.data.frame(p.adjust(c(pvals_aDMCs_Xfemales_Reynolds_Mono$pval_females,pvals_aDMCs_Xmales_Reynolds_Mono$pval_males),method="BH"))
colnames(padjs_aDMCs_Reynolds_Mono)<-"bacon_p.adj_linear"

# check males (how many of 1837) and females (how many of 80) separately how many were FDR significant 
# 34 female aDMCs are still FDR significant in Reynolds Mono data.
padjs_aDMCs_females_Reynolds_Mono<-as.data.frame(padjs_aDMCs_Reynolds_Mono[1:80,])
colnames(padjs_aDMCs_females_Reynolds_Mono)<-"bacon_p.adj_linear"
df_aDMCs_females_Reynolds_Mono$bacon_p.adj_linear<-padjs_aDMCs_females_Reynolds_Mono$bacon_p.adj_linear

df_aDMCs_females_Reynolds_Mono<-filter(df_aDMCs_females_Reynolds_Mono,bacon_p.adj_linear<0.05)
df_aDMCs_females_effectsize_up_Reynolds<-filter(df_aDMCs_females_Reynolds_Mono,bacon_effectsize_linear>0) # 10
df_aDMCs_females_effectsize_down_Reynolds<-filter(df_aDMCs_females_Reynolds_Mono,bacon_effectsize_linear<0) # 24

# 340 male aDMCs are still FDR significant in Reynolds Mono data.
padjs_aDMCs_males_Reynolds_Mono<-as.data.frame(padjs_aDMCs_Reynolds_Mono[81:1917,])
colnames(padjs_aDMCs_males_Reynolds_Mono)<-"bacon_p.adj_linear"
df_aDMCs_males_Reynolds_Mono$bacon_p.adj_linear<-padjs_aDMCs_males_Reynolds_Mono$bacon_p.adj_linear

df_aDMCs_males_Reynolds_Mono<-filter(df_aDMCs_males_Reynolds_Mono,bacon_p.adj_linear<0.05)
df_aDMCs_males_effectsize_up_Reynolds<-filter(df_aDMCs_males_Reynolds_Mono,bacon_effectsize_linear>0) # 145
df_aDMCs_males_effectsize_down_Reynolds<-filter(df_aDMCs_males_Reynolds_Mono,bacon_effectsize_linear<0) # 195

save(df_aDMCs_females_effectsize_up_Reynolds,df_aDMCs_females_effectsize_down_Reynolds,df_aDMCs_females_Reynolds_Mono,df_aDMCs_males_effectsize_up_Reynolds,df_aDMCs_males_effectsize_down_Reynolds,df_aDMCs_males_Reynolds_Mono,file = "Output_05_Reynolds_Catogue of aDMCs.RData")
```

```{r}
# Effect size of aDMCs 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_aDMCs_validation<-rbind(df_Johansson_validation_females,df_Johansson_validation_males,df_Reynolds_validation_females,df_Reynolds_validation_males)
df_aDMCs_validation$tag<-factor(df_aDMCs_validation$tag,levels = c("Johansson Blood","Reynolds Monocytes"))
df_aDMCs_validation$sex<-factor(df_aDMCs_validation$sex,
                                levels = c("Males","Females"))
df_aDMCs_validation$aDMCs<-factor(df_aDMCs_validation$aDMCs,levels = c("female-specific","male-specific","both-sex"))

tiff("Output_05_Validation of aDMCs.tiff", units="in", width=16, height=13, res=600,compression = 'lzw')
ggplot(df_aDMCs_validation) +
    aes(x = BIOS,y = External.datasets,color=aDMCs) +
      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),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 aDMCs between BIOS data and two external datasets.
# Females:33
Overlap_females<-Reduce(intersect, list(rownames(dat_DGLM_Xfemales_bacon_linear),rownames(df_aDMCs_females_Johansson_Blood),rownames(df_aDMCs_females_Reynolds_Mono)))
dat_aDMCs_females_validated<-dat_DGLM_Xfemales_bacon_linear[Overlap_females,]

# Males: 316
Overlap_males<-Reduce(intersect, list(rownames(dat_DGLM_Xmales_bacon_linear),rownames(df_aDMCs_males_Johansson_Blood),rownames(df_aDMCs_males_Reynolds_Mono)))
dat_aDMCs_males_validated<-dat_DGLM_Xmales_bacon_linear[Overlap_males,]

# Statistically significant aDMCs both in females and males.
both_sex<-intersect(rownames(dat_aDMCs_females_validated),rownames(dat_aDMCs_males_validated))
df_aDMCs_both_sex_validated<-dat_DGLM_Xfemales_bacon_linear[both_sex,]

# Statistically significant aDMCs only in females.
df_aDMCs_female_specific_validated<-subset(dat_aDMCs_females_validated,!rownames(dat_aDMCs_females_validated)%in% both_sex)

# Statistically significant aDMCs only in males.
df_aDMCs_male_specific_validated<-subset(dat_aDMCs_males_validated,!rownames(dat_aDMCs_males_validated)%in% both_sex)

save(df_aDMCs_both_sex_validated,df_aDMCs_female_specific_validated,df_aDMCs_male_specific_validated,dat_aDMCs_females_validated,dat_aDMCs_males_validated,file = "Output_05_Catalogue_aDMCs_validated.RData")
```