Comparison escapees vs non-escapees
Gaja Matassa
Date: July 24, 2025
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(bsseq)
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## 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 objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:plyr':
##
## 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:plyr':
##
## desc
## 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
## The following object is masked from 'package:plyr':
##
## 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(dmrseq)
library(ggplot2)
source("../Methylation_helper.R")
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
##
## expandSet random seed and number of permutations
set.seed(123)
num_permutations <- params$num_permutations1. Data loading and fast exploration by PCA
Loading bsseq objects
bsseq_obj <- readRDS("~/EGCLCclones/2.PreliminaryAnalysis/bsseq_obj.rds")
bsseq_obj_escapees <- readRDS(paste0(params$InputFolder, "bsseq_obj_escapees.rds"))
bsseq_obj_NonEscapees <- readRDS(paste0(params$InputFolder, "bsseq_obj_NonEscapees.rds"))- Our bsseq object contain 3734728 CpG loci.
- Among them 8941 CpG loci belong to escapees regions (escapee CpGs).
- While 3635237 CpG loci do not belong to escapees regions (non-escapee CpGs).
Colors for figures
cellTypes_colors <- c(hiPSCs ="#dd1c77", iMeLCs ="#377eb8", hPGCLCs ="#4daf4a", hEGCLCs_p5 = "#ff7f00", hEGCLCs_p10 = "#ffc400")PCA only on escapee CpGs
pca_res <- do_PCA(getMeth(bsseq_obj_escapees, type = "raw"))
sample_anno <- bsseq::pData(bsseq_obj_escapees)
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = NULL, custom_colors = cellTypes_colors, point_size = 5, show_labels = TRUE)
PCA only on non-escapee CpGs
pca_res <- do_PCA(getMeth(bsseq_obj_NonEscapees, type = "raw"))
sample_anno <- bsseq::pData(bsseq_obj_NonEscapees)
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = NULL, custom_colors = cellTypes_colors, point_size = 5, show_labels = TRUE)
As expected it resembles the PCA done on the complete one, since we removed only 8941 out of a total of 3734728 CpG loci.
2.Generation of 500 random sets of 8941 non-escapee CpGs and distributions plots
Let’s generate random indices to randomly select 8941 CpG loci among the non-escapee ones.
random_sets <- list()
for (i in 1:num_permutations) {
random_indices <- sample(1:dim(bsseq_obj_NonEscapees)[1], dim(bsseq_obj_escapees)[1], replace = FALSE)
bsseq_obj_randomNonEscapees <- bsseq_obj_NonEscapees[random_indices, ]
random_sets[[i]] <- bsseq_obj_randomNonEscapees
}Head of BSSEQ with escapee CpGs
getMeth(bsseq_obj_escapees, type = "raw") %>% head()## hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 hiPSC_rep5 iMeLC_rep1
## [1,] 0.9545455 0.9473684 0.9482759 0.9821429 0.9473684 0.9767442
## [2,] 1.0000000 0.9803922 0.9615385 0.9800000 0.9797980 0.9714286
## [3,] 0.8684211 0.9298246 0.7313433 0.8985507 0.8620690 0.9090909
## [4,] 0.8888889 0.9166667 0.8666667 0.9090909 0.8875000 0.9600000
## [5,] 0.7500000 0.8431373 0.8163265 0.7083333 0.8260870 0.8372093
## [6,] 0.8000000 0.8095238 0.8260870 0.8260870 0.8333333 0.8750000
## iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,] 0.9813084 0.8909091 0.9534884 0.9062500 1.0000000 0.9428571
## [2,] 0.9622642 0.9230769 0.9032258 0.9600000 0.9545455 0.9677419
## [3,] 0.7758621 0.7608696 0.8666667 0.8108108 0.6969697 0.7948718
## [4,] 0.8198198 0.8297872 0.8666667 0.8437500 0.7741935 0.8837209
## [5,] 0.8461538 0.7894737 0.7333333 0.7307692 0.7272727 0.8378378
## [6,] 0.6551724 0.9000000 0.8500000 0.8750000 0.8947368 0.6944444
## hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,] 0.9230769 0.9081633 0.9318182 0.9425287
## [2,] 0.9850746 0.9484536 0.9594595 0.9827586
## [3,] 0.8426966 0.8750000 0.8400000 0.8939394
## [4,] 0.8571429 0.8450704 0.8311688 0.8928571
## [5,] 0.7333333 0.6949153 0.8730159 0.8896552
## [6,] 0.7812500 0.7692308 0.8888889 0.8358209
## hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,] 0.9607843 0.9487179
## [2,] 0.9775281 1.0000000
## [3,] 0.9021739 0.8895349
## [4,] 0.9210526 0.8742857
## [5,] 0.8226950 0.8500000
## [6,] 0.8367347 0.9333333Escapee CpG methylation distribution plot
plotEmpiricalDistribution(bsseq_obj_escapees,
bySample = FALSE,
testCovariate = "Type",
adj = 3) + guides(linetype="none") +
labs(title = "Escapee CpG methylation distribution by cell type") + ggplot2::scale_color_manual(values = cellTypes_colors)
plotEmpiricalDistribution(bsseq_obj_escapees,
bySample = FALSE,
testCovariate = "Line",
adj = 3) + guides(linetype="none") +
labs(title = "Escapee CpG methylation distribution collapsing cell types")
Non-Escapee CpG methylation distribution plot
plotEmpiricalDistribution(bsseq_obj_NonEscapees,
bySample = FALSE,
testCovariate = "Type",
adj = 3) + guides(linetype="none") +
labs(title = "Non-escapee CpG methylation distribution by cell type") + ggplot2::scale_color_manual(values = cellTypes_colors)
plotEmpiricalDistribution(bsseq_obj_NonEscapees,
bySample = FALSE,
testCovariate = "Line",
adj = 3) + guides(linetype="none") +
labs(title = "Non-escapee CpG methylation distribution collapsing cell types")
Head of first BSSEQ with randomly selected non-escapee CpGs
getMeth(random_sets[[1]], type = "raw") %>% head()## hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 hiPSC_rep5 iMeLC_rep1
## [1,] 0.01694915 0.03125000 0.02272727 0.01162791 0.031413613 0.01785714
## [2,] 0.05405405 0.02173913 0.02439024 0.00000000 0.022875817 0.04255319
## [3,] 0.00000000 0.00000000 0.00000000 0.00000000 0.002808989 0.00000000
## [4,] 1.00000000 1.00000000 1.00000000 1.00000000 1.000000000 1.00000000
## [5,] 0.52777778 0.48333333 0.50000000 0.61111111 0.704545455 0.53488372
## [6,] 0.05555556 0.00000000 0.01923077 0.00000000 0.006172840 0.00000000
## iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,] 0.007462687 0.01666667 0.04000000 0.02127660 0.0000000 0.05769231
## [2,] 0.022471910 0.00000000 0.01886792 0.02083333 0.0000000 0.00000000
## [3,] 0.000000000 0.00000000 0.00000000 0.00000000 0.0000000 0.01515152
## [4,] 1.000000000 0.95238095 0.95833333 0.90476190 1.0000000 1.00000000
## [5,] 0.510869565 0.63157895 0.57894737 0.42857143 0.5555556 0.59701493
## [6,] 0.000000000 0.00000000 0.00000000 0.00000000 0.0000000 0.03174603
## hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,] 0.01600000 0.06194690 0.03361345 0.047904192
## [2,] 0.03488372 0.01587302 0.01470588 0.032110092
## [3,] 0.00000000 0.00000000 0.00000000 0.000000000
## [4,] 0.92105263 0.94444444 0.96428571 1.000000000
## [5,] 0.36842105 0.48351648 0.41025641 0.700854701
## [6,] 0.01282051 0.00000000 0.00000000 0.004016064
## hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,] 0.02681992 0.028469751
## [2,] 0.01731602 0.009009009
## [3,] 0.03103448 0.004048583
## [4,] 0.94444444 1.000000000
## [5,] 0.81779661 0.693023256
## [6,] 0.04291845 0.009302326First set of randomly selected non-escapee CpG methylation distribution plot
plotEmpiricalDistribution(random_sets[[1]],
bySample = FALSE,
testCovariate = "Type",
adj = 3) + guides(linetype="none") +
labs(title = "Non-escapee CpG methylation distribution by cell type (1)") + ggplot2::scale_color_manual(values = cellTypes_colors)
Head of second BSSEQ with randomly selected non-escapee CpGs
getMeth(random_sets[[2]], type = "raw") %>% head()## hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 hiPSC_rep5 iMeLC_rep1
## [1,] 0.9245283 0.9831933 0.8739496 0.9490446 0.95702479 0.96581197
## [2,] 0.0000000 0.0000000 0.0000000 0.0000000 0.00000000 0.02272727
## [3,] 0.9600000 0.9032258 0.8928571 0.9142857 0.96078431 0.96296296
## [4,] 0.0000000 0.0000000 0.0000000 0.0000000 0.01285347 0.00000000
## [5,] 0.9333333 0.9189189 0.9148936 0.8688525 0.93072289 0.98333333
## [6,] 0.9503546 0.9669604 0.9559748 0.9633333 0.94202899 0.96467391
## iMeLC_rep2 iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3
## [1,] 0.9322917 0.98290598 0.9558824 0.93333333 0.9142857 0.89333333
## [2,] 0.0000000 0.00000000 0.0000000 0.00000000 0.0000000 0.01923077
## [3,] 0.8367347 0.82142857 0.9047619 0.85185185 0.9090909 0.90322581
## [4,] 0.0000000 0.01111111 0.0000000 0.01694915 0.0000000 0.00000000
## [5,] 0.9107143 0.81034483 0.8709677 0.95652174 0.8965517 0.95348837
## [6,] 0.9367429 0.96167883 0.9443378 0.95693780 0.9655172 0.92073171
## hEGCLC_p5_rep1 hEGCLC_p5_rep2 hEGCLC_p5_rep3 hEGCLC_p10_rep1
## [1,] 0.89784946 0.89171975 0.916666667 0.943025540
## [2,] 0.01428571 0.01204819 0.037037037 0.004608295
## [3,] 0.86792453 0.82692308 0.879310345 0.934782609
## [4,] 0.01398601 0.01459854 0.008849558 0.000000000
## [5,] 0.85869565 0.84523810 0.840909091 0.948529412
## [6,] 0.91095890 0.90218978 0.909433962 0.950000000
## hEGCLC_p10_rep2 hEGCLC_p10_rep3
## [1,] 0.970183486 0.958236659
## [2,] 0.009950249 0.005555556
## [3,] 0.933333333 0.914285714
## [4,] 0.023529412 0.000000000
## [5,] 0.950226244 0.940476190
## [6,] 0.981132075 0.945454545Second set of randomly selected non-escapee CpG methylation distribution plot
plotEmpiricalDistribution(random_sets[[2]],
bySample = FALSE,
testCovariate = "Type",
adj = 3) + guides(linetype="none") +
labs(title = "Non-escapee CpG methylation distribution by cell type (2)") + ggplot2::scale_color_manual(values = cellTypes_colors)
3. Kolmogorov-Smirnov test
Example of KS test result comparing distributions in case of hEGCLC_p10_rep1 sample: escapee CpG methylation levels vs first randomly selected non-escapee set.
ks_test_result <- ks.test(getMeth(bsseq_obj_escapees, type = "raw")[, "hEGCLC_p10_rep1"], getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_p10_rep1"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(bsseq_obj_escapees, type = "raw")[, "hEGCLC_p10_rep1"] and getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_p10_rep1"]
## D = 0.50617, p-value < 2.2e-16
## alternative hypothesis: two-sidedExample of KS test result comparing distributions in case of hEGCLC_p10_rep1 sample: first randomly selected non-escapee set vs second randomly selected non-escapee set.
ks_test_result <- ks.test(getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_p10_rep1"], getMeth(random_sets[[2]], type = "raw")[, "hEGCLC_p10_rep1"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_p10_rep1"] and getMeth(random_sets[[2]], type = "raw")[, "hEGCLC_p10_rep1"]
## D = 0.010185, p-value = 0.7498
## alternative hypothesis: two-sidedExample of KS test result comparing distributions: escapee CpG methylation levels (hiPSC_rep5 sample) vs first randomly selected non-escapee set (hEGCLC_p10_rep3 sample).
ks_test_result <- ks.test(getMeth(random_sets[[1]], type = "raw")[, "hiPSC_rep5"], getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_p10_rep3"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(random_sets[[1]], type = "raw")[, "hiPSC_rep5"] and getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_p10_rep3"]
## D = 0.055911, p-value = 1.956e-12
## alternative hypothesis: two-sidedExample of KS test result comparing distributions: escapee CpG methylation levels vs first randomly selected non-escapee set considering the mean levels across cell types (rowwise).
ks_test_result <- ks.test(getMeth(bsseq_obj_escapees, type = "raw") %>% rowMeans(na.rm = TRUE), getMeth(random_sets[[1]], type = "raw") %>% rowMeans(na.rm = TRUE))
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(bsseq_obj_escapees, type = "raw") %>% rowMeans(na.rm = TRUE) and getMeth(random_sets[[1]], type = "raw") %>% rowMeans(na.rm = TRUE)
## D = 0.51214, p-value < 2.2e-16
## alternative hypothesis: two-sidedExample of KS test result comparing distributions: first randomly selected non-escapee set vs second randomly selected non-escapee set considering the mean levels across cell types (rowwise).
ks_test_result <- ks.test(getMeth(random_sets[[1]], type = "raw") %>% rowMeans(na.rm = TRUE), getMeth(random_sets[[2]], type = "raw") %>% rowMeans(na.rm = TRUE))
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(random_sets[[1]], type = "raw") %>% rowMeans(na.rm = TRUE) and getMeth(random_sets[[2]], type = "raw") %>% rowMeans(na.rm = TRUE)
## D = 0.011073, p-value = 0.6435
## alternative hypothesis: two-sidedLet’s calculate the mean methylation levels across cell types in all sets of randomly selected non-escapee CpGs.
random_MeanMethLevels <- lapply(random_sets, function(bsseq){
getMeth(bsseq, type = "raw") %>% rowMeans(na.rm = TRUE)})Let’s do KS test comparing escapee CpG methylation levels vs all 500 randomly selected non-escapee sets.
ks_test_results <- lapply(random_MeanMethLevels, function(random_MeanMeth) {
ks.test(getMeth(bsseq_obj_escapees, type = "raw") %>% rowMeans(na.rm = TRUE), random_MeanMeth)
})
pvalues_kstest <- sapply(ks_test_results, function(result) {
result$p.value
})
adjusted_pvalues_kstest <- p.adjust(pvalues_kstest, method = "fdr")
print(adjusted_pvalues_kstest)## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [223] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [260] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [297] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [334] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [371] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [408] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [445] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [482] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0How many adjusted p-values are < 0.05?
table(adjusted_pvalues_kstest < 0.05)##
## TRUE
## 500Let’s do KS test comparing first randomly selected non-escapee set vs all 500 randomly selected non-escapee sets.
ks_test_results <- lapply(random_MeanMethLevels, function(random_MeanMeth) {
ks.test(random_MeanMethLevels[[1]], random_MeanMeth)
})
pvalues_kstest <- sapply(ks_test_results, function(result) {
result$p.value
})
adjusted_pvalues_kstest <- p.adjust(pvalues_kstest, method = "fdr")
print(adjusted_pvalues_kstest)## [1] 1.0000000 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [8] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [15] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [22] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [29] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [36] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [43] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [50] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [57] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [64] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [71] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [78] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [85] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [92] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [99] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [106] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [113] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [120] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [127] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [134] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [141] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [148] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [155] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [162] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [169] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [176] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [183] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [190] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [197] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [204] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [211] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [218] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [225] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [232] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [239] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [246] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [253] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [260] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [267] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [274] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [281] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [288] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [295] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [302] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [309] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [316] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [323] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [330] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [337] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [344] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [351] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [358] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [365] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [372] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [379] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [386] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [393] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [400] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [407] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [414] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [421] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [428] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [435] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [442] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [449] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [456] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [463] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [470] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [477] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [484] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232
## [491] 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987232 0.9987821
## [498] 0.9987232 0.9987232 0.9987232How many adjusted p-values are < 0.05?
table(adjusted_pvalues_kstest < 0.05)##
## FALSE
## 5004.Escapee and non escapee CpG methylation distribution plot
MeanMethLevels_df <- c(unlist(random_MeanMethLevels), getMeth(bsseq_obj_escapees, type = "raw") %>% rowMeans(na.rm = TRUE) %>% as.vector()) %>% as.data.frame()
colnames(MeanMethLevels_df) <- "MeanMeth"
MeanMethLevels_df$Distribution <- c(paste0("Random_", rep(1:num_permutations, each = dim(bsseq_obj_escapees)[1])), rep("Escapees_CpGs", each = dim(bsseq_obj_escapees)[1]))
MeanMethLevels_df$Category <- "Non Escapees"
MeanMethLevels_df$Category[(nrow(MeanMethLevels_df)+1-nrow(bsseq_obj_escapees)):nrow(MeanMethLevels_df)] <- "Escapees"ggplot(MeanMethLevels_df, aes(x = MeanMeth, color = Category, group = Distribution)) +
geom_density(alpha = 0.5) +
labs(
title = "Distribution CpG mean methylation levels in escapees and non-escapees regions",
x = "Mean Methylation Levels across cell types",
y = "Density",
color = "Category"
) +
scale_color_manual(values = c("Escapees" = "#FF7D00", "Non Escapees" = "#15616D")) +
theme_minimal()
5.Session Info
date()
## [1] "Thu Jul 24 16:02:20 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] tidyr_1.3.0 ggplot2_3.4.2
## [3] dmrseq_1.18.1 bsseq_1.34.0
## [5] SummarizedExperiment_1.28.0 Biobase_2.58.0
## [7] MatrixGenerics_1.10.0 matrixStats_1.0.0
## [9] GenomicRanges_1.50.2 GenomeInfoDb_1.34.9
## [11] IRanges_2.32.0 S4Vectors_0.36.2
## [13] BiocGenerics_0.44.0 dplyr_1.1.2
## [15] plyr_1.8.8
##
## loaded via a namespace (and not attached):
## [1] AnnotationHub_3.6.0 BiocFileCache_2.6.1
## [3] splines_4.2.1 BiocParallel_1.32.6
## [5] digest_0.6.33 foreach_1.5.2
## [7] htmltools_0.5.5 fansi_1.0.4
## [9] magrittr_2.0.3 memoise_2.0.1
## [11] BSgenome_1.66.3 tzdb_0.4.0
## [13] limma_3.54.2 Biostrings_2.66.0
## [15] readr_2.1.4 R.utils_2.12.2
## [17] prettyunits_1.1.1 colorspace_2.1-0
## [19] blob_1.2.4 rappdirs_0.3.3
## [21] ggrepel_0.9.3 xfun_0.39
## [23] crayon_1.5.2 RCurl_1.98-1.12
## [25] jsonlite_1.8.7 annotatr_1.24.0
## [27] iterators_1.0.14 glue_1.6.2
## [29] gtable_0.3.3 zlibbioc_1.44.0
## [31] XVector_0.38.0 DelayedArray_0.24.0
## [33] Rhdf5lib_1.20.0 HDF5Array_1.26.0
## [35] scales_1.2.1 DBI_1.1.3
## [37] rngtools_1.5.2 Rcpp_1.0.11
## [39] xtable_1.8-4 progress_1.2.2
## [41] bumphunter_1.40.0 bit_4.0.5
## [43] httr_1.4.6 RColorBrewer_1.1-3
## [45] ellipsis_0.3.2 pkgconfig_2.0.3
## [47] XML_3.99-0.14 R.methodsS3_1.8.2
## [49] farver_2.1.1 sass_0.4.7
## [51] dbplyr_2.3.3 locfit_1.5-9.7
## [53] utf8_1.2.3 tidyselect_1.2.0
## [55] labeling_0.4.2 rlang_1.1.1
## [57] reshape2_1.4.4 later_1.3.1
## [59] AnnotationDbi_1.60.2 munsell_0.5.0
## [61] BiocVersion_3.16.0 tools_4.2.1
## [63] cachem_1.0.8 cli_3.6.1
## [65] generics_0.1.3 RSQLite_2.3.1
## [67] evaluate_0.21 stringr_1.5.0
## [69] fastmap_1.1.1 yaml_2.3.7
## [71] outliers_0.15 knitr_1.43
## [73] bit64_4.0.5 purrr_1.0.1
## [75] KEGGREST_1.38.0 nlme_3.1-162
## [77] doRNG_1.8.6 sparseMatrixStats_1.10.0
## [79] mime_0.12 R.oo_1.25.0
## [81] xml2_1.3.5 biomaRt_2.54.1
## [83] compiler_4.2.1 rstudioapi_0.15.0
## [85] filelock_1.0.2 curl_5.0.1
## [87] png_0.1-8 interactiveDisplayBase_1.36.0
## [89] tibble_3.2.1 bslib_0.5.0
## [91] stringi_1.7.12 highr_0.10
## [93] GenomicFeatures_1.50.4 lattice_0.21-8
## [95] Matrix_1.6-0 permute_0.9-7
## [97] vctrs_0.6.3 pillar_1.9.0
## [99] lifecycle_1.0.3 rhdf5filters_1.10.1
## [101] BiocManager_1.30.20 jquerylib_0.1.4
## [103] data.table_1.14.8 bitops_1.0-7
## [105] httpuv_1.6.11 rtracklayer_1.58.0
## [107] R6_2.5.1 BiocIO_1.8.0
## [109] promises_1.2.0.1 codetools_0.2-19
## [111] gtools_3.9.4 rhdf5_2.42.1
## [113] rjson_0.2.21 withr_2.5.0
## [115] regioneR_1.30.0 GenomicAlignments_1.34.1
## [117] Rsamtools_2.14.0 GenomeInfoDbData_1.2.9
## [119] parallel_4.2.1 hms_1.1.3
## [121] grid_4.2.1 rmarkdown_2.23
## [123] DelayedMatrixStats_1.20.0 shiny_1.7.4.1
## [125] restfulr_0.0.15