1. Environment setting

library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(scales)
library(gridExtra)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(dmrseq)
## Loading required package: bsseq
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:tidyr':
## 
##     expand
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:plotly':
## 
##     rename
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## The following object is masked from 'package:plotly':
## 
##     slice
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
library(annotatr)
library(sechm)
library(GenomicRanges)
source("~/Transgenerational/EMSeq/NewData_fromLaSitua/Methylation_helper.R")

2. Loading data

2.1 BSSeq objects

Loading of bsseq objects with CpGs shared by all samples.

bsseq_obj <- readRDS("/group/testa/Project/PGCLC/emseq_rebuttal/2.PreliminaryAnalysis/bsseq_obj_sharedby75ofall.rds") #getting bsseq object where sample selection and CpG filtering were already performed

2.2 Annotated DMRs

The DMRs here selected were called like with: delta=0.25 and p.threshold=0.0001

DMRs_annotated1 <- readRDS(params$DMRAnnotated1)
DMRs_annotated2 <- readRDS(params$DMRAnnotated2)
DMRs_annotated3 <- readRDS(params$DMRAnnotated3)
DMRs_annotated4 <- readRDS(params$DMRAnnotated4)

From DSS vignette

Q: What’s the meaning of areaStat parameter?

A: We adapt that from the bsseq package. It is the sum of the test statistics of all CpG sites within a DMR. It doesn’t have a direct biological meaning. One can imagine when we try to rank the DMRs, we don’t know whether the height or the width is more important? AreaStat is a combination of the two. This is an ad hoc way to rank DMRs, but larger AreaStat is more likely to be a DMR

3. Get annoTrack with annotatr and dmrseq

# get annotations for hg38
annoTrack <- dmrseq::getAnnot("hg38")
## Building CpG islands...
## Building CpG shores...
## Building CpG shelves...
## Building inter-CpG-islands...
## Download of CpG annotation successful!
## Loading required package: GenomicFeatures
## Loading required package: AnnotationDbi
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:plotly':
## 
##     select
## 
## 'select()' returned 1:1 mapping between keys and columns
## Building cds...
## Download of Gene annotation successful!

4. Exploration of DMR-associated genes

cellTypes_colors <- c(hiPSCs ="#dd1c77", iMeLCs ="#377eb8", hPGCLCs ="#4daf4a", hEGCLCs_p5 = "#ff7f00", hEGCLCs_p10 = "#ffc400")
colData(bsseq_obj)$Type <- factor(colData(bsseq_obj)$Type, levels = c("hiPSCs", "iMeLCs" , "hPGCLCs", "hEGCLCs_p5", "hEGCLCs_p10"))

4.1 hEGCLCs_p10 vs hiPSCs

comparison <- "hEGCLCs_p10vshiPSCs"
DMRs <- DMRs_annotated1[[comparison]]
DMRs$DMR <- paste0("DMR", rownames(DMRs))
DMRs_gr <- makeGRangesFromDataFrame(DMRs, keep.extra.columns = TRUE)
m <- findOverlaps(rowRanges(bsseq_obj), DMRs_gr)

bsseq_obj_subset <- bsseq_obj[queryHits(m),]

rowData(bsseq_obj_subset)$DMR <- DMRs_gr[subjectHits(m), ]$DMR
rowData(bsseq_obj_subset)$GenomicAnn <- DMRs_gr[subjectHits(m), ]$ChipSeekerAnn
rowData(bsseq_obj_subset)$Gene <- DMRs_gr[subjectHits(m), ]$hgnc_symbol
rowData(bsseq_obj_subset)$GeneBiotype <- DMRs_gr[subjectHits(m),]$gene_biotype
rowData(bsseq_obj_subset)$GenomicRegion <- rowData(bsseq_obj_subset)$GenomicAnn

rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Exon")] <- "Exon"
rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Intron")] <- "Intron"

In total we have 5407 CpGs in the significant DMRs.

Calculate average methylation of CpGs in each DMR

assays(bsseq_obj_subset)$Meth <- getMeth(bsseq_obj_subset, type = "raw")

AvgMeth <- aggregate(assays(bsseq_obj_subset)$Meth, list(rowData(bsseq_obj_subset)$DMR), mean, na.rm = T)
rownames(AvgMeth) <- AvgMeth$Group.1
AvgMeth <- AvgMeth[, -1]

rowDataSE <- rowData(bsseq_obj_subset)
rownames(rowDataSE) <- rowData(bsseq_obj_subset)$DMR
rowDataSE <- rowDataSE[rownames(AvgMeth), ] 
SE <- SummarizedExperiment(assays=AvgMeth, colData=colData(bsseq_obj_subset), rowData = rowDataSE, )

assays(SE) <- assays(SE) %>% setNames("AvgMeth")

identical(rownames(assay(SE)), rownames(rowData(SE))) #TRUE
## [1] TRUE
## [1] TRUE

Setting colors

ScaledCols <- c('darkblue', "purple", "white", "lightgoldenrod1", 'goldenrod1')

metadata(SE)$anno_colors <- list(Type = cellTypes_colors,
                                               GenomicRegion = c ("3' UTR" = "#FFADAD", "5' UTR" = "#FFD6A5", "Distal Intergenic" = "#FDFFB6", "Downstream (<=300bp)" = "#FFFFFC", "Exon" = "#9BF6FF", "Intron" = "#A0C4FF", "Promoter (<=1kb)" = "#BDB2FF", "Promoter (1-2kb)" = "#FFC6FF", "Promoter (2-3kb)" = "#CAFFBF"))

sample_colors <- c("hiPSC_rep1" = "#dd1c77", "hiPSC_rep2" = "#dd1c77", "hiPSC_rep3" = "#dd1c77", "hiPSC_rep4" = "#dd1c77", "hiPSC_rep5" = "#dd1c77",
                   "iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
                   "hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
                   "hEGCLC_p5_rep1" = "#ff7f00", "hEGCLC_p5_rep2" = "#ff7f00", "hEGCLC_p5_rep3" = "#ff7f00",
                    "hEGCLC_p10_rep1" = "#ffc400", "hEGCLC_p10_rep2" = "#ffc400", "hEGCLC_p10_rep3" = "#ffc400")

4.1.1 Heatmaps

Without sample clustering

break_coord <- assay(SE) %>% quantile(seq(0,1, by=0.1), na.rm =TRUE)
ExpCols <- viridisLite::plasma(length(seq(head(break_coord,1), tail(break_coord,1), by=0.2)))

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hiPSCs", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0)

With sample clustering

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, cluster_cols = TRUE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hiPSCs", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0)

4.1.2 Violin plots

Dividing by sample

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hiPSCs", fill = "Genomic region")
## Warning: Removed 8 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning: Removed 8 rows containing non-finite values (`stat_summary()`).

Dividing by genomic region

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hiPSCs", fill = "Sample", colors = sample_colors)
## Warning: Removed 8 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
## Warning: Removed 8 rows containing non-finite values (`stat_summary()`).

4.1.3 Top 20 DMR-associated genes

Useful to visualize single DMR with exon and CpG annotations. Here the top 20 significant DMRs are plotted. Rows in DMRs are sorted by “areaStat”, which is the sum of test statistics of all CpG sites in the region. Remember that we do not have p-value for DMRs because in DMR calling, the statistical test is perform for each CpG, and then the significant CpG are merged into regions. It is very difficult to estimate FDR at the region-level, based on the site level test results (p-values or FDR). This is a difficult question for all types of genome-wide assays such as the ChIP-seq.

Here we can filter for only differentially methylated promoters of protein coding genes

DMRs <- DMRs[grep("Promoter", DMRs$ChipSeekerAnn),] %>% filter(gene_biotype == "protein_coding")
nr_of_top <- 20

if (nrow(DMRs) > 20){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs)}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[i,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 500, main = paste0(DMRs[i,]$hgnc_symbol, " - ", DMRs[i,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: HTRA4, PLAGL1, GNAS, CTSF, SLITRK2, COMT, KHDRBS2, TRIM58, ZSCAN30, LRRN4CL, C5orf63, ARHGEF4, GSTM5, ITGBL1, MARCKS, CD1D, ZNF490, CDH24, KLHL40, PLEKHG5

4.1.3 Top 10 DMR-associated genes where there is gain in methylation

Here the top 10 DMRs where methylation is higher in hEGCLCs_p10 (and lower in hiPSCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy<0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy<0, ])}

for(i in 1:nr_plots){
  plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy<0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy<0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy<0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: HTRA4, PLAGL1, GNAS, CTSF, COMT, KHDRBS2, TRIM58, ZSCAN30, LRRN4CL, C5orf63

4.1.4 Top 10 DMR-associated genes where there is loss in methylation

Here the top 10 DMRs where methylation is higher in hiPSCs (and lower in hEGCLCs_p10).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy>0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy>0, ])}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy>0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy>0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy>0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: SLITRK2, CHCHD2, ZNF717, SLC15A4, PABPC4L, SLC30A8, HTR1A, RIMBP2, LYNX1-SLURP2, PRDM16

4.1.5 Exploration of PGC markers and other relevant genes

SelectedGenes <- c("PIWIL2", "NANOS3", "IGF1", "DNMT3L", "PIFO", "SOX15", "CSF3R", "NUMB", "CXCR4")
PresentGenes <- vector()
AbsentGenes <- vector()
for(gene_index in SelectedGenes){
  if (any(DMRs$hgnc_symbol %in% gene_index)) {
    PresentGenes <- c(PresentGenes, gene_index)
    plotDMRs(bsseq_obj, regions=DMRs[DMRs$hgnc_symbol %in% gene_index,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(gene_index," - ", DMRs[DMRs$hgnc_symbol %in% gene_index,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
  } else{ 
    AbsentGenes <- c(AbsentGenes, gene_index)
    next}
}

These genes were found: These genes were not found: PIWIL2, NANOS3, IGF1, DNMT3L, PIFO, SOX15, CSF3R, NUMB, CXCR4

4.2 iMeLCs vs hEGCLCs_p10

comparison <- "hEGCLCs_p10vsiMeLCs"
DMRs <- DMRs_annotated2[[comparison]]
DMRs$DMR <- paste0("DMR", rownames(DMRs))
DMRs_gr <- makeGRangesFromDataFrame(DMRs, keep.extra.columns = TRUE)
m <- findOverlaps(rowRanges(bsseq_obj), DMRs_gr)

bsseq_obj_subset <- bsseq_obj[queryHits(m),]

rowData(bsseq_obj_subset)$DMR <- DMRs_gr[subjectHits(m), ]$DMR
rowData(bsseq_obj_subset)$GenomicAnn <- DMRs_gr[subjectHits(m), ]$ChipSeekerAnn
rowData(bsseq_obj_subset)$Gene <- DMRs_gr[subjectHits(m), ]$hgnc_symbol
rowData(bsseq_obj_subset)$GeneBiotype <- DMRs_gr[subjectHits(m),]$gene_biotype
rowData(bsseq_obj_subset)$GenomicRegion <- rowData(bsseq_obj_subset)$GenomicAnn

rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Exon")] <- "Exon"
rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Intron")] <- "Intron"

In total we have 6127 CpGs in the significant DMRs.

Calculate average methylation of CpGs in each DMR

assays(bsseq_obj_subset)$Meth <- getMeth(bsseq_obj_subset, type = "raw")

AvgMeth <- aggregate(assays(bsseq_obj_subset)$Meth, list(rowData(bsseq_obj_subset)$DMR), mean, na.rm = T)
rownames(AvgMeth) <- AvgMeth$Group.1
AvgMeth <- AvgMeth[, -1]

rowDataSE <- rowData(bsseq_obj_subset)
rownames(rowDataSE) <- rowData(bsseq_obj_subset)$DMR
rowDataSE <- rowDataSE[rownames(AvgMeth), ] 
SE <- SummarizedExperiment(assays=AvgMeth, colData=colData(bsseq_obj_subset), rowData = rowDataSE, )

assays(SE) <- assays(SE) %>% setNames("AvgMeth")

identical(rownames(assay(SE)), rownames(rowData(SE))) #TRUE
## [1] TRUE
## [1] TRUE

Setting colors

ScaledCols <- c('darkblue', "purple", "white", "lightgoldenrod1", 'goldenrod1')

metadata(SE)$anno_colors <- list(Type = cellTypes_colors,
                                               GenomicRegion = c ("3' UTR" = "#FFADAD", "5' UTR" = "#FFD6A5", "Distal Intergenic" = "#FDFFB6", "Downstream (<=300bp)" = "#FFFFFC", "Exon" = "#9BF6FF", "Intron" = "#A0C4FF", "Promoter (<=1kb)" = "#BDB2FF", "Promoter (1-2kb)" = "#FFC6FF", "Promoter (2-3kb)" = "#CAFFBF"))

sample_colors <- c("hiPSC_rep1" = "#dd1c77", "hiPSC_rep2" = "#dd1c77", "hiPSC_rep3" = "#dd1c77", "hiPSC_rep4" = "#dd1c77", "hiPSC_rep5" = "#dd1c77",
                   "iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
                   "hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
                   "hEGCLC_p5_rep1" = "#ff7f00", "hEGCLC_p5_rep2" = "#ff7f00", "hEGCLC_p5_rep3" = "#ff7f00",
                    "hEGCLC_p10_rep1" = "#ffc400", "hEGCLC_p10_rep2" = "#ffc400", "hEGCLC_p10_rep3" = "#ffc400")

4.2.1 Heatmaps

Without sample clustering

break_coord <- assay(SE) %>% quantile(seq(0,1, by=0.1), na.rm =TRUE)
ExpCols <- viridisLite::plasma(length(seq(head(break_coord,1), tail(break_coord,1), by=0.2)))

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs iMeLCs", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0)

With sample clustering

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, cluster_cols = TRUE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs iMeLCs", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0)

4.2.2 Violin plots

Dividing by sample

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs iMeLCs", fill = "Genomic region")
## Warning: Removed 12 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning: Removed 12 rows containing non-finite values (`stat_summary()`).

Dividing by genomic region

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs iMeLCs", fill = "Sample", colors = sample_colors)
## Warning: Removed 12 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
## Warning: Removed 12 rows containing non-finite values (`stat_summary()`).

4.2.3 Top 20 DMR-associated genes

Useful to visualize single DMR with exon and CpG annotations. Here the top 20 significant DMRs are plotted. Rows in DMRs are sorted by “areaStat”, which is the sum of test statistics of all CpG sites in the region. Remember that we do not have p-value for DMRs because in DMR calling, the statistical test is perform for each CpG, and then the significant CpG are merged into regions. It is very difficult to estimate FDR at the region-level, based on the site level test results (p-values or FDR). This is a difficult question for all types of genome-wide assays such as the ChIP-seq.

Here we can filter for only differentially methylated promoters of protein coding genes

DMRs <- DMRs[grep("Promoter", DMRs$ChipSeekerAnn),] %>% filter(gene_biotype == "protein_coding")
nr_of_top <- 20

if (nrow(DMRs) > 20){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs)}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[i,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 500, main = paste0(DMRs[i,]$hgnc_symbol, " - ", DMRs[i,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: HTRA4, PLAGL1, GNAS, CTSF, SLITRK2, IFITM1, COMT, KHDRBS2, C5orf63, ZSCAN30, LRRN4CL, ARHGEF4, ZNF490, CHCHD2, GSTM5, CDH24, CD1D, PLEKHG5, KLHL40, TRIM58

4.2.3 Top 10 DMR-associated genes where there is gain in methylation

Here the top 10 DMRs where methylation is higher in hEGCLCs (and lower in hPGCLCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy<0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy<0, ])}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy<0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy<0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy<0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: HTRA4, PLAGL1, GNAS, CTSF, COMT, KHDRBS2, C5orf63, ZSCAN30, LRRN4CL, ARHGEF4

4.2.4 Top 10 DMR-associated genes where there is loss in methylation

Here the top 10 DMRs where methylation is higher in hPGCLCs (and lower in hEGCLCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy>0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy>0, ])}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy>0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy>0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy>0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: SLITRK2, IFITM1, CHCHD2, FZD5, SLC15A4, FAM83G, TMEM121, ZNF717, HTR1A, PABPC4L

4.2.5 Exploration of PGC markers and other relevant genes

SelectedGenes <- c("PIWIL2", "NANOS3", "IGF1", "DNMT3L", "PIFO", "SOX15", "CSF3R", "NUMB", "CXCR4")
PresentGenes <- vector()
AbsentGenes <- vector()
for(gene_index in SelectedGenes){
  if (any(DMRs$hgnc_symbol %in% gene_index)) {
    PresentGenes <- c(PresentGenes, gene_index)
    plotDMRs(bsseq_obj, regions=DMRs[DMRs$hgnc_symbol %in% gene_index,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(gene_index," - ", DMRs[DMRs$hgnc_symbol %in% gene_index,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
  } else{ 
    AbsentGenes <- c(AbsentGenes, gene_index)
    next}
}

These genes were found: SOX15 These genes were not found: PIWIL2, NANOS3, IGF1, DNMT3L, PIFO, CSF3R, NUMB, CXCR4

4.3 hPGCLCs vs hEGCLCs_p10

comparison <- "hEGCLCs_p10vshPGCLCs"
DMRs <- DMRs_annotated3[[comparison]]
DMRs$DMR <- paste0("DMR", rownames(DMRs))
DMRs_gr <- makeGRangesFromDataFrame(DMRs, keep.extra.columns = TRUE)
m <- findOverlaps(rowRanges(bsseq_obj), DMRs_gr)

bsseq_obj_subset <- bsseq_obj[queryHits(m),]

rowData(bsseq_obj_subset)$DMR <- DMRs_gr[subjectHits(m), ]$DMR
rowData(bsseq_obj_subset)$GenomicAnn <- DMRs_gr[subjectHits(m), ]$ChipSeekerAnn
rowData(bsseq_obj_subset)$Gene <- DMRs_gr[subjectHits(m), ]$hgnc_symbol
rowData(bsseq_obj_subset)$GeneBiotype <- DMRs_gr[subjectHits(m),]$gene_biotype
rowData(bsseq_obj_subset)$GenomicRegion <- rowData(bsseq_obj_subset)$GenomicAnn

rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Exon")] <- "Exon"
rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Intron")] <- "Intron"

In total we have 25651 CpGs in the significant DMRs.

Calculate average methylation of CpGs in each DMR

assays(bsseq_obj_subset)$Meth <- getMeth(bsseq_obj_subset, type = "raw")

AvgMeth <- aggregate(assays(bsseq_obj_subset)$Meth, list(rowData(bsseq_obj_subset)$DMR), mean, na.rm = T)
rownames(AvgMeth) <- AvgMeth$Group.1
AvgMeth <- AvgMeth[, -1]

rowDataSE <- rowData(bsseq_obj_subset)
rownames(rowDataSE) <- rowData(bsseq_obj_subset)$DMR
rowDataSE <- rowDataSE[rownames(AvgMeth), ] 
SE <- SummarizedExperiment(assays=AvgMeth, colData=colData(bsseq_obj_subset), rowData = rowDataSE, )

assays(SE) <- assays(SE) %>% setNames("AvgMeth")

identical(rownames(assay(SE)), rownames(rowData(SE))) #TRUE
## [1] TRUE
## [1] TRUE

Setting colors

ScaledCols <- c('darkblue', "purple", "white", "lightgoldenrod1", 'goldenrod1')

metadata(SE)$anno_colors <- list(Type = cellTypes_colors,
                                               GenomicRegion = c ("3' UTR" = "#FFADAD", "5' UTR" = "#FFD6A5", "Distal Intergenic" = "#FDFFB6", "Downstream (<=300bp)" = "#FFFFFC", "Exon" = "#9BF6FF", "Intron" = "#A0C4FF", "Promoter (<=1kb)" = "#BDB2FF", "Promoter (1-2kb)" = "#FFC6FF", "Promoter (2-3kb)" = "#CAFFBF"))

sample_colors <- c("hiPSC_rep1" = "#dd1c77", "hiPSC_rep2" = "#dd1c77", "hiPSC_rep3" = "#dd1c77", "hiPSC_rep4" = "#dd1c77", "hiPSC_rep5" = "#dd1c77",
                   "iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
                   "hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
                   "hEGCLC_p5_rep1" = "#ff7f00", "hEGCLC_p5_rep2" = "#ff7f00", "hEGCLC_p5_rep3" = "#ff7f00",
                    "hEGCLC_p10_rep1" = "#ffc400", "hEGCLC_p10_rep2" = "#ffc400", "hEGCLC_p10_rep3" = "#ffc400")

4.3.1 Heatmaps

Without sample clustering

break_coord <- assay(SE) %>% quantile(seq(0,1, by=0.1), na.rm =TRUE)
ExpCols <- viridisLite::plasma(length(seq(head(break_coord,1), tail(break_coord,1), by=0.2)))

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hPGCLCs", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0) 
## `use_raster` is automatically set to TRUE for a matrix with more than
## 2000 rows. You can control `use_raster` argument by explicitly setting
## TRUE/FALSE to it.
## 
## Set `ht_opt$message = FALSE` to turn off this message.
## 'magick' package is suggested to install to give better rasterization.
## 
## Set `ht_opt$message = FALSE` to turn off this message.

With sample clustering

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, cluster_cols = TRUE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hPGCLCs", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0) 
## `use_raster` is automatically set to TRUE for a matrix with more than
## 2000 rows. You can control `use_raster` argument by explicitly setting
## TRUE/FALSE to it.
## 
## Set `ht_opt$message = FALSE` to turn off this message.
## 'magick' package is suggested to install to give better rasterization.
## 
## Set `ht_opt$message = FALSE` to turn off this message.

4.9.2 Violin plots

Dividing by sample

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hPGCLCs", fill = "Genomic region")
## Warning: Removed 57 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning: Removed 57 rows containing non-finite values (`stat_summary()`).

Dividing by genomic region

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hPGCLCs", fill = "Sample", colors = sample_colors)
## Warning: Removed 57 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
## Warning: Removed 57 rows containing non-finite values (`stat_summary()`).

4.3.3 Top 20 DMR-associated genes

Useful to visualize single DMR with exon and CpG annotations. Here the top 20 significant DMRs are plotted. Rows in DMRs are sorted by “areaStat”, which is the sum of test statistics of all CpG sites in the region. Remember that we do not have p-value for DMRs because in DMR calling, the statistical test is perform for each CpG, and then the significant CpG are merged into regions. It is very difficult to estimate FDR at the region-level, based on the site level test results (p-values or FDR). This is a difficult question for all types of genome-wide assays such as the ChIP-seq.

Here we can filter for only differentially methylated promoters of protein coding genes

DMRs <- DMRs[grep("Promoter", DMRs$ChipSeekerAnn),] %>% filter(gene_biotype == "protein_coding")
nr_of_top <- 20

if (nrow(DMRs) > 20){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs)}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[i,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 500, main = paste0(DMRs[i,]$hgnc_symbol, " - ", DMRs[i,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: IFITM1, HTRA4, PLAGL1, PPP1R13L, PLEKHG5, CTSF, MCF2L, TRIO, SLITRK2, TNRC18, COL9A2, C5orf63, COMT, FAM83G, ZSCAN30, ARHGEF4, SLC15A4, CDH24, CHCHD2, GNAS

4.3.3 Top 10 DMR-associated genes where there is gain in methylation

Here the top 10 DMRs where methylation is higher in hEGCLCs (and lower in hiPSCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy<0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy<0, ])}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy<0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy<0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy<0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: HTRA4, PLAGL1, PPP1R13L, PLEKHG5, CTSF, TRIO, TNRC18, COL9A2, C5orf63, COMT

4.3.4 Top 10 DMR-associated genes where there is loss in methylation

Here the top 10 DMRs where methylation is higher in hiPSCs (and lower in hEGCLCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy>0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy>0, ])}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy>0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy>0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy>0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: IFITM1, MCF2L, SLITRK2, FAM83G, SLC15A4, CHCHD2, LSR, TMEM121, AXIN2, PABPC4L

4.3.5 Exploration of PGC markers and other relevant genes

SelectedGenes <- c("PIWIL2", "NANOS3", "IGF1", "DNMT3L", "PIFO", "SOX15", "CSF3R", "NUMB", "CXCR4")
PresentGenes <- vector()
AbsentGenes <- vector()
for(gene_index in SelectedGenes){
  if (any(DMRs$hgnc_symbol %in% gene_index)) {
    PresentGenes <- c(PresentGenes, gene_index)
    plotDMRs(bsseq_obj, regions=DMRs[DMRs$hgnc_symbol %in% gene_index,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(gene_index," - ", DMRs[DMRs$hgnc_symbol %in% gene_index,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
  } else{ 
    AbsentGenes <- c(AbsentGenes, gene_index)
    next}
}

These genes were found: IGF1, DNMT3L, PIFO, SOX15, CSF3R, NUMB These genes were not found: PIWIL2, NANOS3, CXCR4

4.4 hEGCLCs_p10 vs hEGCLCs_p5

comparison <- "hEGCLCs_p10vshEGCLCs_p5"
DMRs <- DMRs_annotated4[[comparison]]
DMRs$DMR <- paste0("DMR", rownames(DMRs))
DMRs_gr <- makeGRangesFromDataFrame(DMRs, keep.extra.columns = TRUE)
m <- findOverlaps(rowRanges(bsseq_obj), DMRs_gr)

bsseq_obj_subset <- bsseq_obj[queryHits(m),]

rowData(bsseq_obj_subset)$DMR <- DMRs_gr[subjectHits(m), ]$DMR
rowData(bsseq_obj_subset)$GenomicAnn <- DMRs_gr[subjectHits(m), ]$ChipSeekerAnn
rowData(bsseq_obj_subset)$Gene <- DMRs_gr[subjectHits(m), ]$hgnc_symbol
rowData(bsseq_obj_subset)$GeneBiotype <- DMRs_gr[subjectHits(m),]$gene_biotype
rowData(bsseq_obj_subset)$GenomicRegion <- rowData(bsseq_obj_subset)$GenomicAnn

rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Exon")] <- "Exon"
rowData(bsseq_obj_subset)$GenomicRegion[grep(rowData(bsseq_obj_subset)$GenomicRegion, pattern = "Intron")] <- "Intron"

In total we have 20901 CpGs in the significant DMRs.

Calculate average methylation of CpGs in each DMR

assays(bsseq_obj_subset)$Meth <- getMeth(bsseq_obj_subset, type = "raw")

AvgMeth <- aggregate(assays(bsseq_obj_subset)$Meth, list(rowData(bsseq_obj_subset)$DMR), mean, na.rm = T)
rownames(AvgMeth) <- AvgMeth$Group.1
AvgMeth <- AvgMeth[, -1]

rowDataSE <- rowData(bsseq_obj_subset)
rownames(rowDataSE) <- rowData(bsseq_obj_subset)$DMR
rowDataSE <- rowDataSE[rownames(AvgMeth), ] 
SE <- SummarizedExperiment(assays=AvgMeth, colData=colData(bsseq_obj_subset), rowData = rowDataSE, )

assays(SE) <- assays(SE) %>% setNames("AvgMeth")

identical(rownames(assay(SE)), rownames(rowData(SE))) #TRUE
## [1] TRUE
## [1] TRUE

Setting colors

ScaledCols <- c('darkblue', "purple", "white", "lightgoldenrod1", 'goldenrod1')

metadata(SE)$anno_colors <- list(Type = cellTypes_colors,
                                               GenomicRegion = c ("3' UTR" = "#FFADAD", "5' UTR" = "#FFD6A5", "Distal Intergenic" = "#FDFFB6", "Downstream (<=300bp)" = "#FFFFFC", "Exon" = "#9BF6FF", "Intron" = "#A0C4FF", "Promoter (<=1kb)" = "#BDB2FF", "Promoter (1-2kb)" = "#FFC6FF", "Promoter (2-3kb)" = "#CAFFBF"))

sample_colors <- c("hiPSC_rep1" = "#dd1c77", "hiPSC_rep2" = "#dd1c77", "hiPSC_rep3" = "#dd1c77", "hiPSC_rep4" = "#dd1c77", "hiPSC_rep5" = "#dd1c77",
                   "iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
                   "hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
                   "hEGCLC_p5_rep1" = "#ff7f00", "hEGCLC_p5_rep2" = "#ff7f00", "hEGCLC_p5_rep3" = "#ff7f00",
                    "hEGCLC_p10_rep1" = "#ffc400", "hEGCLC_p10_rep2" = "#ffc400", "hEGCLC_p10_rep3" = "#ffc400")

4.4.1 Heatmaps

Without sample clustering

break_coord <- assay(SE) %>% quantile(seq(0,1, by=0.1), na.rm =TRUE)
ExpCols <- viridisLite::plasma(length(seq(head(break_coord,1), tail(break_coord,1), by=0.2)))

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hEGCLCs_p5", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0)

With sample clustering

sechm(SE, features = rownames(SE), assayName = "AvgMeth", hmcols=ExpCols, breaks=seq(head(break_coord,1), tail(break_coord,1), by=0.2), show_colnames=TRUE, do.scale=FALSE, cluster_cols = TRUE, column_title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hEGCLCs_p5", top_annotation = c("Type"), left_annotation = c("GenomicRegion"), gaps_row = "GenomicRegion", row_title_rot = 0)

4.4.2 Violin plots

Dividing by sample

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hEGCLCs_p5", fill = "Genomic region")
## Warning: Removed 14 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning: Removed 14 rows containing non-finite values (`stat_summary()`).

Dividing by genomic region

DMRviolin(bsseq_obj = SE, assayName = "AvgMeth", title = "Methylation Levels in DMRs (averaging CpGs) - hEGCLCs_p10 vs hEGCLCs_p5", fill = "Sample", colors = sample_colors)
## Warning: Removed 14 rows containing non-finite values (`stat_ydensity()`).
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
## Warning: Removed 14 rows containing non-finite values (`stat_summary()`).

4.4.3 Top 20 DMR-associated genes

Useful to visualize single DMR with exon and CpG annotations. Here the top 20 significant DMRs are plotted. Rows in DMRs are sorted by “areaStat”, which is the sum of test statistics of all CpG sites in the region. Remember that we do not have p-value for DMRs because in DMR calling, the statistical test is perform for each CpG, and then the significant CpG are merged into regions. It is very difficult to estimate FDR at the region-level, based on the site level test results (p-values or FDR). This is a difficult question for all types of genome-wide assays such as the ChIP-seq.

Here we can filter for only differentially methylated promoters of protein coding genes

DMRs <- DMRs[grep("Promoter", DMRs$ChipSeekerAnn),] %>% filter(gene_biotype == "protein_coding")
nr_of_top <- 20

if (nrow(DMRs) > 20){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs)}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[i,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 500, main = paste0(DMRs[i,]$hgnc_symbol, " - ", DMRs[i,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: IFITM1, SLITRK2, SPACA1, GCSAML, PLPP2, RBM46, RHBDF1, NUPR2, USP6, CCDC17, NKAPL, MAL, MARCKS, TLE7, NCL, DMGDH, FAM228A, GRID2IP, FZD5, SH2B3

4.4.3 Top 10 DMR-associated genes where there is gain in methylation

Here the top 10 DMRs where methylation is higher in iMeLCs (and lower in hiPSCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy<0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy<0, ])}

# for(i in 1:nr_plots){
# plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy<0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy<0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy<0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
# }

The genes are: ITGB2, GPR132, AGAP3, CASZ1, MCF2L, MEIS3, NA, NA, NA, NA

4.4.4 Top 10 DMR-associated genes where there is loss in methylation

Here the top 10 DMRs where methylation is higher in hiPSCs (and lower in iMeLCs).

nr_of_top <- 10

if (nrow(DMRs[DMRs$diff.Methy>0, ]) > 10){nr_plots <- nr_of_top}else{nr_plots <- nrow(DMRs[DMRs$diff.Methy>0, ])}

for(i in 1:nr_plots){
plotDMRs(bsseq_obj, regions=DMRs[DMRs$diff.Methy>0, ][i, ], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(DMRs[DMRs$diff.Methy>0, ][i, ]$hgnc_symbol, " - ", DMRs[DMRs$diff.Methy>0, ][i, ]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
}

The genes are: IFITM1, SLITRK2, SPACA1, GCSAML, PLPP2, RBM46, RHBDF1, NUPR2, USP6, CCDC17

4.4.5 Exploration of PGC markers and other relevant genes

SelectedGenes <- c("PIWIL2", "NANOS3", "IGF1", "DNMT3L", "PIFO", "SOX15", "CSF3R", "NUMB", "CXCR4")
PresentGenes <- vector()
AbsentGenes <- vector()
for(gene_index in SelectedGenes){
  if (any(DMRs$hgnc_symbol %in% gene_index)) {
    PresentGenes <- c(PresentGenes, gene_index)
    plotDMRs(bsseq_obj, regions=DMRs[DMRs$hgnc_symbol %in% gene_index,], testCovariate=3, annoTrack=annoTrack, qval = FALSE, stat = FALSE, extend = 5000, main = paste0(gene_index," - ", DMRs[DMRs$hgnc_symbol %in% gene_index,]$ChipSeekerAnn), verbose = FALSE, horizLegend = TRUE, col = c(rep(cellTypes_colors[1], each = 5), rep(cellTypes_colors[2], each = 4),rep(cellTypes_colors[3:5], each = 3)))
  } else{ 
    AbsentGenes <- c(AbsentGenes, gene_index)
    next}
}

These genes were found: PIWIL2, SOX15 These genes were not found: NANOS3, IGF1, DNMT3L, PIFO, CSF3R, NUMB, CXCR4

5. Date and Session Info

date()
## [1] "Tue Aug 19 12:54:02 2025"
sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] org.Hs.eg.db_3.16.0                     
##  [2] TxDb.Hsapiens.UCSC.hg38.knownGene_3.16.0
##  [3] GenomicFeatures_1.50.4                  
##  [4] AnnotationDbi_1.60.2                    
##  [5] sechm_1.6.0                             
##  [6] annotatr_1.24.0                         
##  [7] dmrseq_1.18.1                           
##  [8] bsseq_1.34.0                            
##  [9] SummarizedExperiment_1.28.0             
## [10] Biobase_2.58.0                          
## [11] MatrixGenerics_1.10.0                   
## [12] matrixStats_1.0.0                       
## [13] GenomicRanges_1.50.2                    
## [14] GenomeInfoDb_1.34.9                     
## [15] IRanges_2.32.0                          
## [16] S4Vectors_0.36.2                        
## [17] BiocGenerics_0.44.0                     
## [18] tidyr_1.3.0                             
## [19] dplyr_1.1.2                             
## [20] gridExtra_2.3                           
## [21] scales_1.2.1                            
## [22] plotly_4.10.2                           
## [23] ggplot2_3.4.2                           
## 
## loaded via a namespace (and not attached):
##   [1] circlize_0.4.15               AnnotationHub_3.6.0          
##   [3] BiocFileCache_2.6.1           plyr_1.8.8                   
##   [5] lazyeval_0.2.2                splines_4.2.1                
##   [7] BiocParallel_1.32.6           digest_0.6.33                
##   [9] ca_0.71.1                     foreach_1.5.2                
##  [11] htmltools_0.5.5               fansi_1.0.4                  
##  [13] magrittr_2.0.3                memoise_2.0.1                
##  [15] BSgenome_1.66.3               cluster_2.1.4                
##  [17] doParallel_1.0.17             tzdb_0.4.0                   
##  [19] limma_3.54.2                  ComplexHeatmap_2.14.0        
##  [21] Biostrings_2.66.0             readr_2.1.4                  
##  [23] vroom_1.6.3                   R.utils_2.12.2               
##  [25] prettyunits_1.1.1             colorspace_2.1-0             
##  [27] blob_1.2.4                    rappdirs_0.3.3               
##  [29] xfun_0.39                     crayon_1.5.2                 
##  [31] RCurl_1.98-1.12               jsonlite_1.8.7               
##  [33] iterators_1.0.14              glue_1.6.2                   
##  [35] registry_0.5-1                gtable_0.3.3                 
##  [37] zlibbioc_1.44.0               XVector_0.38.0               
##  [39] V8_4.3.0                      GetoptLong_1.0.5             
##  [41] DelayedArray_0.24.0           Rhdf5lib_1.20.0              
##  [43] shape_1.4.6                   HDF5Array_1.26.0             
##  [45] rngtools_1.5.2                DBI_1.1.3                    
##  [47] randomcoloR_1.1.0.1           Rcpp_1.0.11                  
##  [49] viridisLite_0.4.2             xtable_1.8-4                 
##  [51] progress_1.2.2                clue_0.3-64                  
##  [53] bumphunter_1.40.0             bit_4.0.5                    
##  [55] htmlwidgets_1.6.2             httr_1.4.6                   
##  [57] RColorBrewer_1.1-3            ellipsis_0.3.2               
##  [59] farver_2.1.1                  pkgconfig_2.0.3              
##  [61] XML_3.99-0.14                 R.methodsS3_1.8.2            
##  [63] sass_0.4.7                    dbplyr_2.3.3                 
##  [65] locfit_1.5-9.7                utf8_1.2.3                   
##  [67] labeling_0.4.2                tidyselect_1.2.0             
##  [69] rlang_1.1.1                   reshape2_1.4.4               
##  [71] later_1.3.1                   munsell_0.5.0                
##  [73] BiocVersion_3.16.0            tools_4.2.1                  
##  [75] cachem_1.0.8                  cli_3.6.1                    
##  [77] generics_0.1.3                RSQLite_2.3.1                
##  [79] evaluate_0.21                 stringr_1.5.0                
##  [81] fastmap_1.1.1                 yaml_2.3.7                   
##  [83] outliers_0.15                 knitr_1.43                   
##  [85] bit64_4.0.5                   purrr_1.0.1                  
##  [87] KEGGREST_1.38.0               doRNG_1.8.6                  
##  [89] nlme_3.1-162                  sparseMatrixStats_1.10.0     
##  [91] mime_0.12                     R.oo_1.25.0                  
##  [93] xml2_1.3.5                    biomaRt_2.54.1               
##  [95] compiler_4.2.1                rstudioapi_0.15.0            
##  [97] filelock_1.0.2                curl_5.0.1                   
##  [99] png_0.1-8                     interactiveDisplayBase_1.36.0
## [101] tibble_3.2.1                  bslib_0.5.0                  
## [103] stringi_1.7.12                highr_0.10                   
## [105] lattice_0.21-8                Matrix_1.6-0                 
## [107] permute_0.9-7                 vctrs_0.6.3                  
## [109] pillar_1.9.0                  lifecycle_1.0.3              
## [111] rhdf5filters_1.10.1           BiocManager_1.30.20          
## [113] GlobalOptions_0.1.2           jquerylib_0.1.4              
## [115] data.table_1.14.8             bitops_1.0-7                 
## [117] seriation_1.4.2               httpuv_1.6.11                
## [119] rtracklayer_1.58.0            R6_2.5.1                     
## [121] BiocIO_1.8.0                  TSP_1.2-4                    
## [123] promises_1.2.0.1              codetools_0.2-19             
## [125] gtools_3.9.4                  rhdf5_2.42.1                 
## [127] rjson_0.2.21                  withr_2.5.0                  
## [129] regioneR_1.30.0               GenomicAlignments_1.34.1     
## [131] Rsamtools_2.14.0              GenomeInfoDbData_1.2.9       
## [133] parallel_4.2.1                hms_1.1.3                    
## [135] grid_4.2.1                    rmarkdown_2.23               
## [137] DelayedMatrixStats_1.20.0     Cairo_1.6-0                  
## [139] Rtsne_0.16                    shiny_1.7.4.1                
## [141] restfulr_0.0.15