Project Summary

This script contains the code for all figures and data visualizations.

##Data setup You will need, in the same directory as this file:

  1. A folder titled Manifest containing infinium-methylationepic-v-1-0-b5-manifest-file.csv
  2. A data file objectsforpreplot.RData, containing meta_m_final, betas.n.bs, and betas.n.ox, from EPIC_QC_MLE
  3. A data file postQC_mc_hmc_meta.RData, the output of EPIC_QC_MLE
  4. A data file containing the output from Chunk 14 in EPIC_QC_MLE, titled “data_chunk_14.Rdata”
  5. The data file containing the output fit_md from DiffMeth
  6. A CSV containing the results from annotation with ChromHMM

R Markdown

Density Plots

#Uses files meta_m_final and betas.n.bs and betas.n.ox, from EPIC_QC_MLE

load(file="objectsforpreplot.RData")

#Create meta files for data object filtering
meta_bs_c <- meta_m_final[which(meta_m_final$Assay == "1" & meta_m_final$Disease_coded == 1),]
meta_ox_c <- meta_m_final[which(meta_m_final$Assay == "2" & meta_m_final$Disease_coded == 1),]

#Create data objects for plotting
betas.bs.c<-betas.n.bs[,colnames(betas.n.bs) %in% meta_bs_c$Sample_Name]
betas.ox.c<-betas.n.ox[,colnames(betas.n.ox) %in% meta_ox_c$Sample_Name]

#Change shape of objects for plotting
samples.c.bs <- melt(betas.bs.c)
samples.c.bs$Var1<-NULL
colnames(samples.c.bs)[colnames(samples.c.bs) == "Var2"] <- "variable"
samples.c.ox <- melt(betas.ox.c)

samples.c.bs$assay <- "Control BS"
samples.c.ox$assay <- "Control oxBS"

# Combine raw beta values into one data object
raw.beta.plot <-rbind(samples.c.ox, samples.c.bs)

# Plot beta values by assay
plot.raw.beta <-ggplot(raw.beta.plot, aes(x=value)) +
                       geom_density(aes(group=assay, colour=assay)) +
                       scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
                       scale_y_continuous(limits = c(0, 5), expand = c(0,0)) +
                       theme(axis.line = element_line(size = 0.5, colour = "grey70"),
                       panel.grid = element_blank(), 
                       panel.background = element_rect(fill = "transparent"),
    plot.background = element_rect(fill='transparent', color=NA), 
          legend.background = element_rect(fill='transparent'), #transparent legend bg
          legend.box.background = element_rect(fill='transparent')) +
                       theme(legend.title = element_blank ()) +
                       xlab("beta") + ylab("density") +
                       stat_density(aes(x=value, colour=assay),
                       geom="line", position="identity") +
                       scale_color_brewer(palette = "Set1")

ggsave("densityplotpremle_c_only.png", plot.raw.beta, bg="transparent")

plot.raw.beta

#Uses MLE results (mc.f, hmc.f) and meta_m_final from EPIC_QC_MLE
load("./postQC_mc_hmc_meta.RData")


# Remove extra columns
mc.f$probe <- NULL
mc.f$mean_5mC <- NULL

# Remove extra columns
hmc.f$probe <- NULL
hmc.f$mean_5hmC <- NULL

meta_bs_c <- meta_m_final[which(meta_m_final$Assay == "1" & meta_m_final$Disease_coded == 1),]
meta_ox_c <- meta_m_final[which(meta_m_final$Assay == "2" & meta_m_final$Disease_coded == 1),]

#Create data objects for plotting
mc_C <- mc.f[,c(meta_bs_c$Sample_Name)]
hmc_C <- hmc.f[,c(meta_ox_c$Sample_Name)]

#Change shape of data objects for plotting
mle.c.mc <- melt(mc_C)
mle.c.hmc <- melt(hmc_C)

# Add assay identifier column for grouping in plot
mle.c.mc$assay <- "Control mC"
mle.c.hmc$assay <- "Control hmC"

# Combine raw beta values into one object for plotting
mle.plot <-rbind(mle.c.mc, mle.c.hmc)

# Plot beta values by assay
plot.mle <-ggplot(mle.plot, aes(x=value)) +
                       geom_density(aes(group=assay, colour=assay)) +
                       scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
                       scale_y_continuous(limits = c(0, 5), expand = c(0,0)) +
                       theme(axis.line = element_line(linewidth = 0.5, colour = "grey70"),
                       panel.grid = element_blank(), 
                       panel.background = element_rect(fill = "transparent"), plot.background=element_rect(fill="transparent"), 
          legend.background = element_rect(fill='transparent'), #transparent legend bg
          legend.box.background = element_rect(fill='transparent') )+
                       theme(legend.title = element_blank ()) +
                       xlab("beta") + ylab("density") +
                       stat_density(aes(x=value, colour=assay),
                       geom="line", position="identity") +
                       scale_color_brewer(palette = "Set1")

ggsave("mleplot_c_only.png", plot.mle, width=7.28, height=4.51, bg="transparent")

plot.mle

##Volcano plot

# Read in manifest (specify correct file for 450K or EPIC), remove headers, and
#label probe ID as "probe".
Illumina.EPIC.manifest  <- read.csv("./Manifest/infinium-methylationepic-v-1-0-b5-manifest-file.csv", skip = 7)
colnames(Illumina.EPIC.manifest)[colnames(Illumina.EPIC.manifest) == "Name"]<- "probe"

myvars <- c("probe", "CHR", "MAPINFO", "UCSC_RefGene_Name", "UCSC_RefGene_Accession", 
  "UCSC_RefGene_Group", "UCSC_CpG_Islands_Name", "Relation_to_UCSC_CpG_Island", 
  "Phantom4_Enhancers", "Phantom5_Enhancers", "DMR","X450k_Enhancer", "HMM_Island",
  "Regulatory_Feature_Name", "Regulatory_Feature_Group", "DNase_Hypersensitivity_NAME",
  "DNase_Hypersensitivity_Evidence_Count", "OpenChromatin_NAME", "OpenChromatin_Evidence_Count",
  "TFBS_NAME", "TFBS_Evidence_Count")

#Limit manifest to only variables of interest (see code directly above):
Illumina.EPIC.manifest.lim <- Illumina.EPIC.manifest[myvars]

#Now create a smaller version of the manifest that ONLY includes probes with 
#annotated gene.
#First, must change all blanks to NA for omission in later step:
Illumina.EPIC.manifest.lim$UCSC_RefGene_Name <- as.character(Illumina.EPIC.manifest.lim$UCSC_RefGene_Name)
Illumina.EPIC.manifest.lim$UCSC_RefGene_Name[Illumina.EPIC.manifest.lim$UCSC_RefGene_Name==""] <- NA
Illumina.EPIC.manifest.lim$UCSC_RefGene_Name <- as.factor(Illumina.EPIC.manifest.lim$UCSC_RefGene_Name)

#Now subset by complete cases for variable 4 (UCSC Refgene name)
Illumina.EPIC.manifest.lim.1 <- Illumina.EPIC.manifest.lim[complete.cases(Illumina.EPIC.manifest.lim[ , 4]),]
head(Illumina.EPIC.manifest.lim.1)

Clean and annotate data sets–convert the results to a data frame, sort to remove probes with a high standard error, and annotate with gene names, etc

load("./fit_md_pre_fdr.RData")

fit_m.se <- fit_md[(fit_md$Std..Error < 10), ] #Merge data frames of gamlss 
# interaction term model outputs with the EPIC manifest by "probe" variable:
colnames(fit_m.se)[colnames(fit_m.se) == ".id"] <- "probe" #rename column
fit_m.merged <- merge(fit_m.se,Illumina.EPIC.manifest.lim.1,by=c("probe")) 
## Saving 7 x 5 in image

##Plots of most significant CpGs positive and negative, raw BS/oxBS

#Theme for plots: font 28 for figure output
My_Theme = theme(
  axis.title.x = element_text(size = 14),
  axis.text.x = element_text(size = 14),
  axis.text.y = element_text(size = 14),
  axis.title.y = element_text(size = 14),
  plot.title = element_text(hjust = 0.5, size = 16),
  legend.title = element_text(size=14),
  legend.text=element_text(size=14))

#Sort by order of most extreme coefficients and pull most negative and most positive
fit_m.se.f <- fit_m.se[(fit_m.se$fdr_m < 0.05), ]
fit_m.se.term.sort <- fit_m.se.f[order(fit_m.se.f$Estimate),]
extremenegative<-fit_m.se.term.sort[1,]
extremepositive<-fit_m.se.term.sort[length(fit_m.se.term.sort$Estimate),]

extremeneground<-extremenegative %>% mutate(across(where(is.numeric), ~ round(.,5)))
extremeposround<-extremepositive %>% mutate(across(where(is.numeric), ~ round(.,5)))

#Sort by order of FDR and separate into positive and negative
fit_m.se.fdr.sort <- fit_m.se[order(fit_m.se$fdr_m),]
positive<-fit_m.se.fdr.sort[which(fit_m.se.fdr.sort$Estimate>0),]
negative<-fit_m.se.fdr.sort[which(fit_m.se.fdr.sort$Estimate<0),]

#Take out most significant negative and positive, round numbers to five digits. Create lists of probe names
lowestten<-negative[1:10,]
highestten<-positive[1:10,]
lowest<-lowestten %>% mutate(across(where(is.numeric), ~ round(.,5)))
highest<-highestten %>% mutate(across(where(is.numeric), ~round(.,5)))

lowestcgs<-c()
highestcgs<-c()
lowestcgs<-append(lowestcgs, lowest$probe)
highestcgs<-append(highestcgs, highest$probe)

#Load data from EPIC_QC_MLE chunk 14
load("./chunk_14_data.RData")

bs.fp_t<-as.data.frame(t(bs.fp_m_final))
oxbs.fp_t<-as.data.frame(t(oxbs.fp_m_final))

raw_betas<-bind_rows(bs.fp_t, oxbs.fp_t)

#Function to pull betas
pull <- function (data, x, outframe, direction, type){
  z<<-z+1
  outframe<-data.frame(data[,x])
  index<-which(grepl(x, outframe))
  names(outframe)<- paste(direction,z,type, sep="_")
  return(outframe)
  }

#pull from ten lowest
z=0
dir<-"lowest"
type<-"raw"
lowestframeraw<-data.frame(sapply(lowestcgs, function(x) pull (raw_betas, x, lowestframeraw, dir, type)))

type<-"corrected"
z=0
lowestbs<-data.frame(sapply(lowestcgs, function(x) pull(bs.fp_t, x, lowestbs, dir, type)))

#Pull from ten highest
z=0
dir<-"highest"
type<-"raw"
highestframeraw<-data.frame(sapply(highestcgs, function(x) pull (raw_betas, x, highestframeraw, dir, type)))

z=0
type<-"corrected"
highestbs<-data.frame(sapply(highestcgs, function(x) pull(bs.fp_t, x, highestbs, dir, type)))

meta_m_extreme<-meta_m_final
meta_extremebs<-meta_m_final[1:56,]
meta_m_plot<-cbind(meta_m_final, lowestframeraw)
meta_m_plot<-cbind(meta_m_plot, highestframeraw)
meta_bs<-meta_m_final[1:56,]

lowesttitles<-colnames(lowestframeraw)
lowest$title<-lowesttitles
highesttitles<-colnames(highestframeraw)
highest$title<-highesttitles

#Pull BS for extremeneg and extremepos
mostnegative<-bs.fp_t[,extremeneground$probe]
mostpositive<-bs.fp_t[,extremeposround$probe]

meta_extremebs$negative<-mostnegative
meta_extremebs$positive<-mostpositive

#Lowest cg plots
for (yay in 10:1){
  name<-lowest[yay, 8]
lowestcgplotraw<-ggplot(data=meta_m_plot, aes(x=Disease, y = meta_m_plot[,lowesttitles[c(yay)]])) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=meta_m_plot[,lowesttitles[c(yay)]], group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("Raw betas for ",yay, " most significant CpG with negative coefficient", " \n ", lowest[yay,1],", beta coefficient = ", lowest[yay, 2], ", FDR = ", lowest[yay, 6])) +
  xlab(" ") +  ylab("Raw BS/oxBS beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave(filename=str_c("raw_bs_oxbs_lowest_", yay, "_cpg_correct.png"), lowestcgplotraw, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

lowestcgplotraw

#Highest cg plots
for (yay in 10:1){
highestcgplotraw<-ggplot(data=meta_m_plot, aes(x=Disease, y = meta_m_plot[,highesttitles[c(yay)]])) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=meta_m_plot[,highesttitles[c(yay)]], group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("Raw betas for ",yay, " most significant CpG with positive coefficient ", " \n ", highest[yay,1]," beta coefficient = ", highest[yay, 2], ", FDR = ", highest[yay, 6])) +
  xlab(" ") +  ylab("Raw BS/oxBS beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent')) +
  My_Theme

ggsave(filename=str_c("raw_bs_oxbs_highest_", yay, "_cpg_correct.png"), highestcgplotraw, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

highestcgplotraw

#DNAJC6
dnajc6raw<-raw_betas$cg01781382
meta_m_final$DNAJC6raw<- dnajc6raw

DNAJC6plotraw<-ggplot(data=meta_m_final, aes(x=Disease, y = DNAJC6raw)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=DNAJC6raw, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle("Raw BS/oxBS betas for DNAJC6 (cg01781382) \n beta coefficient = -0.10709") +
  xlab(" ") +  ylab("Raw BS/oxBS beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("raw_bx_oxbs_DNAJC6.png", DNAJC6plotraw, width=11.8, height= 8.68, dpi=700, bg="transparent")

DNAJC6plotraw

#PTPRN2
ptprn2raw<-raw_betas$cg06423822
meta_m_final$PTPRN2raw<- ptprn2raw

PTPRN2plotraw<-ggplot(data=meta_m_final, aes(x=Disease, y=PTPRN2raw)) +
geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=PTPRN2raw, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle("Raw BS/oxBS betas for PTPRN2 (cg06423822) \n beta coefficient = 0.16870") +
  xlab(" ") +  ylab("Raw BS/oxBS beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("raw_bx_oxbs_PTPRN2.png", PTPRN2plotraw, width=11.8, height= 8.68, dpi=700, bg="transparent")  

PTPRN2plotraw

##Plots of most significant CpGs both negative and positive mC hmC

#Modify meta file for purposes of plotting
meta_graph<-replace(meta_m_final$Assay_1, meta_m_final$Assay_1=="oxBS", "hmC")
meta_m_final$Assay_1<-meta_graph
meta_graph<-replace(meta_m_final$Assay_1, meta_m_final$Assay_1=="BS", "mC")
meta_m_final$Assay_1<-meta_graph

meta_graph<-replace(meta_m_extreme$Assay_1, meta_m_extreme$Assay_1=="oxBS", "hmC")
meta_m_extreme$Assay_1<-meta_graph
meta_graph<-replace(meta_m_extreme$Assay_1, meta_m_extreme$Assay_1=="BS", "mC")
meta_m_extreme$Assay_1<-meta_graph

#prep bvals for pulling
bvals.t <- as.data.frame(t(bvals.complete))

negative_adjusted<-bvals.t[,extremeneground$probe]
positive_adjusted<-bvals.t[,extremeposround$probe]
meta_m_extreme$negative<-negative_adjusted
meta_m_extreme$positive<-positive_adjusted

#Pull adjusted bvals for lowest
z=0
dir<-"lowest"
type<-"corrected"
lowestframeadjusted<-data.frame(sapply(lowestcgs, function(x) pull (bvals.t, x, lowestframeadjusted, dir, type)))

#Pull adjusted bvals for highest
z=0
dir<-"highest"
highestframeadjusted<-data.frame(sapply(highestcgs, function(x) pull (bvals.t, x, highestframeadjusted, dir, type)))

#Bind data to meta
meta_m_adj<-cbind(meta_m_final, lowestframeadjusted)
meta_m_adj<-cbind(meta_m_adj, highestframeadjusted)

lowesttitles<-colnames(lowestframeadjusted)
lowest$title<-lowesttitles
highesttitles<-colnames(highestframeadjusted)
highest$title<-highesttitles


#Plot lowest
for (yay in 10:1){
  name<-lowest[yay, 8]
lowestcgplot<-ggplot(data=meta_m_adj, aes(x=Disease, y = meta_m_adj[,lowesttitles[c(yay)]])) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=meta_m_adj[,lowesttitles[c(yay)]], group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("mC/hmC betas for ",yay, " most significant CpG with negative coefficient ", " \n ", lowest[yay,1]," beta coefficient = ", lowest[yay, 2], ", FDR = ", lowest[yay, 6])) +
  xlab(" ") +  ylab("MLE-adjusted values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave(filename=str_c("MLE-adjusted_lowest_", yay, "_cpg_correct.png"), lowestcgplot, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

lowestcgplot

#Plot highest
for (yay in 10:1){
highestcgplot<-ggplot(data=meta_m_adj, aes(x=Disease, y = meta_m_adj[,highesttitles[c(yay)]])) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=meta_m_adj[,highesttitles[c(yay)]], group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("mC/hmC betas for ",yay, " most significant CpG with positive coefficient ", " \n ", highest[yay,1]," beta coefficient = ", highest[yay, 2], ", FDR = ", highest[yay, 6])) +
  xlab(" ") +  ylab("MLE-adjusted beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent')) +
  My_Theme

ggsave(filename=str_c("MLE-adjusted_highest_", yay, "_cpg_correct.png"), highestcgplot, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

highestcgplot

#DNAJC6
dnajc6adjusted<-bvals.t$cg01781382

meta_m_final$DNAJC6adjusted <- dnajc6adjusted

DNAJC6plot<-ggplot(data=meta_m_final, aes(x=Disease, y = DNAJC6adjusted)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=DNAJC6adjusted, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle("mC/hmC betas for DNAJC6 (cg01781382) \n beta coefficient = -0.10709") +
  xlab(" ") +  ylab("MLE-adjusted beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("MLE-adjusted_DNAJC6.png", DNAJC6plot, width=11.8, height= 8.68, dpi=700, bg="transparent")

DNAJC6plot

#PTPRN2
ptprn2adjusted<-bvals.t$cg06423822

meta_m_final$PTPRN2adjusted <- ptprn2adjusted

PTPRN2plot<-ggplot(data=meta_m_final, aes(x=Disease, y = PTPRN2adjusted)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=PTPRN2adjusted, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle("mC/hmC betas for PTPRN2 (cg06423822) \n beta coefficient = 0.16870") +
  xlab(" ") +  ylab("MLE-adjusted beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("MLE-adjusted_PTPRN2.png", PTPRN2plot, width=11.8, height= 8.68, dpi=700, bg="transparent")

PTPRN2plot

##Plot most significant CpGs both positive and negative, mC hmC BS

#Make three-level plots
meta_bs_plot<-cbind(meta_bs, lowestbs)
meta_bs_plot<-cbind(meta_bs_plot, highestbs)

dnajc6bs<-bs.fp_t$cg01781382
meta_bs_plot$DNAJC6adjusted <- dnajc6bs
meta_m_adj$DNAJC6raw<-NULL
meta_m_adj$DNAJC6adjusted<-dnajc6adjusted

ptprn2bs<-bs.fp_t$cg06423822
meta_bs_plot$PTPRN2adjusted <- ptprn2bs
meta_m_adj$PTPRN2raw<-NULL
meta_m_adj$PTPRN2adjusted<-ptprn2adjusted

meta_three<-rbind(meta_m_adj, meta_bs_plot)
meta_for<-rbind(meta_m_adj, meta_bs_plot)

meta_plot_extreme<-rbind(meta_m_extreme, meta_extremebs)


for (yay in 10:1){
  name<-lowest[yay, 8]
lowestcgplot<-ggplot(data=meta_three, aes(x=Disease, y = meta_three[,lowesttitles[c(yay)]])) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=meta_three[,lowesttitles[c(yay)]], group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("Betas for ",yay, " most significant CpG with negative coefficient ", " \n ", lowest[yay,1]," beta coefficient = ", lowest[yay, 2], ", FDR = ", lowest[yay, 6])) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave(filename=str_c("three-layer-lowest_", yay, "_cpg.png"), lowestcgplot, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

lowestcgplot

for (yay in 10:1){
highestcgplot<-ggplot(data=meta_three, aes(x=Disease, y = meta_three[,highesttitles[c(yay)]])) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=meta_three[,highesttitles[c(yay)]], group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("Betas for ",yay, " most significant CpG with positive coefficient ", " \n ", highest[yay,1]," beta coefficient = ", highest[yay, 2], ", FDR = ", highest[yay, 6])) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent')) +
  My_Theme

ggsave(filename=str_c("three-layer_highest_", yay, "_cpg.png"), highestcgplot, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

highestcgplot

DNAJC6plot<-ggplot(data=meta_for, aes(x=Disease, y = DNAJC6adjusted)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=DNAJC6adjusted, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle("Betas for DNAJC6 (cg01781382) \n beta coefficient = -0.10709") +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("three-level_DNAJC6.png", DNAJC6plot, width=11.8, height= 8.68, dpi=700, bg="transparent")

DNAJC6plot

PTPRN2plot<-ggplot(data=meta_for, aes(x=Disease, y = PTPRN2adjusted)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=PTPRN2adjusted, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle("Betas for PTPRN2 (cg06423822) \n beta coefficient = 0.16870") +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme
 
ggsave("three-level_PTPRN2.png", PTPRN2plot, width=11.8, height= 8.68, dpi=700, bg="transparent")

PTPRN2plot

extremenegativeplot<-ggplot(data=meta_plot_extreme, aes(x=Disease, y = negative)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=negative, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("Betas for most negative beta coefficient ", extremeneground$probe, " \n beta coefficient = ", extremeneground$Estimate, ", FDR = ", extremeneground$fdr_m)) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("three-level_mostnegative.png", extremenegativeplot, width=11.8, height= 8.68, dpi=700, bg="transparent")

extremenegativeplot

extremepositiveplot<-ggplot(data=meta_plot_extreme, aes(x=Disease, y = positive)) + 
  geom_point(aes(color=Assay_1)) +
  geom_smooth(mapping = aes(x=Disease, y=positive, group=Assay_1, color=Assay_1, fill= Assay_1), method = lm) +
  ggtitle(paste0("Betas for most positive beta coefficient ", extremeposround$probe, " \n beta coefficient = ", extremeposround$Estimate, ", FDR = ", extremeposround$fdr_m)) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Assay", fill= "Assay") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("three-level_mostpositive.png", extremepositiveplot, width=11.8, height= 8.68, dpi=700, bg="transparent")

extremepositiveplot

##Plot most significant CpGs with modification on the X

for (yay in 10:1){
  name<-lowest[yay, 8]
lowestcgplot<-ggplot(data=meta_m_adj, aes(x=Assay_1, y = meta_m_adj[,lowesttitles[c(yay)]])) + 
  geom_point(aes(color=Disease)) +
  geom_smooth(mapping = aes(x=Assay_1, y=meta_m_adj[,lowesttitles[c(yay)]], group=Disease, color=Disease, fill= Disease), method = lm) +
  ggtitle(paste0("Betas for ",yay, " most significant CpG with negative coefficient ", " \n ", lowest[yay,1]," beta coefficient = ", lowest[yay, 2], ", FDR = ", lowest[yay, 6])) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Group", fill= "Group") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave(filename=str_c("modonx_lowest_", yay, "_cpg.png"), lowestcgplot, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

lowestcgplot

for (yay in 10:1){
  name<-highest[yay, 8]
highestcgplot<-ggplot(data=meta_m_adj, aes(x=Assay_1, y = meta_m_adj[,highesttitles[c(yay)]])) + 
  geom_point(aes(color=Disease)) +
  geom_smooth(mapping = aes(x=Assay_1, y=meta_m_adj[,highesttitles[c(yay)]], group=Disease, color=Disease, fill= Disease), method = lm) +
  ggtitle(paste0("Betas for ",yay, " most significant CpG with positive coefficient ", " \n ", highest[yay,1]," beta coefficient = ", highest[yay, 2], ", FDR = ", highest[yay, 6])) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Group", fill= "Group") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave(filename=str_c("modonx_highest_", yay, "_cpg.png"), highestcgplot, width=11.8, height = 8.68, dpi=700, bg="transparent")
}

highestcgplot

DNAJC6plot<-ggplot(data=meta_m_adj, aes(x=Assay_1, y = DNAJC6adjusted)) + 
  geom_point(aes(color=Disease)) +
  geom_smooth(mapping = aes(x=Assay_1, y=DNAJC6adjusted, group=Disease, color=Disease, fill= Disease), method = lm) +
  ggtitle("Beta values for DNAJC6 (cg01781382) \n beta coefficient = -0.10709") +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Group", fill= "Group") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("modonx_DNAJC6.png", DNAJC6plot, width=11.8, height= 8.68, dpi=700, bg="transparent")

DNAJC6plot

PTPRN2plot<-ggplot(data=meta_m_adj, aes(x=Assay_1, y = PTPRN2adjusted)) + 
  geom_point(aes(color=Disease)) +
  geom_smooth(mapping = aes(x=Assay_1, y=PTPRN2adjusted, group=Disease, color=Disease, fill= Disease), method = lm) +
  ggtitle("Beta values for PTPRN2 (cg06423822) \n beta coefficient = 0.16870") +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Group", fill= "Group") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("modonx_PTPRN2.png", PTPRN2plot, width=11.8, height= 8.68, dpi=700, bg="transparent")

PTPRN2plot

extremenegativeplot<-ggplot(data=meta_m_extreme, aes(x=Assay_1, y = negative)) + 
  geom_point(aes(color=Disease)) +
  geom_smooth(mapping = aes(x=Assay_1, y=negative, group=Disease, color=Disease, fill= Disease), method = lm) +
  ggtitle(paste0("Betas for most negative beta coefficient ", extremeneground$probe, " \n beta coefficient = ", extremeneground$Estimate, ", FDR = ", extremeneground$fdr_m)) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Group", fill= "Group") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("mostnegative_modonx.png", extremenegativeplot, width=11.8, height= 8.68, dpi=700, bg="transparent")

extremenegativeplot

extremepositiveplot<-ggplot(data=meta_m_extreme, aes(x=Assay_1, y = positive)) + 
  geom_point(aes(color=Disease)) +
  geom_smooth(mapping = aes(x=Assay_1, y=positive, group=Disease, color=Disease, fill= Disease), method = lm) +
  ggtitle(paste0("Betas for most positive beta coefficient ", extremeposround$probe, " \n beta coefficient = ", extremeposround$Estimate, ", FDR = ", extremeposround$fdr_m)) +
  xlab(" ") +  ylab("Beta values") +
  labs(color="Group", fill= "Group") +
  ylim(0,1) +
  theme(plot.background = element_rect(fill='transparent'))+
  My_Theme

ggsave("mostpositive_modonx.png", extremepositiveplot, width=11.8, height= 8.68, dpi=700, bg="transparent")

extremepositiveplot

##Histogram of chromatin annotations

My_Theme2 = theme(
  axis.title.x = element_text(size = 14),
  axis.text.x = element_text(size = 14),
  axis.text.y = element_text(size = 14),
  axis.title.y = element_text(size = 14),
  plot.title = element_text(hjust = 0.5, size = 16),
  legend.title = element_text(size=14),
  legend.text=element_text(size=14),
            panel.background = element_rect(fill='transparent'),
            plot.background = element_rect(fill='transparent', color=NA),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.background = element_rect(fill='transparent'),
            legend.box.background = element_rect(fill='transparent'))

fit_full_merged<-read.csv("./240515_Interaction_Male_DMCs_FDR_with_PMI_Manifest_ChromHMM_Annotated.csv")

annot.dmcs <- fit_full_merged %>% dplyr::select(annot.symbol)

annot.dmcs<-annot.dmcs %>% drop_na(annot.symbol)


clean_df <- map_df(annot.dmcs, ~ gsub(".*TxWk.*", "weak transcription", .x))
clean_df <- map_df(clean_df, ~ gsub(".*GapArtf.*", "assembly gaps & artifacts", .x))
clean_df <- map_df(clean_df, ~ gsub(".*BivProm.*", "bivalent promoters", .x))
clean_df <- map_df(clean_df, ~ gsub(".*Quies.*", "quiescent", .x))
clean_df <- map_df(clean_df, ~ gsub(".*Acet.*", "acetylation marks", .x))
clean_df <- map_df(clean_df, ~ gsub(".*ReprPC.*", "polycomb repressed & \n open chromatin", .x))
clean_df <- map_df(clean_df, ~ gsub(".*HET.*", "heterochromatin", .x))
clean_df <- map_df(clean_df, ~ gsub(".*EnhA.*", "active enhancers", .x))
clean_df <- map_df(clean_df, ~ gsub(".*Tx.*", "transcription", .x))
clean_df <- map_df(clean_df, ~ gsub(".*TxEnh.*", "transcribed enhancers", .x))
clean_df <- map_df(clean_df, ~ gsub(".*EnhWk.*", "weak enhancer", .x))
clean_df <- map_df(clean_df, ~ gsub(".*TxEx.*", "transcription and exons", .x))
clean_df <- map_df(clean_df, ~ gsub(".*znf.*", "zinc finger genes", .x))
clean_df <- map_df(clean_df, ~ gsub(".*PromF.*", "promoter flank", .x))
clean_df <- map_df(clean_df, ~ gsub(".*DNase.*", "DNase", .x))
clean_df <- map_df(clean_df, ~ gsub(".*TSS.*", "transcription start sites", .x))

annotation_plot <- ggplot() + 
  aes(clean_df$annot.symbol)+
      geom_bar(fill = '#003546') + 
      labs(title="Chromatin Annotation Types", x ="Annotation Type", y = "Frequency") +
      theme(axis.title = element_text(size=30), axis.text=element_text(size=24), plot.title = element_text(size = 26))+ 
  My_Theme2+
      coord_flip() 

ggsave(filename="featureannotations.png", annotation_plot,width=15, height= 10, dpi=700, bg="transparent")

annotation_plot

##Histogram of intragenic annotations

annot.dmcs.intra <- fit_m.merged %>% dplyr::select(UCSC_RefGene_Group)

annot.dmcs.intra<-annot.dmcs.intra %>% drop_na(UCSC_RefGene_Group)

featurelist<- annot.dmcs.intra %>% mutate(rn=row_number()) %>%
  separate_rows(UCSC_RefGene_Group, sep=";") %>%
  distinct %>%
  group_by(rn) %>%
  summarize(col=str_c(UCSC_RefGene_Group, collapse = ";"), .groups = "drop")

featurelist<-select(featurelist, col)

#separate genes into single-gene rows. Note that the separator is ";"
annot.dmcs.unique <- separate_rows(featurelist, col, sep=";")

#remove NA, if applicable
annot.dmcs.unique.n <- drop_na(annot.dmcs.unique, col)

#remove whitespace, if applicable
annot.dmcs.unique.clean <- as.data.frame(apply(annot.dmcs.unique.n, 1, str_remove_all, " "))

#names the columns for later merge
colnames(annot.dmcs.unique.clean) = c('c1')


annot_features1<-annot.dmcs.unique.clean %>% drop_na(c1)

annot_features <- annot_features1$c1

annotationcount<-ggplot() +
  aes(annot_features) +
  geom_bar(fill = '#003546') +
  xlab("Annotation type") +  ylab("Frequency") +
  theme(plot.background = element_rect(fill='transparent'), panel.background= element_rect(fill='transparent'))+ My_Theme2 + coord_flip()

ggsave(filename="featureannotations_intragenic.png", annotationcount,width=11.8, height= 8.6, dpi=700, bg="transparent")

annotationcount