Vignettes:

N.B. delta=0.25, p.threshold=0.0001

1. Environment setting

library(DSS)
## Loading required package: Biobase
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## 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
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: BiocParallel
## Loading required package: bsseq
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## 
## 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
## The following object is masked from 'package:Biobase':
## 
##     rowMedians
## Loading required package: parallel
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:IRanges':
## 
##     slice
## The following object is masked from 'package:S4Vectors':
## 
##     rename
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(ggplot2)
library(scales)

source("~/Transgenerational/EMSeq/NewData_fromLaSitua/Methylation_helper.R")
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
## 
##     expand
InputFolder <- params$InputFolder
OutputFolder <- params$OutputFolder

2. Loading of BSSeq object

Loading of bsseq object with CpGs shared by at least 75% of 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 
bsseq_obj
## An object of type 'BSseq' with
##   3644178 methylation loci
##   18 samples
## has not been smoothed
## All assays are in-memory

PCA plot

This PCA is done on the new filtered BSseq object (different from the one in the preliminary exploration script!)

cellTypes_colors <- c(hiPSCs ="#dd1c77", iMeLCs ="#377eb8", hPGCLCs ="#4daf4a", hEGCLCs_p5 = "#ff7f00", hEGCLCs_p10 = "#ffc400")
pca_res <- do_PCA(getMeth(BSseq = bsseq_obj, type = "raw"))
plot_PCA(pca_res = pca_res, anno = pData(bsseq_obj), col_anno = "Type", shape_anno = NULL, custom_colors = cellTypes_colors, point_size = 5)

plot_PCA(pca_res = pca_res, anno = pData(bsseq_obj), col_anno = "Type", shape_anno = NULL, pc_x = "PC2", pc_y = "PC3", custom_colors = cellTypes_colors, point_size = 5)

3. Differential Methylation analysis

Loading results obtained with DSS::DMLtest function

dml_tests_files <-  list.files(path = InputFolder, pattern = '*\\.rds', full.names = TRUE) 
#dml_tests_files <- dml_tests_files[-1]
dml_tests_files
##  [1] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hEGCLCs_p10vshEGCLCs_p5.rds"
##  [2] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hEGCLCs_p10vshiPSCs.rds"    
##  [3] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hEGCLCs_p10vshPGCLCs.rds"   
##  [4] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hEGCLCs_p10vsiMeLCs.rds"    
##  [5] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hEGCLCs_p5vshiPSCs.rds"     
##  [6] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hEGCLCs_p5vshPGCLCs.rds"    
##  [7] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hiPSCsvshPGCLCs.rds"        
##  [8] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/hiPSCsvsiMeLCs.rds"         
##  [9] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/iMeLCsvshEGCLCs_p5.rds"     
## [10] "/group/testa/Project/PGCLC/Rebuttal_paperEG2025/Output_DM/iMeLCsvshPGCLCs.rds"

Storing DMRs and DMLs Plotting number of DMRs

nr_tests <- length(dml_tests_files)

DM_tests_Res <- list()
DMR_Res <- list()
DML_Res <- list()

for(i in 1:nr_tests){
  
  DM_tests_Res[[i]] <- readRDS(dml_tests_files[i])
  DMR_Res[[i]] <- list()
  DML_Res[[i]] <- list()
  
  names(DM_tests_Res)[i] <- sub(".rds", "", basename(dml_tests_files))[i]
  names(DMR_Res)[i] <- sub(".rds", "", basename(dml_tests_files))[i]
  names(DML_Res)[i] <- sub(".rds", "", basename(dml_tests_files))[i]
  
  DMR_Res[[i]] <- callDMR(DM_tests_Res[[i]], delta=0.25, p.threshold=0.0001)
  DML_Res[[i]] <- callDML(DM_tests_Res[[i]], delta=0.25, p.threshold=0.0001)
  
  }
DMRbarplot(DMR=DMR_Res, comparison=1, y_lim = c(0,3300), interactive=TRUE)
DMRbarplot_length(DMR=DMR_Res, comparison=1)
DMRbarplot(DMR=DMR_Res, comparison=2, y_lim = c(0,3300), interactive=TRUE)
DMRbarplot_length(DMR=DMR_Res, comparison=2)
DMRbarplot(DMR=DMR_Res, comparison=3, y_lim = c(0,3300), interactive=TRUE)
DMRbarplot_length(DMR=DMR_Res, comparison=3)
DMRbarplot(DMR=DMR_Res, comparison=4, y_lim = c(0,3300), interactive=TRUE)
DMRbarplot_length(DMR=DMR_Res, comparison=4)

Barplot with all together

comparisons <- c("hEGCLCs_p10vshiPSCs", "hEGCLCs_p10vsiMeLCs", "hEGCLCs_p10vshPGCLCs", "hEGCLCs_p10vshEGCLCs_p5") #ordered with criteria

new_ylabel <- c("hEGCLCs_p10vshiPSCs" = "hiPSCs vs hEGCLCs_p10", 
                "hEGCLCs_p10vsiMeLCs" = "iMeLCs vs hEGCLCs_p10",
                "hEGCLCs_p10vshPGCLCs" = "hPGCLCs vs hEGCLCs_p10",
                "hEGCLCs_p10vshEGCLCs_p5" = "hEGCLCs_p5 vs hEGCLCs_p10")

DMRbarplot_all(DMR=DMR_Res, comparisons_ordered=comparisons, label_y=new_ylabel, y_lim = 3300, interactive=TRUE)
comparisons <- c("hEGCLCs_p10vshiPSCs", "hEGCLCs_p10vshPGCLCs", "hEGCLCs_p10vshEGCLCs_p5") #ordered with criteria

new_ylabel <- c("hEGCLCs_p10vshiPSCs" = "hiPSCs vs hEGCLCs_p10", 
                "hEGCLCs_p10vshPGCLCs" = "hPGCLCs vs hEGCLCs_p10",
                "hEGCLCs_p10vshEGCLCs_p5" = "hEGCLCs_p5 vs hEGCLCs_p10")

DMRbarplot_all(DMR=DMR_Res, comparisons_ordered=comparisons, label_y=new_ylabel, y_lim = 3300, interactive=TRUE)
comparisons <- c("hEGCLCs_p10vshiPSCs", "hEGCLCs_p10vshEGCLCs_p5") #ordered with criteria

new_ylabel <- c("hEGCLCs_p10vshiPSCs" = "hiPSCs vs hEGCLCs_p10", 
                "hEGCLCs_p10vshEGCLCs_p5" = "hEGCLCs_p5 vs hEGCLCs_p10")

DMRbarplot_all(DMR=DMR_Res, comparisons_ordered=comparisons, label_y=new_ylabel, y_lim = 1000, interactive=TRUE)

4. Interesting FAQ

N.B: Q: Why doesn’t callDMR function return FDR for identified DMRs?

A: 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. So I rather not to report a DMR-level FDR if it’s not accurate.

5. Saving DMRs and DMLs and Session Info

saveRDS(DM_tests_Res, paste0(OutputFolder, '/', 'DML_test_results.rds'))
saveRDS(DMR_Res, paste0(OutputFolder, '/', 'DMRs.rds'))
saveRDS(DML_Res, paste0(OutputFolder, '/', 'DMLs.rds'))
SessionInfo <- sessionInfo()
Date <- date()
Date
## [1] "Tue Aug 19 11:24:19 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] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] tidyr_1.3.0                 scales_1.2.1               
##  [3] plotly_4.10.2               ggplot2_3.4.2              
##  [5] DSS_2.46.0                  bsseq_1.34.0               
##  [7] SummarizedExperiment_1.28.0 MatrixGenerics_1.10.0      
##  [9] matrixStats_1.0.0           GenomicRanges_1.50.2       
## [11] GenomeInfoDb_1.34.9         IRanges_2.32.0             
## [13] S4Vectors_0.36.2            BiocParallel_1.32.6        
## [15] Biobase_2.58.0              BiocGenerics_0.44.0        
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.6                sass_0.4.7               
##  [3] viridisLite_0.4.2         jsonlite_1.8.7           
##  [5] splines_4.2.1             DelayedMatrixStats_1.20.0
##  [7] R.utils_2.12.2            gtools_3.9.4             
##  [9] bslib_0.5.0               highr_0.10               
## [11] BSgenome_1.66.3           GenomeInfoDbData_1.2.9   
## [13] Rsamtools_2.14.0          yaml_2.3.7               
## [15] pillar_1.9.0              lattice_0.21-8           
## [17] glue_1.6.2                limma_3.54.2             
## [19] digest_0.6.33             RColorBrewer_1.1-3       
## [21] XVector_0.38.0            colorspace_2.1-0         
## [23] htmltools_0.5.5           Matrix_1.6-0             
## [25] R.oo_1.25.0               XML_3.99-0.14            
## [27] pkgconfig_2.0.3           zlibbioc_1.44.0          
## [29] purrr_1.0.1               HDF5Array_1.26.0         
## [31] tibble_3.2.1              farver_2.1.1             
## [33] generics_0.1.3            ellipsis_0.3.2           
## [35] withr_2.5.0               cachem_1.0.8             
## [37] lazyeval_0.2.2            cli_3.6.1                
## [39] magrittr_2.0.3            crayon_1.5.2             
## [41] evaluate_0.21             R.methodsS3_1.8.2        
## [43] fansi_1.0.4               tools_4.2.1              
## [45] data.table_1.14.8         BiocIO_1.8.0             
## [47] lifecycle_1.0.3           Rhdf5lib_1.20.0          
## [49] munsell_0.5.0             locfit_1.5-9.7           
## [51] DelayedArray_0.24.0       Biostrings_2.66.0        
## [53] compiler_4.2.1            jquerylib_0.1.4          
## [55] rlang_1.1.1               rhdf5_2.42.1             
## [57] grid_4.2.1                RCurl_1.98-1.12          
## [59] rhdf5filters_1.10.1       rstudioapi_0.15.0        
## [61] htmlwidgets_1.6.2         rjson_0.2.21             
## [63] crosstalk_1.2.0           labeling_0.4.2           
## [65] bitops_1.0-7              rmarkdown_2.23           
## [67] restfulr_0.0.15           gtable_0.3.3             
## [69] codetools_0.2-19          R6_2.5.1                 
## [71] GenomicAlignments_1.34.1  knitr_1.43               
## [73] dplyr_1.1.2               rtracklayer_1.58.0       
## [75] fastmap_1.1.1             utf8_1.2.3               
## [77] permute_0.9-7             Rcpp_1.0.11              
## [79] vctrs_0.6.3               tidyselect_1.2.0         
## [81] xfun_0.39                 sparseMatrixStats_1.10.0