NanoImprint_2_newbedfiles.Rmd

---
title: "NanoImprint"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output: html_document
---

reference genome T2T-CHM13v2.0

```{r include=FALSE}
library(readxl)
library(tidyverse, warning(F))
require(gridExtra)
library(ggplot2, warning(F))
library(edgeR, warning(F))
library(knitr)
library(kableExtra)
```


```{r echo=TRUE}
#import data
sample_df <- read.table("Your_file.methyl.filtered.bed", quote="\"", comment.char="")
hap1 <- read.table("Your_file_hap_1.methyl.filtered.bed", quote="\"", comment.char="")
hap2 <- read.table("Your_file_hap_2.methyl.filtered.bed", quote="\"", comment.char="")

#filter to only 5mC
sample_df<- sample_df[sample_df$V4== "m",]
hp1<- hap1[hap1$V4== "m",]
hp2<- hap2[hap2$V4== "m",]
```


```{r include=FALSE}
#import ctrls
ctrls= read_excel("ctrls_2.xlsx",1)
ctrls<- ctrls [, c(2:4)]


#dataset with selected columns + filtering
sorted_df <- sample_df[,c(2,5,11)]
colnames(sorted_df) <- c("start_position", "coverage", "methylation_frequency")


#dataset (hp1) with selected cloumns + filtering
hp1 <- hp1[,c(2,5,11)]
colnames(hp1) <- c("start_position", "coverage", "methylation_frequency")

#dataset (hp1) with selected cloumns + filtering
hp2 <- hp2[,c(2,5,11)]
colnames(hp2)<- c("start_position", "coverage", "methylation_frequency")




#Imprinted regions with T2T coordinates

PLAGL1_alt_TSS_DMR <- sorted_df %>% filter(between(start_position,145200200,145201432)) # MS-MLPA probe
PLAGL1_ctrl <- ctrls%>%filter(between(start_position,145200200,145201432))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
PLAGL1_ctrl_filtered <- PLAGL1_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
PLAGL1_alt_TSS_DMR_filtered <- PLAGL1_alt_TSS_DMR %>% semi_join(PLAGL1_ctrl_filtered, by = "start_position")

hp1_PLAGL1_alt_TSS_DMR <- hp1 %>% filter(between(start_position,145200200,145201432)) 
hp1_PLAGL1_alt_TSS_DMR_filtered<- hp1_PLAGL1_alt_TSS_DMR %>% semi_join(PLAGL1_ctrl_filtered, by = "start_position")

hp2_PLAGL1_alt_TSS_DMR <- hp2 %>% filter(between(start_position,145200200,145201432)) 
hp2_PLAGL1_alt_TSS_DMR_filtered<- hp2_PLAGL1_alt_TSS_DMR %>% semi_join(PLAGL1_ctrl_filtered, by = "start_position")



GRB10_alt_TSS_DMR<- sorted_df %>% filter(between(start_position,50943355,50944182)) # alders og vævs specifik
GRB10_ctrl <- ctrls %>% filter(between(start_position,50943355,50944182))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
GRB10_ctrl_filtered <- GRB10_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
GRB10_alt_TSS_DMR_filtered <- GRB10_alt_TSS_DMR %>% semi_join(GRB10_ctrl_filtered, by = "start_position")

hp1_GRB10_alt_TSS_DMR <- hp1 %>% filter(between(start_position,50943355,50944182)) 
hp1_GRB10_alt_TSS_DMR_filtered<- hp1_GRB10_alt_TSS_DMR %>% semi_join(GRB10_ctrl_filtered, by = "start_position")

hp2_GRB10_alt_TSS_DMR <- hp2 %>% filter(between(start_position,50943355,50944182)) 
hp2_GRB10_alt_TSS_DMR_filtered<- hp2_GRB10_alt_TSS_DMR %>% semi_join(GRB10_ctrl_filtered, by = "start_position")


MEST_alt_TSS_DMR<- sorted_df %>% filter(between(start_position,131804587,131807000))# MS-MLPA probe
MEST_ctrl <- ctrls %>% filter(between(start_position,131804587,131807000))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
MEST_ctrl_filtered <- MEST_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
MEST_alt_TSS_DMR_filtered <- MEST_alt_TSS_DMR %>% semi_join(MEST_ctrl_filtered, by = "start_position")

hp1_MEST_alt_TSS_DMR<- hp1 %>% filter(between(start_position,131804587,131807000))
hp1_MEST_alt_TSS_DMR_filtered <- hp1_MEST_alt_TSS_DMR %>% semi_join(MEST_ctrl_filtered, by = "start_position")

hp2_MEST_alt_TSS_DMR<- hp2 %>% filter(between(start_position,131804587,131807000))
hp2_MEST_alt_TSS_DMR_filtered <- hp2_MEST_alt_TSS_DMR %>% semi_join(MEST_ctrl_filtered, by = "start_position")



H19_IGF2_IG_DMR<- sorted_df %>% filter(between(start_position,2085438,2090906))# MS-MLPA probe
H19_ctrl <- ctrls%>% filter(between(start_position,2085438,2090906))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
H19_ctrl_filtered <- H19_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
H19_IGF2_IG_DMR_filtered <- H19_IGF2_IG_DMR %>% semi_join(H19_ctrl_filtered, by = "start_position")

hp1_H19_IGF2_IG_DMR<- hp1 %>% filter(between(start_position,2085438,2090906))
hp1_H19_IGF2_IG_DMR_filtered <- hp1_H19_IGF2_IG_DMR %>% semi_join(H19_ctrl_filtered, by = "start_position")

hp2_H19_IGF2_IG_DMR<- hp2 %>% filter(between(start_position,2085438,2090906))
hp2_H19_IGF2_IG_DMR_filtered <- hp2_H19_IGF2_IG_DMR %>% semi_join(H19_ctrl_filtered, by = "start_position")

KCNQ1OT1_TSS_DMR<- sorted_df %>% filter(between(start_position,2788106,2790240))# MS-MLPA probe
KCNQ_ctrl <- ctrls%>% filter(between(start_position,2788106,2790240))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
KCNQ_ctrl_filtered <- KCNQ_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
KCNQ1OT1_TSS_DMR_filtered <- KCNQ1OT1_TSS_DMR %>% semi_join(KCNQ_ctrl_filtered, by = "start_position")

hp1_KCNQ1OT1_TSS_DMR<- hp1 %>% filter(between(start_position,2788106,2790240))
hp1_KCNQ1OT1_TSS_DMR_filtered <- hp1_KCNQ1OT1_TSS_DMR %>% semi_join(KCNQ_ctrl_filtered, by = "start_position")

hp2_KCNQ1OT1_TSS_DMR<- hp2 %>% filter(between(start_position,2788106,2790240))
hp2_KCNQ1OT1_TSS_DMR_filtered <- hp2_KCNQ1OT1_TSS_DMR %>% semi_join(KCNQ_ctrl_filtered, by = "start_position")


#MEG3_DLK_IG_DMR<- sorted_df %>% filter(between(start_position,95044223,95046285))# 
#MEG_DLK_ctrl <- ctrls%>% filter(between(start_position,95044223,95046285))

MEG3_TSS_DMR<- sorted_df %>% filter(between(start_position,95059075,95062529))# MS-MLPA probe
MEG3_ctrl <- ctrls%>% filter(between(start_position,95059075,95062529))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
MEG3_ctrl_filtered <- MEG3_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
MEG3_TSS_DMR_filtered <- MEG3_TSS_DMR %>% semi_join(MEG3_ctrl_filtered, by = "start_position")

hp1_MEG3_TSS_DMR<- hp1 %>% filter(between(start_position,95059075,95062529))
hp1_MEG3_TSS_DMR_filtered <- hp1_MEG3_TSS_DMR %>% semi_join(MEG3_ctrl_filtered, by = "start_position")

hp2_MEG3_TSS_DMR<- hp2 %>% filter(between(start_position,95059075,95062529))
hp2_MEG3_TSS_DMR_filtered <- hp2_MEG3_TSS_DMR %>% semi_join(MEG3_ctrl_filtered, by = "start_position")


MEG8_Int2_DMR <- sorted_df %>% filter(between(start_position,95140000,95140384))# MS-MLPA probe
MEG8_ctrl <- ctrls%>% filter(between(start_position,95140000,95140384))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
MEG8_ctrl_filtered <- MEG8_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
MEG8_Int2_DMR_filtered <- MEG8_Int2_DMR %>% semi_join(MEG8_ctrl_filtered, by = "start_position")

hp1_MEG8_Int2_DMR <- hp1 %>% filter(between(start_position,95140000,95140384))
hp1_MEG8_Int2_DMR_filtered <- hp1_MEG8_Int2_DMR %>% semi_join(MEG8_ctrl_filtered, by = "start_position")

hp2_MEG8_Int2_DMR <- hp2 %>% filter(between(start_position,95140000,95140384))
hp2_MEG8_Int2_DMR_filtered <- hp2_MEG8_Int2_DMR %>% semi_join(MEG8_ctrl_filtered, by = "start_position")


MAGEL2_TSS_DMR<- sorted_df %>% filter(between(start_position,21381784,21383056))# MS-MLPA probe
MAGEL2_ctrl <- ctrls%>% filter(between(start_position,21381784,21383056))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
MAGEL2_ctrl_filtered <- MAGEL2_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
MAGEL2_TSS_DMR_filtered <- MAGEL2_TSS_DMR %>% semi_join(MAGEL2_ctrl_filtered, by = "start_position")

hp1_MAGEL2_TSS_DMR<- hp1 %>% filter(between(start_position,21381784,21383056))
hp1_MAGEL2_TSS_DMR_filtered <- hp1_MAGEL2_TSS_DMR %>% semi_join(MAGEL2_ctrl_filtered, by = "start_position")

hp2_MAGEL2_TSS_DMR<- hp2 %>% filter(between(start_position,21381784,21383056))
hp2_MAGEL2_TSS_DMR_filtered <- hp2_MAGEL2_TSS_DMR %>% semi_join(MAGEL2_ctrl_filtered, by = "start_position")

SNURF_TSS_DMR<- sorted_df %>% filter(between(start_position,22691523,22693493))#MS-MLPA probe
SNURF_ctrl <- ctrls  %>% filter(between(start_position,22691523,22693493))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
SNURF_ctrl_filtered <- SNURF_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
SNURF_TSS_DMR_filtered <- SNURF_TSS_DMR %>% semi_join(SNURF_ctrl_filtered, by = "start_position")

hp1_SNURF_TSS_DMR<- hp1 %>% filter(between(start_position,22691523,22693493))
hp1_SNURF_TSS_DMR_filtered <- hp1_SNURF_TSS_DMR %>% semi_join(SNURF_ctrl_filtered, by = "start_position")

hp2_SNURF_TSS_DMR<- hp2 %>% filter(between(start_position,22691523,22693493))
hp2_SNURF_TSS_DMR_filtered <- hp2_SNURF_TSS_DMR %>% semi_join(SNURF_ctrl_filtered, by = "start_position")

PEG3_TSS_DMR <- sorted_df %>% filter(between(start_position,59932193,59936667)) #MS-MLPA probe
PEG3_ctrl <- ctrls%>% filter(between(start_position,59932193,59936667))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
PEG3_ctrl_filtered <- PEG3_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
PEG3_TSS_DMR_filtered <- PEG3_TSS_DMR %>% semi_join(PEG3_ctrl_filtered, by = "start_position")

hp1_PEG3_TSS_DMR <- hp1 %>% filter(between(start_position,59932193,59936667))
hp1_PEG3_TSS_DMR_filtered <- hp1_PEG3_TSS_DMR %>% semi_join(PEG3_ctrl_filtered, by = "start_position")

hp2_PEG3_TSS_DMR <- hp2 %>% filter(between(start_position,59932193,59936667))
hp2_PEG3_TSS_DMR_filtered <- hp2_PEG3_TSS_DMR %>% semi_join(PEG3_ctrl_filtered, by = "start_position")

GNAS_NESP_TSS_DMR <- sorted_df %>% filter(between(start_position,60622124,60626697))#MS-MLPA probe
GNAS_NESP_ctrl <- ctrls %>% filter(between(start_position,60622124,60626697))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
GNAS_NESP_ctrl_filtered <- GNAS_NESP_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
GNAS_NESP_TSS_DMR_filtered <- GNAS_NESP_TSS_DMR %>% semi_join(GNAS_NESP_ctrl_filtered, by = "start_position")

hp1_GNAS_NESP_TSS_DMR <- hp1 %>% filter(between(start_position,60622124,60626697))
hp1_GNAS_NESP_TSS_DMR_filtered <- hp1_GNAS_NESP_TSS_DMR %>% semi_join(GNAS_NESP_ctrl_filtered, by = "start_position")

hp2_GNAS_NESP_TSS_DMR <- hp2 %>% filter(between(start_position,60622124,60626697))
hp2_GNAS_NESP_TSS_DMR_filtered <- hp2_GNAS_NESP_TSS_DMR %>% semi_join(GNAS_NESP_ctrl_filtered, by = "start_position")

GNAS_AS1_TSS_DMR <- sorted_df %>% filter(between(start_position,60633713,60636097))#MS-MLPA probe
GNAS_AS_ctrl <- ctrls %>% filter(between(start_position,60633713,60636097))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
GNAS_AS_ctrl_filtered <- GNAS_AS_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
GNAS_AS1_TSS_DMR_filtered <- GNAS_AS1_TSS_DMR %>% semi_join(GNAS_AS_ctrl_filtered, by = "start_position")

hp1_GNAS_AS1_TSS_DMR <- hp1 %>% filter(between(start_position,60633713,60636097))
hp1_GNAS_AS1_TSS_DMR_filtered <- hp1_GNAS_AS1_TSS_DMR %>% semi_join(GNAS_AS_ctrl_filtered, by = "start_position")

hp2_GNAS_AS1_TSS_DMR <- hp2 %>% filter(between(start_position,60633713,60636097))
hp2_GNAS_AS1_TSS_DMR_filtered <- hp2_GNAS_AS1_TSS_DMR %>% semi_join(GNAS_AS_ctrl_filtered, by = "start_position")

GNASXL_Ex1_DMR <- sorted_df %>% filter(between(start_position,60637080,60639527))#MS-MLPA probe
GNASXL_ctrl <- ctrls %>% filter(between(start_position,60637080,60639527))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
GNASXL_ctrl_filtered <- GNASXL_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
GNASXL_Ex1_DMR_filtered <- GNASXL_Ex1_DMR %>% semi_join(GNASXL_ctrl_filtered, by = "start_position")

hp1_GNASXL_Ex1_DMR <- hp1 %>% filter(between(start_position,60637080,60639527))
hp1_GNASXL_Ex1_DMR_filtered <- hp1_GNASXL_Ex1_DMR %>% semi_join(GNASXL_ctrl_filtered, by = "start_position")

hp2_GNASXL_Ex1_DMR <- hp2 %>% filter(between(start_position,60637080,60639527))
hp2_GNASXL_Ex1_DMR_filtered <- hp2_GNASXL_Ex1_DMR %>% semi_join(GNASXL_ctrl_filtered, by = "start_position")


GNAS_AB_TSS_DMR<- sorted_df %>% filter(between(start_position,60671509,60673270))#MS-MLPA probe
GNASAB_ctrl <- ctrls %>% filter(between(start_position,60671509,60673270))
#filter ctrl dataframe to rows where mean methylation frequency is between 40 and 60
GNASAB_ctrl_filtered <- GNASAB_ctrl %>% filter(between(methylation_frequency, 40, 60))
#filter dataframe to only include positions where the ctrl mean methylation is between 40 and 60 since these are most likely "imprinted" CpG positions
GNAS_AB_TSS_DMR_filtered <- GNAS_AB_TSS_DMR %>% semi_join(GNASAB_ctrl_filtered, by = "start_position")

hp1_GNAS_AB_TSS_DMR<- hp1 %>% filter(between(start_position,60671509,60673270))
hp1_GNAS_AB_TSS_DMR_filtered <- hp1_GNAS_AB_TSS_DMR %>% semi_join(GNASAB_ctrl_filtered, by = "start_position")

hp2_GNAS_AB_TSS_DMR<- hp2 %>% filter(between(start_position,60671509,60673270))
hp2_GNAS_AB_TSS_DMR_filtered <- hp2_GNAS_AB_TSS_DMR %>% semi_join(GNASAB_ctrl_filtered, by = "start_position")


#make table

#mean methylation is only calculated from positions where the controls have a mean methylation between 40 and 60
df<- data.frame(
  ID = c("TNDM","SRS","SRS","SRS/BWS", "SRS/BWS", "TS14/KOS14", "TS14/KOS14", "PWS/AS", "PWS/AS","MLID", "PHP", "PHP", "PHP", "PHP"),
  Chr = c("6q24","7p12","7q32","11p15", "11p15","14q32","14q32", "15q11", "15q11","19q13", "20q13","20q13","20q13","20q13"),
  DMR.official.name = c("PLAGL1:alt-TSS-DMR","GRB10:alt-TSS-DMR","MEST:alt-TSS-DMR"
                        ,"H19/IGF2:IG-DMR", "KCNQ1OT1:TSS-DMR","MEG3:TSS-DMR","MEG8:Int2-DMR",
                        "MAGEL2:TSS-DMR", "SNURF:TSS-DMR","PEG3:TSS-DMR","GNAS-NESP:TSS-DMR","GNAS-AS1:TSS-DMR","GNASXL:Ex1-DMR","GNAS A/B:TSS-DMR"),
  Genomic.coordinates = c("145199654-145201464",
                          "50942206-50944802",
                          "131804002-131808268",
                          "2085259-2091187",
                          "2788057-2790367",
                          "95059075-95062529",
                          "95139836-95140514",
                          "21381699-21383303",
                          "22691523-22693495",
                          "59932193-59936971",
                          "60622124-60626697",
                          "60633713-60636097",
                          "60636969-60639527",
                          "60671329-60673270"),
  CpGs = c(nrow(PLAGL1_alt_TSS_DMR_filtered), nrow(GRB10_alt_TSS_DMR_filtered), nrow(MEST_alt_TSS_DMR_filtered),
                    nrow(H19_IGF2_IG_DMR_filtered),nrow(KCNQ1OT1_TSS_DMR_filtered),
                    nrow(MEG3_TSS_DMR_filtered),nrow(MEG8_Int2_DMR_filtered), nrow(MAGEL2_TSS_DMR_filtered),
                    nrow(SNURF_TSS_DMR_filtered),nrow(PEG3_TSS_DMR_filtered),nrow(GNAS_NESP_TSS_DMR_filtered),
                    nrow(GNAS_AS1_TSS_DMR_filtered),nrow(GNASXL_Ex1_DMR_filtered), nrow(GNAS_AB_TSS_DMR_filtered)),
  Mean.coverage= c(mean(PLAGL1_alt_TSS_DMR$coverage),
                    mean(GRB10_alt_TSS_DMR$coverage),
                    mean(MEST_alt_TSS_DMR$coverage),
                    mean(H19_IGF2_IG_DMR$coverage),
                    mean(KCNQ1OT1_TSS_DMR$coverage),
                    mean(MEG3_TSS_DMR$coverage),
                    mean(MEG8_Int2_DMR$coverage),
                    mean(MAGEL2_TSS_DMR$coverage),
                    mean(SNURF_TSS_DMR$coverage),
                    mean(PEG3_TSS_DMR$coverage),
                    mean(GNAS_NESP_TSS_DMR$coverage),
                    mean(GNAS_AS1_TSS_DMR$coverage),
                    mean(GNASXL_Ex1_DMR$coverage),
                    mean(GNAS_AB_TSS_DMR$coverage)),
  hp1.mean.coverage= c(mean(hp1_PLAGL1_alt_TSS_DMR$coverage),
                    mean(hp1_GRB10_alt_TSS_DMR$coverage),
                    mean(hp1_MEST_alt_TSS_DMR$coverage),
                    mean(hp1_H19_IGF2_IG_DMR$coverage),
                    mean(hp1_KCNQ1OT1_TSS_DMR$coverage),
                    mean(hp1_MEG3_TSS_DMR$coverage),
                    mean(hp1_MEG8_Int2_DMR$coverage),
                    mean(hp1_MAGEL2_TSS_DMR$coverage),
                    mean(hp1_SNURF_TSS_DMR$coverage),
                    mean(hp1_PEG3_TSS_DMR$coverage),
                    mean(hp1_GNAS_NESP_TSS_DMR$coverage),
                    mean(hp1_GNAS_AS1_TSS_DMR$coverage),
                    mean(hp1_GNASXL_Ex1_DMR$coverage),
                    mean(hp1_GNAS_AB_TSS_DMR$coverage)),
  hp2.mean.coverage= c(mean(hp2_PLAGL1_alt_TSS_DMR$coverage),
                    mean(hp2_GRB10_alt_TSS_DMR$coverage),
                    mean(hp2_MEST_alt_TSS_DMR$coverage),
                    mean(hp2_H19_IGF2_IG_DMR$coverage),
                    mean(hp2_KCNQ1OT1_TSS_DMR$coverage),
                    mean(hp2_MEG3_TSS_DMR$coverage),
                    mean(hp2_MEG8_Int2_DMR$coverage),
                    mean(hp2_MAGEL2_TSS_DMR$coverage),
                    mean(hp2_SNURF_TSS_DMR$coverage),
                    mean(hp2_PEG3_TSS_DMR$coverage),
                    mean(hp2_GNAS_NESP_TSS_DMR$coverage),
                    mean(hp2_GNAS_AS1_TSS_DMR$coverage),
                    mean(hp2_GNASXL_Ex1_DMR$coverage),
                    mean(hp2_GNAS_AB_TSS_DMR$coverage)),
  Mean.methylation= c(mean(PLAGL1_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(GRB10_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(MEST_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(H19_IGF2_IG_DMR_filtered$methylation_frequency),
                        mean(KCNQ1OT1_TSS_DMR_filtered$methylation_frequency),
                        mean(MEG3_TSS_DMR_filtered$methylation_frequency),
                        mean(MEG8_Int2_DMR_filtered$methylation_frequency),
                        mean(MAGEL2_TSS_DMR_filtered$methylation_frequency),
                        mean(SNURF_TSS_DMR_filtered$methylation_frequency),
                        mean(PEG3_TSS_DMR_filtered$methylation_frequency),
                        mean(GNAS_NESP_TSS_DMR_filtered$methylation_frequency),
                        mean(GNAS_AS1_TSS_DMR_filtered$methylation_frequency),
                        mean(GNASXL_Ex1_DMR_filtered$methylation_frequency),
                        mean(GNAS_AB_TSS_DMR_filtered$methylation_frequency)),
  Mean.hp1=c(mean(hp1_PLAGL1_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_GRB10_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_MEST_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_H19_IGF2_IG_DMR_filtered$methylation_frequency),
                        mean(hp1_KCNQ1OT1_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_MEG3_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_MEG8_Int2_DMR_filtered$methylation_frequency),
                        mean(hp1_MAGEL2_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_SNURF_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_PEG3_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_GNAS_NESP_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_GNAS_AS1_TSS_DMR_filtered$methylation_frequency),
                        mean(hp1_GNASXL_Ex1_DMR_filtered$methylation_frequency),
                        mean(hp1_GNAS_AB_TSS_DMR_filtered$methylation_frequency)),
  Mean.hp2=c(mean(hp2_PLAGL1_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_GRB10_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_MEST_alt_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_H19_IGF2_IG_DMR_filtered$methylation_frequency),
                        mean(hp2_KCNQ1OT1_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_MEG3_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_MEG8_Int2_DMR_filtered$methylation_frequency),
                        mean(hp2_MAGEL2_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_SNURF_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_PEG3_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_GNAS_NESP_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_GNAS_AS1_TSS_DMR_filtered$methylation_frequency),
                        mean(hp2_GNASXL_Ex1_DMR_filtered$methylation_frequency),
                        mean(hp2_GNAS_AB_TSS_DMR_filtered$methylation_frequency)))


colnames(df) <- c("Imprinting disorder","Chr", "DMR Official name", "Genomic coordinates", "Number CpGs","Mean coverage","HP1 mean coverage",  "HP2 mean coverage", "Mean methylation (%)", "HP1 mean methylation (%)","HP2 mean methylation (%)")

df$"Mean coverage" <- round(df$"Mean coverage", 0)
df$"HP1 mean coverage" <- round(df$"HP1 mean coverage", 0)
df$"HP2 mean coverage" <- round(df$"HP2 mean coverage", 0)
df$"Mean methylation (%)" <- round(df$"Mean methylation", 1)
df$"HP1 mean methylation (%)" <- round(df$"HP1 mean methylation", 1)
df$"HP2 mean methylation (%)" <- round(df$"HP2 mean methylation", 1)



df <- df %>%
  mutate(across(everything(), as.character)) %>%
  kable("html", escape = FALSE) %>%
  kable_styling() %>%
  column_spec(9, color = ifelse(df$"Mean methylation" < 30 | df$"Mean methylation" > 75, "red", "black"),
              bold = ifelse(df$"Mean methylation" < 30 | df$"Mean methylation" > 75, TRUE, FALSE))


```

```{r echo=FALSE}
df 

```


```{r include=FALSE}
plagl1 <- ggplot(data = PLAGL1_alt_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=PLAGL1_ctrl, aes(ymin=methylation_frequency - sd, 
                                    ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=PLAGL1_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                    ymax= methylation_frequency + 0.1*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr6 position") +
  ggtitle ("Unphased: PLAGL1:alt-TSS-DMR") +
  coord_cartesian(xlim= c(145200200,145201432), ylim = c(-5,105), expand=F)+
  theme_light()+
  theme(text = element_text(size = 15))


grb10 <- ggplot(data = GRB10_alt_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=GRB10_ctrl, aes(ymin=methylation_frequency - sd, 
                                    ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=GRB10_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                    ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr7 position") +
  ggtitle ("Unphased: GRB10:alt-TSS-DMR") +
  coord_cartesian(xlim= c(50943355,50944182), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

mest <- ggplot(data = MEST_alt_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=MEST_ctrl, aes(ymin=methylation_frequency - sd, 
                                    ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=MEST_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                    ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr7 position") +
  ggtitle ("Unphased: MEST:alt-TSS-DMR") +
  coord_cartesian(xlim= c(131804587,131807000), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

h19 <- ggplot(data = H19_IGF2_IG_DMR,aes(x= start_position, y=methylation_frequency, group = 1)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=H19_ctrl, aes(ymin=methylation_frequency - sd, 
                                 ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=H19_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                 ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr11 position") +
  ggtitle ("Unphased: H19/IGF2:IG-DMR") +
  coord_cartesian(xlim= c(2085438,2090906), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

kcnq1 <- ggplot(data = KCNQ1OT1_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=KCNQ_ctrl, aes(ymin=methylation_frequency - sd, 
                                  ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=KCNQ_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                  ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr11 position") +
  ggtitle ("Unphased: KCNQ1OT1:TSS-DMR")+
  coord_cartesian(xlim= c(2788106,2790240), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))


meg3 <- ggplot(data = MEG3_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=MEG3_ctrl, aes(ymin=methylation_frequency - sd, 
                                  ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=MEG3_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                  ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr14 position") +
  ggtitle ("Unphased: MEG3:TSS-DMR")+
  coord_cartesian(xlim= c(95059075,95062529), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

meg8 <- ggplot(data = MEG8_Int2_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=MEG8_ctrl, aes(ymin=methylation_frequency - sd, 
                                  ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=MEG8_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                  ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr14 position") +
  ggtitle ("Unphased: MEG8:Int2-DMR")+
  coord_cartesian(xlim= c(95140029,95140384), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

magel2 <- ggplot(data = MAGEL2_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=MAGEL2_ctrl, aes(ymin=methylation_frequency - sd, 
                                  ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=MAGEL2_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                  ymax= methylation_frequency + 1.5*sd), alpha=0.1)+
  ylab("Methylation frequency (%)") +
  xlab("Chr15 position") +
  ggtitle ("Unphased: MAGEL2:TSS-DMR")+
  coord_cartesian(xlim= c(21381784,21383056), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))


snurf <- ggplot(data = SNURF_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=SNURF_ctrl, aes(ymin=methylation_frequency - sd, 
                                  ymax= methylation_frequency + sd), alpha=0.2)+
  
  geom_ribbon(data=SNURF_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                   ymax= methylation_frequency + 1.5*sd), alpha=0.1)+ 
  ylab("Methylation frequency (%)") +
  xlab("Chr15 position") +
  ggtitle ("Unphased: SNURF:TSS-DMR")+
  coord_cartesian(xlim= c(22691523,22693493), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

peg3 <- ggplot(data = PEG3_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=PEG3_ctrl, aes(ymin=methylation_frequency - sd, 
                                   ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=PEG3_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                   ymax= methylation_frequency + 1.5*sd), alpha=0.1)+ 
  ylab("Methylation frequency (%)") +
  xlab("Chr19 position") +
  ggtitle ("Unphased: PEG3:TSS-DMR")+
  coord_cartesian(xlim= c(59932193,59936667), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

gnasnesp <- ggplot(data = GNAS_NESP_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=GNAS_NESP_ctrl, aes(ymin=methylation_frequency - sd, 
                                          ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=GNAS_NESP_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                          ymax= methylation_frequency + 1.5*sd), alpha=0.1)+ 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Unphased: GNAS-NESP:TSS-DMR")+
  coord_cartesian(xlim= c(60622124,60626697), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

gnas1 <- ggplot(data = GNAS_AS1_TSS_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=GNAS_AS_ctrl, aes(ymin=methylation_frequency - sd, 
                                          ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=GNAS_AS_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                          ymax= methylation_frequency + 1.5*sd), alpha=0.1)+ 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Unphased: GNAS-AS1:TSS-DMR")+
  coord_cartesian(xlim= c(60633713,60636097), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

gnasxl <- ggplot(data = GNASXL_Ex1_DMR,aes(x= start_position, group = 1, se=F)) + 
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=GNASXL_ctrl, aes(ymin=methylation_frequency - sd, 
                                          ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=GNASXL_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                          ymax= methylation_frequency + 1.5*sd), alpha=0.1)+ 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Unphased: GNASXL:Ex1-DMR")+
  coord_cartesian(xlim= c(60637080,60639527), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

gnasab <- ggplot(data = GNAS_AB_TSS_DMR,aes(x= start_position, group = 1, se=F)) +
  geom_smooth(aes(y=methylation_frequency), colour="black", se=F) + 
  geom_point(aes(y = methylation_frequency), colour = "black", alpha=0.3) +
  geom_ribbon(data=GNASAB_ctrl, aes(ymin=methylation_frequency - sd, 
                                          ymax= methylation_frequency + sd), alpha=0.2)+
  geom_ribbon(data=GNASAB_ctrl, aes(ymin=methylation_frequency - 1.5*sd, 
                                          ymax= methylation_frequency + 1.5*sd), alpha=0.1)+ 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Unphased: GNAS A/B:TSS-DMR")+
  coord_cartesian(xlim= c(60671509,60673270), ylim = c(-5,105), expand=F)+
  theme_light()+
    theme(text = element_text(size = 15))

# Check if any data frames are empty
if(nrow(hp1_PLAGL1_alt_TSS_DMR) > 0 && nrow(hp2_PLAGL1_alt_TSS_DMR) > 0) {
  # Create the plot
  plagl1_hp <- ggplot() + 
    geom_smooth(data = hp1_PLAGL1_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
    geom_point(data = hp1_PLAGL1_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
    geom_smooth(data = hp2_PLAGL1_alt_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
    geom_point(data = hp2_PLAGL1_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
    ylab("Methylation frequency (%)") +
    xlab("Chr6 position") +
    ggtitle("Phased: PLAGL1:alt-TSS-DMR ") +
    coord_cartesian(xlim = c(145200200, 145201432), ylim = c(-5, 105), expand = F) +
    theme_light()+
    theme(text = element_text(size = 15))
} else {
# Create the plot without data
  plagl1_hp <- ggplot() + 
    ylab("Methylation frequency (%)") +
    xlab("Chr6 position") +
    ggtitle("Phased: PLAGL1:alt-TSS-DMR ") +
    coord_cartesian(xlim = c(145200200, 145201432), ylim = c(-5, 105), expand = F) +
    theme_light() +
    theme(text = element_text(size = 15)) +
    geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_GRB10_alt_TSS_DMR) > 0 && nrow(hp2_GRB10_alt_TSS_DMR) > 0) {
grb10_hp <- ggplot() + 
  geom_smooth(data = hp1_GRB10_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_GRB10_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_GRB10_alt_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_GRB10_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr7 position") +
  ggtitle ("Phased: GRB10:alt-TSS-DMR") +
  coord_cartesian(xlim= c(50943355,50944182), ylim = c(-5,105), expand=F)+
  theme_light()+
  theme(text = element_text(size = 15))
} else {
# Create the plot without data
  grb10_hp <- ggplot() + 
    ylab("Methylation frequency (%)") +
  xlab("Chr7 position") +
  ggtitle ("Phased: GRB10:alt-TSS-DMR") +
  coord_cartesian(xlim= c(50943355,50944182), ylim = c(-5,105), expand=F)+
  theme_light() +
  theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_MEST_alt_TSS_DMR) > 0 && nrow(hp2_MEST_alt_TSS_DMR) > 0) {
mest_hp <- ggplot() + 
  geom_smooth(data = hp1_MEST_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_MEST_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_MEST_alt_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_MEST_alt_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr7 position") +
  ggtitle ("Phased: MEST:alt-TSS-DMR") +
  coord_cartesian(xlim= c(131804587,131807000), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  mest_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr7 position") +
  ggtitle ("Phased: MEST:alt-TSS-DMR") +
  coord_cartesian(xlim= c(131804587,131807000), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_H19_IGF2_IG_DMR) > 0 && nrow(hp2_H19_IGF2_IG_DMR) > 0) {
h19_hp <- ggplot() + 
  geom_smooth(data = hp1_H19_IGF2_IG_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_H19_IGF2_IG_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_H19_IGF2_IG_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_H19_IGF2_IG_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr11 position") +
  ggtitle ("Phased: H19/IGF2:IG-DMR") +
  coord_cartesian(xlim= c(2085438,2090906), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  h19_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr11 position") +
  ggtitle ("Phased: H19/IGF2:IG-DMR") +
  coord_cartesian(xlim= c(2085438,2090906), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_KCNQ1OT1_TSS_DMR) > 0 && nrow(hp2_KCNQ1OT1_TSS_DMR) > 0) {
kcnq1_hp <- ggplot() + 
  geom_smooth(data = hp1_KCNQ1OT1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_KCNQ1OT1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_KCNQ1OT1_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_KCNQ1OT1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr11 position") +
  ggtitle ("Phased: KCNQ1OT1:TSS-DMR")+
  coord_cartesian(xlim= c(2788106,2790240), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  kcnq1_hp <- ggplot() + 
ylab("Methylation frequency (%)") +
  xlab("Chr11 position") +
  ggtitle ("Phased: KCNQ1OT1:TSS-DMR")+
  coord_cartesian(xlim= c(2788106,2790240), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_MEG3_TSS_DMR) > 0 && nrow(hp2_MEG3_TSS_DMR) > 0) {
meg3_hp <- ggplot() + 
  geom_smooth(data = hp1_MEG3_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_MEG3_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_MEG3_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_MEG3_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr14 position") +
  ggtitle ("Phased: MEG3:TSS-DMR")+
  coord_cartesian(xlim= c(95059075,95062529), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  meg3_hp <- ggplot() + 
 ylab("Methylation frequency (%)") +
  xlab("Chr14 position") +
  ggtitle ("Phased: MEG3:TSS-DMR")+
  coord_cartesian(xlim= c(95059075,95062529), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}


# Check if any data frames are empty
if(nrow(hp1_MEG8_Int2_DMR) > 0 && nrow(hp2_MEG8_Int2_DMR) > 0) {
meg8_hp <- ggplot() + 
  geom_smooth(data = hp1_MEG8_Int2_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_MEG8_Int2_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_MEG8_Int2_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_MEG8_Int2_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr14 position") +
  ggtitle ("Phased: MEG8:Int2-DMR")+
  coord_cartesian(xlim= c(95140029,95140384), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  meg8_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr14 position") +
  ggtitle ("Phased: MEG8:Int2-DMR")+
  coord_cartesian(xlim= c(95140029,95140384), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_MAGEL2_TSS_DMR) > 0 && nrow(hp2_MAGEL2_TSS_DMR) > 0) {
magel2_hp <- ggplot() + 
  geom_smooth(data = hp1_MAGEL2_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_MAGEL2_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_MAGEL2_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_MAGEL2_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr15 position") +
  ggtitle ("Phased: MAGEL2:TSS-DMR")+
  coord_cartesian(xlim= c(21381784,21383056), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  magel2_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr15 position") +
  ggtitle ("Phased: MAGEL2:TSS-DMR")+
  coord_cartesian(xlim= c(21381784,21383056), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_SNURF_TSS_DMR) > 0 && nrow(hp2_SNURF_TSS_DMR) > 0) {
snurf_hp <- ggplot() + 
  geom_smooth(data = hp1_SNURF_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_SNURF_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_SNURF_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_SNURF_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr15 position") +
  ggtitle ("Phased: SNURF:TSS-DMR")+
  coord_cartesian(xlim= c(22691523,22693493), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  snurf_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr15 position") +
  ggtitle ("Phased: SNURF:TSS-DMR")+
  coord_cartesian(xlim= c(22691523,22693493), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_PEG3_TSS_DMR) > 0 && nrow(hp2_PEG3_TSS_DMR) > 0) {
peg3_hp <- ggplot() + 
  geom_smooth(data = hp1_PEG3_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_PEG3_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_PEG3_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_PEG3_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr19 position") +
  ggtitle ("Phased: PEG3:TSS-DMR")+
  coord_cartesian(xlim= c(59932193,59936667), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  peg3_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr19 position") +
  ggtitle ("Phased: PEG3:TSS-DMR")+
  coord_cartesian(xlim= c(59932193,59936667), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_GNAS_NESP_TSS_DMR) > 0 && nrow(hp2_GNAS_NESP_TSS_DMR) > 0) {
gnasnesp_hp <- ggplot() + 
  geom_smooth(data = hp1_GNAS_NESP_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_GNAS_NESP_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_GNAS_NESP_TSS_DMR , aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_GNAS_NESP_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNAS-NESP:TSS-DMR")+
  coord_cartesian(xlim= c(60622124,60626697), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  gnasnesp_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNAS-NESP:TSS-DMR")+
  coord_cartesian(xlim= c(60622124,60626697), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_GNAS_AS1_TSS_DMR ) > 0 && nrow(hp2_GNAS_AS1_TSS_DMR) > 0) {
gnas1_hp <- ggplot() + 
  geom_smooth(data = hp1_GNAS_AS1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_GNAS_AS1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_GNAS_AS1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_GNAS_AS1_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNAS-AS1:TSS-DMR")+
  coord_cartesian(xlim= c(60633713,60636097), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  gnas1_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNAS-AS1:TSS-DMR")+
  coord_cartesian(xlim= c(60633713,60636097), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_GNASXL_Ex1_DMR ) > 0 && nrow(hp2_GNASXL_Ex1_DMR) > 0) {
gnasxl_hp <- ggplot() + 
  geom_smooth(data = hp1_GNASXL_Ex1_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_GNASXL_Ex1_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_GNASXL_Ex1_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_GNASXL_Ex1_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNASXL:Ex1-DMR")+
  coord_cartesian(xlim= c(60637080,60639527), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  gnasxl_hp <- ggplot() + 
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNASXL:Ex1-DMR")+
  coord_cartesian(xlim= c(60637080,60639527), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}

# Check if any data frames are empty
if(nrow(hp1_GNAS_AB_TSS_DMR ) > 0 && nrow(hp2_GNAS_AB_TSS_DMR) > 0) {
gnasab_hp <- ggplot() + 
  geom_smooth(data = hp1_GNAS_AB_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", se = F) + 
  geom_point(data = hp1_GNAS_AB_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "red", alpha = 0.3) +
  geom_smooth(data = hp2_GNAS_AB_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", se = F) + 
  geom_point(data = hp2_GNAS_AB_TSS_DMR, aes(x = start_position, y = methylation_frequency), colour = "blue", alpha = 0.3) +
  ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNAS A/B:TSS-DMR")+
  coord_cartesian(xlim= c(60671509,60673270), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15))
} else {
# Create the plot without data
  gnasab_hp <- ggplot() + 
   ylab("Methylation frequency (%)") +
  xlab("Chr20 position") +
  ggtitle ("Phased: GNAS A/B:TSS-DMR")+
  coord_cartesian(xlim= c(60671509,60673270), ylim = c(-5,105), expand=F)+
  theme_light()+
theme(text = element_text(size = 15)) +
  geom_blank()
}



```

## Visualization of methylation at imprinted regions
In the left side each imprinted region is shown with methylation frequency (%). Grey shaded regions mark mean $\pm$ 1 SD (darker grey) and $\pm$ 1.5 SD (lighter grey) for control samples. Haplotype phased methylation frequencies for the regions are seen to the right. HP1 is **<span style="color:red">red</span>** and HP2 is **<span style="color:blue">blue</span>**. If plots are empty phasing was not successful. 


### TNDM
```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
plagl1
plagl1_hp
```

### SRS
```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
grb10
grb10_hp
mest
mest_hp
```

### SRS/BWS
```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
h19
h19_hp
kcnq1
kcnq1_hp
```

### TS14/KOS14
```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
meg3
meg3_hp
meg8
meg8_hp
```

### PWS/AS
```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
magel2
magel2_hp
snurf
snurf_hp
```

### MLID
```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
peg3
peg3_hp
```

### PHP

```{r, fig.show="hold", out.width="50%", echo=FALSE, message=FALSE}
gnasnesp
gnasnesp_hp
gnas1
gnas1_hp
gnasxl
gnasxl_hp
gnasab
gnasab_hp
```