methylation_analysis.R

## ---------------------------
##
## Script name: DNA Methylation analysis for T21-XIST iPSCs
##
## Purpose of script:
##
## Author: Prakhar Bansal
##
## Date Created: 2021-08-03
##
## Copyright (c) Prakhar Bansal, 2021
## Email: pbansal@uchc.edu
##
## ---------------------------
##
## Notes:
##   
##
## ---------------------------

## set working directory
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

## ---------------------------

## load up the packages we will need:  (uncomment as required)

require(tidyverse); theme_set(theme_classic() + theme(plot.title=element_text(hjust=0.5)))
require(ggsci)
devtools::install_github("achilleasNP/IlluminaHumanMethylationEPICmanifest") 
devtools::install_github("achilleasNP/IlluminaHumanMethylationEPICanno.ilm10b5.hg38")
## ---------------------------

library(minfi)
library(IlluminaHumanMethylationEPICmanifest)

sample_info_df <- tibble(sample_name=c("c7 6wk -dox", "c7 3+/3- dox", "c7 6wk +dox"),
                         clone="c7",
                         dox_treatment=c("6wk-", "3+/3-", "6wk+"),
                         array_id="205103220024",
                         position=6:8,
                         location="input_data/dna_methylation_arrays/205103220024") %>%
  mutate(Sentrix_Position=sprintf("R0%dC01", position),
         Basename=sprintf("%s/%s_%s", location, array_id, Sentrix_Position))
write_csv(sample_info_df, "dname_sample_info.csv")

RGset <- read.metharray.exp(targets = sample_info_df)
RGset@annotation = c(array = "IlluminaHumanMethylationEPIC", annotation = "ilm10b5.hg38")
gmset <- preprocessIllumina(RGset)
gmset <- mapToGenome(gmset) 

beta_mat <- getBeta(gmset) %>%
  as.data.frame() %>%
  rownames_to_column(var="rn")

annotation_df <- getAnnotation(gmset) %>% as.data.frame()

beta_df <- beta_mat |> 
  left_join(annotation_df |> 
              dplyr::select(Name, chr=CHR_hg38, pos=Start_hg38, strand, gene_name=UCSC_RefGene_Name, refgene_group=UCSC_RefGene_Group)|> 
              mutate(pos=as.numeric(pos)), 
            by=c("rn"="Name") ) |> 
  relocate(chr, pos, strand, gene_name, refgene_group, .after=rn) %>%
  dplyr::rename(nodox=`205103220024_R06C01`,
                withdrawal=`205103220024_R07C01`,
                dox=`205103220024_R08C01`) %>%
  mutate(`dox-nodox`=dox-nodox,
         `withdrawal-nodox`=withdrawal-nodox,
         `withdrawal-dox`=withdrawal-dox)



beta_df |> 
  filter(chr=="chr21" & str_detect(refgene_group, "TSS1500|TSS200|5_UTR")) |> 
  pivot_longer(cols = c(`nodox`, `dox`,`withdrawal`), names_to="comparison", values_to="value") |> 
  mutate(comparison = case_when(comparison == "nodox" ~ "-dox",
                                comparison == "dox" ~ "+dox",
                                comparison == "withdrawal" ~ "w/d")) |> 
  filter(value > 0) |> 
  ggplot(aes(x = comparison, y = value, fill=comparison, color=comparison)) + 
  see::geom_violinhalf(position = position_nudge(x=.2)) +
  geom_point(
    size = .005,
    alpha = .15,
    position = position_jitter(
      seed = 1, width = .05
    )
  ) + 
  geom_boxplot(
    width = .1, 
    outlier.shape = NA,
    alpha=.2,
    color="black",
    size=.3
  ) +
  annotate("text", x=1.4, size=2.0, y=0.0655, label = "0.07") +
  annotate("text", x=1.4, size=2.0, y=0.799, label = "0.80") +
  annotate("text", x=2.4, size=2.0, y=0.412, label = "0.41") +
  annotate("text", x=2.4, size=2.0, y=0.817, label = "0.82") +
  annotate("text", x=3.4, size=2.0, y=0.391, label = "0.39") +
  annotate("text", x=3.4, size=2.0, y=0.790, label = "0.79") +
  coord_cartesian(xlim = c(1.2, NA),ylim=c(0, 1), clip = "on") +
  scale_fill_manual(values=c("#BF5AC9", "#0D7174", "#FE7F01"), limits=c("-dox", "+dox", "w/d"), guide="none") +
  scale_color_manual(values=c("#BF5AC9", "#0D7174", "#FE7F01"), limits=c("-dox", "+dox", "w/d"), guide="none") +
  scale_x_discrete(limits=c("-dox", "+dox", "w/d")) +
  labs(title="DNAme of chr21 Promoters",
       y="DNAme") +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle=20, vjust=1, hjust=1),
        text = element_text(size=7),
        axis.title = element_text(size=7),
        plot.title = element_text(size=7, face="bold", margin=margin(t=2, 0)),
        axis.text.y = element_text(size=5),
        line=element_line(size=.3),
        panel.spacing = unit(.2, "mm"),
        strip.text = element_text(size=6, margin=margin(t=0,r=1,b=0,l=1)),
        strip.background = element_rect(size = .3),
        plot.margin = margin(0.2))
ggsave("figs/dname_hsa21_promoter.pdf", width=3.5, height=1.5, units="in")


beta_df |> 
  filter(chr=="chr21" & str_detect(refgene_group, "TSS1500|TSS200|5_UTR")) |> 
  pivot_longer(cols = c(`dox-nodox`, `withdrawal-nodox`,`withdrawal-dox`), names_to="comparison", values_to="value") |> 
  mutate(comparison = factor(comparison, levels=c("dox-nodox", "withdrawal-dox", "withdrawal-nodox"))) |> 
  ggplot(aes(x=pos, y=value, color=comparison)) +
  geom_point(size=1, alpha=.3) +
  facet_grid(rows=vars(comparison), labeller = labeller(comparison=c("dox-nodox"="dox-nd", "withdrawal-nodox"="wd-nd", "withdrawal-dox"="wd-dox"))) +
  scale_color_manual(values=c("#000000", "#0D7174", "#FE7F01"), limits=c("withdrawal-dox", "dox-nodox", "withdrawal-nodox"), guide="none") +
  geom_hline(yintercept = 0, linetype="solid", size=.3) +
  coord_cartesian(ylim=c(-.5, .5)) +
  labs(x="Chr21 Position",
       y="Change in DNAme",
       title="Change in DNAme at 6 weeks")
ggsave("figs/dname_chr21_promoter_change.png", width=4, height=3, units="in")

#WRB gene methylation
wrb_methylation_df <- beta_df |> 
  filter(str_detect(gene_name, "WRB")) |> 
  mutate(pos_str=as.character(pos)) |> 
  pivot_longer(cols = c(`nodox`, `dox`,`withdrawal`), names_to="condition", values_to="value") |> 
  mutate(pos_str=fct_reorder(pos_str, pos))
cpgs <- tibble(island_name = c("CPGI", "CPGII", "CPGIII"),
               start = c(39380162, 39385677, 39388701),
               stop = c(39380899, 39385974, 39388903),
               start_cpg = c("39380170", "39385764", "39388704"),
               stop_cpg = c("39380828", "39385972", "39388869"))
ggplot(wrb_methylation_df) +
  geom_point(aes(x=pos_str, y=value, color=condition)) +
  geom_rect(data = cpgs,
            mapping = aes(xmin=factor(start_cpg, levels = levels(wrb_methylation_df$pos_str)) |> as.numeric() - .5, 
                          xmax=factor(stop_cpg, levels = levels(wrb_methylation_df$pos_str)) |> as.numeric() + .5,
                          ymin=0, ymax=1, fill = island_name), alpha = .4) +
  geom_point(aes(x=pos_str, y=value, color=condition)) +
  scale_color_manual(values=c("#BF5AC9", "#0D7174", "#FE7F01"), limits=c("nodox", "dox", "withdrawal")) +
  labs(title="WRB methylation") +
  geom_hline(yintercept = c(1/3, 2/3), linetype = "dashed") +
  theme(axis.text.x = element_text(angle=45, hjust=1))

beta_df |> 
  filter(rn %in% c("cg26710963", "cg09916765")) |> 
  pivot_longer(cols = c(`nodox`, `dox`,`withdrawal`), names_to="condition", values_to="value") |> 
  ggplot(aes(x=fct_reorder(rn, pos), y=value, color=condition)) +
  geom_point() +
  scale_color_manual(values=c("#BF5AC9", "#0D7174", "#FE7F01"), limits=c("nodox", "dox", "withdrawal")) +
  coord_cartesian(ylim=c(0,1)) +
  labs(title="WRB methylation")


# plot changes in promoter methylation for each chromosome
beta_df |> 
  filter(!is.na(chr) & 
           str_detect(refgene_group, "TSS1500|TSS200|5_UTR")) |>
  group_by(chr) |> 
  summarise(
    gain=sum(dox-nodox > .1, na.rm=T),
            nochange=sum(between(dox-nodox, -0.1, .1), na.rm=T),
            loss=sum(dox-nodox < -.1, na.rm=T)) |>
  pivot_longer(cols = -chr, names_to = "direction", values_to="count") |> 
  mutate(direction=factor(direction, levels=c("gain", "nochange", "loss")),
         chr=factor(str_replace(chr, "chr", ""), levels = c("1", "2", "3", "4", "5", 
                                    "6", "7", "8", "9", "10", 
                                    "11", "12", "13", "14", "15", 
                                    "16", "17", "18", "19", "20", 
                                    "21", "22", "X", "Y"))) |> 
  filter(chr %in% c("1", "2", "3", "4", "5", 
                    "6", "7", "8", "9", "10", 
                    "11", "12", "13", "14", "15", 
                    "16", "17", "18", "19", "20", 
                    "21", "22", "X", "Y")) |> 
  ggplot(aes(x=chr, y=count, fill=direction)) +
  geom_col(position=position_fill()) +
  labs(title="Global Changes in DNAme of Promoters",
       x="Chromosome",
       y="% Probes",
       fill="Direction") +
  geom_hline(yintercept = .5, linetype="dashed", size=.3) +
  scale_fill_manual(values=pal_jco()(4)[c(1,3,4)],
                    limits=c("gain", "nochange", "loss"),
                    labels=c("Gain", "~0", "Loss")) +
  theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1, size=8),
        plot.margin = margin(t=.5,0,0,0, "mm"),
        legend.margin = margin(l=-3, unit = "mm"),
        axis.title.x = element_blank())
+
  theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1),
        text = element_text(size=7),
        axis.title = element_text(size=7),
        plot.title = element_text(size=7, face="bold", margin=margin(t=2, 0)),
        axis.text.y = element_text(size=5),
        line=element_line(size=.3),
        panel.spacing = unit(.2, "mm"),
        strip.text = element_text(size=6, margin=margin(t=0,r=1,b=0,l=1)),
        strip.background = element_rect(size = .3),
        plot.margin = margin(0,0,0,0,"mm"),
        legend.key.size = unit(2, "mm"),
        legend.box.margin = margin(l=-10),
        legend.margin = margin(0))
ggsave("figs/global_changes_in_promoter_dname.pdf", width=4, height=1.3, units="in")