Comparison escapees vs non-escapees - Hsu et al.
Gaja Matassa
Date: July 09, 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_escapees <- readRDS(paste0(params$InputFolder, "bsseq_obj_escapees_Hsu.rds"))
bsseq_obj_NonEscapees <- readRDS(paste0(params$InputFolder, "bsseq_obj_NonEscapees_Hsu.rds"))- 302502 CpG loci belong to escapees regions (escapee CpGs).
- 27513050 CpG loci do not belong to escapees regions (non-escapee CpGs).
Colors for figures
cellTypes_colors <- c(hESC ="#dd1c77", D4hPGCLC ="#4daf4a")PCA only on escapee CpGs - For now focus on WGBS
pca_res <- do_PCA(assays(bsseq_obj_escapees)$Meth_WGBS)
sample_anno <- bsseq::pData(bsseq_obj_escapees)
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = "Line", custom_colors = cellTypes_colors, point_size = 5)
PCA only on non-escapee CpGs
pca_res <- do_PCA(assays(bsseq_obj_NonEscapees)$Meth_WGBS)
sample_anno <- bsseq::pData(bsseq_obj_NonEscapees)
plot_PCA(pca_res = pca_res, anno = sample_anno, col_anno = "Type", shape_anno = "Line", custom_colors = cellTypes_colors, point_size = 5)
2.Generation of 500 random sets of 302502 non-escapee CpGs and distributions plots
Let’s generate random indices to randomly select 302502 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
assays(bsseq_obj_escapees)$Meth_WGBS %>% head()## UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## [1,] 0.8750000 0.80 1.0000000 1.0000000
## [2,] 1.0000000 1.00 1.0000000 1.0000000
## [3,] 1.0000000 0.40 1.0000000 1.0000000
## [4,] 0.9090909 1.00 1.0000000 0.8333333
## [5,] 1.0000000 0.75 1.0000000 0.7777778
## [6,] 1.0000000 1.00 0.9166667 0.9166667
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## [1,] 0.8571429 1.0000000 1.0
## [2,] 1.0000000 1.0000000 1.0
## [3,] 0.8888889 0.7272727 1.0
## [4,] 0.7777778 0.8000000 1.0
## [5,] 0.8125000 1.0000000 0.8
## [6,] 0.8125000 0.8125000 1.0
## UCLA2_D4hPGCLC_rep3
## [1,] 0.8888889
## [2,] 0.8888889
## [3,] 1.0000000
## [4,] 0.8888889
## [5,] 1.0000000
## [6,] 1.0000000Escapee CpG methylation distribution plot
hist(assays(bsseq_obj_escapees)$Meth_WGB, breaks = 100, main = "Distribution of Methylation (Meth_WGBS)",
xlab = "Methylation Level", col = "skyblue", border = "white")
assays(bsseq_obj_escapees)$M <- assays(bsseq_obj_escapees)$M_WGBS
assays(bsseq_obj_escapees)$Cov <- assays(bsseq_obj_escapees)$Cov_WGBS
plotEmpiricalDistribution(bsseq_obj_escapees,
bySample = TRUE,
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
assays(bsseq_obj_NonEscapees)$M <- assays(bsseq_obj_NonEscapees)$M_WGBS
assays(bsseq_obj_NonEscapees)$Cov <- assays(bsseq_obj_NonEscapees)$Cov_WGBS
plotEmpiricalDistribution(bsseq_obj_NonEscapees,
bySample = TRUE,
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
assays(random_sets[[1]])$Meth_WGBS %>% head()## UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## [1,] 1.0000000 0.6000000 0.8000000 0.9333333
## [2,] 0.8947368 0.7500000 0.8181818 0.9166667
## [3,] 1.0000000 0.8750000 1.0000000 0.8571429
## [4,] 0.8235294 0.1666667 0.1428571 0.7142857
## [5,] 0.9000000 0.6000000 0.6666667 0.8461538
## [6,] 0.1000000 0.0000000 0.0000000 0.0000000
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## [1,] 1.0000000 0.9310345 0.9166667
## [2,] 0.6666667 0.8421053 1.0000000
## [3,] 0.9090909 0.9444444 1.0000000
## [4,] 0.6875000 0.2500000 0.4375000
## [5,] 0.8333333 0.7727273 0.6666667
## [6,] 0.0000000 0.0000000 0.0000000
## UCLA2_D4hPGCLC_rep3
## [1,] 0.8750000
## [2,] 1.0000000
## [3,] 0.9285714
## [4,] 0.2222222
## [5,] 0.8571429
## [6,] 0.0000000First set of randomly selected non-escapee CpG methylation distribution plot
assays(random_sets[[1]])$M <- assays(random_sets[[1]])$M_WGBS
assays(random_sets[[1]])$Cov <- assays(random_sets[[1]])$Cov_WGBS
plotEmpiricalDistribution(random_sets[[1]],
bySample = TRUE,
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
assays(random_sets[[2]])$Meth_WGBS %>% head()## UCLA1_hESC_rep3 UCLA1_D4hPGCLC_rep2 UCLA1_D4hPGCLC_rep3 UCLA2_hESC_rep2
## [1,] 0.8000000 0.7500000 0.6666667 1.0000000
## [2,] 0.6250000 0.6250000 0.5833333 1.0000000
## [3,] 0.8571429 0.8000000 0.8000000 0.9333333
## [4,] 0.9230769 0.5714286 0.7142857 0.6666667
## [5,] 0.8571429 1.0000000 0.8333333 0.8000000
## [6,] 0.5000000 1.0000000 0.7142857 1.0000000
## UCLA2_hESC_rep3 UCLA2_D4hPGCLC_rep1 UCLA2_D4hPGCLC_rep2
## [1,] 0.8750000 0.7142857 0.7500000
## [2,] 0.6666667 0.5714286 0.5555556
## [3,] 0.9166667 0.8800000 0.9230769
## [4,] 0.6250000 0.4545455 0.2500000
## [5,] 0.8000000 1.0000000 1.0000000
## [6,] 0.9687500 0.9000000 1.0000000
## UCLA2_D4hPGCLC_rep3
## [1,] 0.8000000
## [2,] 0.7500000
## [3,] 0.7857143
## [4,] 0.6153846
## [5,] 0.7500000
## [6,] 0.8571429Second set of randomly selected non-escapee CpG methylation distribution plot
assays(random_sets[[2]])$M <- assays(random_sets[[2]])$M_WGBS
assays(random_sets[[2]])$Cov <- assays(random_sets[[2]])$Cov_WGBS
plotEmpiricalDistribution(random_sets[[2]],
bySample = TRUE,
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 UCLA1_hESC_rep3 sample: escapee CpG methylation levels vs first randomly selected non-escapee set.
ks_test_result <- ks.test(assays(bsseq_obj_escapees)$Meth_WGBS[, "UCLA1_hESC_rep3"], assays(random_sets[[1]])$Meth_WGBS[, "UCLA1_hESC_rep3"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: assays(bsseq_obj_escapees)$Meth_WGBS[, "UCLA1_hESC_rep3"] and assays(random_sets[[1]])$Meth_WGBS[, "UCLA1_hESC_rep3"]
## D = 0.14213, p-value < 2.2e-16
## alternative hypothesis: two-sidedExample of KS test result comparing distributions in case of UCLA1_hESC_rep3 sample: first randomly selected non-escapee set vs second randomly selected non-escapee set.
ks_test_result <- ks.test(assays(random_sets[[1]])$Meth_WGBS[, "UCLA1_hESC_rep3"], assays(random_sets[[2]])$Meth_WGBS[, "UCLA1_hESC_rep3"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: assays(random_sets[[1]])$Meth_WGBS[, "UCLA1_hESC_rep3"] and assays(random_sets[[2]])$Meth_WGBS[, "UCLA1_hESC_rep3"]
## D = 0.0017454, p-value = 0.7461
## alternative hypothesis: two-sidedExample of KS test result comparing distributions: escapee CpG methylation levels (UCLA2_hESC_rep2 sample) vs first randomly selected non-escapee set (UCLA1_D4hPGCLC_rep3 sample).
ks_test_result <- ks.test(assays(bsseq_obj_escapees)$Meth_WGBS[, "UCLA2_hESC_rep2"], assays(random_sets[[1]])$Meth_WGBS[, "UCLA1_D4hPGCLC_rep3"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: assays(bsseq_obj_escapees)$Meth_WGBS[, "UCLA2_hESC_rep2"] and assays(random_sets[[1]])$Meth_WGBS[, "UCLA1_D4hPGCLC_rep3"]
## D = 0.3007, p-value < 2.2e-16
## 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(assays(bsseq_obj_escapees)$Meth_WGBS %>% rowMeans(na.rm = TRUE), assays(random_sets[[1]])$Meth_WGBS %>% rowMeans(na.rm = TRUE))
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: assays(bsseq_obj_escapees)$Meth_WGBS %>% rowMeans(na.rm = TRUE) and assays(random_sets[[1]])$Meth_WGBS %>% rowMeans(na.rm = TRUE)
## D = 0.2318, 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(assays(random_sets[[1]])$Meth_WGBS %>% rowMeans(na.rm = TRUE), assays(random_sets[[2]])$Meth_WGBS %>% rowMeans(na.rm = TRUE))
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: assays(random_sets[[1]])$Meth_WGBS %>% rowMeans(na.rm = TRUE) and assays(random_sets[[2]])$Meth_WGBS %>% rowMeans(na.rm = TRUE)
## D = 0.0024463, p-value = 0.3258
## 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){
assays(bsseq)$Meth_WGBS %>% 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(assays(bsseq_obj_escapees)$Meth_WGBS %>% 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.8387421 0.9800503 0.8210148 0.9715843 0.8210148 0.9940842
## [8] 0.9887370 0.8210148 0.8210148 0.8210148 0.9940842 0.8210148 0.9285689
## [15] 0.8210148 0.9800503 0.8210148 0.8210148 0.8210148 0.9940842 0.9940842
## [22] 0.9715843 0.8210148 0.8210148 0.9715843 0.9285689 0.8210148 0.9800503
## [29] 0.8210148 0.9800503 0.9800503 0.8210148 0.8210148 0.9887370 0.8654371
## [36] 0.8210148 0.9800503 0.8210148 0.9800503 0.9940842 0.8210148 0.8210148
## [43] 0.8210148 0.9887370 0.9285689 0.9940842 0.8210148 0.8210148 0.9285689
## [50] 0.9887370 0.9940842 0.9285689 0.9285689 0.9100922 0.9940842 0.9800503
## [57] 0.9800503 0.8210148 0.8210148 0.9531270 0.9889643 0.9889643 0.9336101
## [64] 0.9940842 0.9940842 0.9800503 0.8210148 0.9940842 0.8210148 0.9887370
## [71] 0.9940842 0.8210148 0.9940842 0.8210148 0.8210148 0.9940842 0.9940842
## [78] 0.8210148 0.8210148 0.8210148 0.8210148 0.9887370 0.9672040 0.9285689
## [85] 0.8210148 0.9940842 0.8210148 0.9285689 0.9800503 0.9800503 0.8210148
## [92] 0.8210148 0.9800503 0.9285689 0.9940842 0.8210148 0.9940842 0.9940842
## [99] 0.9800503 0.9940842 0.9889643 0.9940842 0.9887370 0.9940842 0.8210148
## [106] 0.8210148 0.9940842 0.8210148 0.9887370 0.9715843 0.8210148 0.9940842
## [113] 0.9715843 0.9436085 0.8210148 0.8210148 0.8210148 0.8210148 0.8210148
## [120] 0.9285689 0.9285689 0.8210148 0.8210148 0.9800503 0.8210148 0.9285689
## [127] 0.8210148 0.9940842 0.8574956 0.9436085 0.9285689 1.0000000 0.9940842
## [134] 0.8210148 0.9285689 0.9285689 0.9285689 0.8210148 0.8210148 0.8210148
## [141] 0.9940842 0.9940842 0.8210148 0.9285689 0.9889643 0.8210148 0.9940842
## [148] 0.8210148 0.9887370 0.8210148 0.8210148 0.9754522 0.9800503 0.9940842
## [155] 0.8210148 0.9285689 0.9800503 0.8210148 0.8210148 0.8210148 0.9436085
## [162] 0.9940842 0.8387421 0.8210148 0.9672040 0.9887370 0.9940842 0.8210148
## [169] 0.8210148 0.9285689 0.9940842 0.8210148 0.9800503 0.9889643 0.8210148
## [176] 0.8210148 0.9715843 1.0000000 0.9940842 0.9285689 0.8210148 0.9715843
## [183] 0.9889643 0.8210148 0.8387421 0.9887370 0.9976991 0.9285689 0.8210148
## [190] 0.8734122 0.9940842 0.9285689 0.9715843 0.8210148 0.9947971 0.9940842
## [197] 1.0000000 0.9940842 0.9441774 0.9800503 0.8210148 0.8210148 0.8210148
## [204] 0.9940842 0.8210148 0.8210148 0.9531270 0.8210148 0.8210148 0.9940842
## [211] 0.9940842 0.8210148 0.9436085 1.0000000 0.9887370 0.9940842 0.9715843
## [218] 0.9940842 0.9940842 0.9800503 0.8210148 0.8210148 0.8210148 0.8210148
## [225] 0.9285689 0.9800503 0.9583313 0.9285689 0.9887370 0.9800503 0.8210148
## [232] 0.8210148 0.8210148 0.9800503 0.9887370 0.9800503 0.8210148 0.8210148
## [239] 0.8210148 0.9940842 0.9800503 0.9800503 0.9800503 0.8210148 0.8210148
## [246] 0.9436085 0.9940842 0.8210148 0.8210148 0.9840050 0.9800503 0.8210148
## [253] 0.8362395 0.8210148 0.8210148 0.9757683 0.9940842 0.9800503 0.8210148
## [260] 0.8210148 0.9940842 0.9940842 0.8387421 0.9837700 0.9940842 0.9715843
## [267] 0.9800503 0.9887370 0.9887370 0.8210148 0.8210148 0.8210148 0.9887370
## [274] 0.9889643 0.8498827 0.8210148 0.8210148 0.9940842 0.9887370 0.9940842
## [281] 0.9940842 0.9940842 0.9887370 0.8210148 0.9887370 0.8465947 0.9800503
## [288] 0.9800503 0.9543261 0.8210148 0.9940842 0.8210148 0.9940842 0.9800503
## [295] 0.9285689 0.9285689 0.9940842 0.8210148 0.9887370 0.9800503 0.8210148
## [302] 0.8210148 0.8210148 0.8498827 0.9586344 0.9285689 0.9285689 0.9285689
## [309] 0.8210148 0.8210148 0.9436085 0.9940842 0.8210148 0.9889643 0.9940842
## [316] 0.9800503 0.9800503 0.8574956 0.9940842 0.8210148 0.8210148 0.9800503
## [323] 0.8210148 0.8498827 1.0000000 0.8210148 0.9285689 0.8210148 0.8210148
## [330] 0.9285689 0.8210148 0.8210148 0.8210148 0.9485111 0.8210148 0.8210148
## [337] 0.9800503 0.8210148 0.8210148 0.8210148 0.9482238 0.8210148 0.9940842
## [344] 0.9338633 0.9285689 0.9887370 0.8210148 0.8210148 0.8210148 0.9543261
## [351] 0.8210148 0.8210148 0.8210148 0.9940842 0.9285689 0.9889643 0.9940842
## [358] 1.0000000 0.9800503 0.9800503 0.8210148 0.9835307 0.9285689 0.9441774
## [365] 0.8210148 0.9800503 0.9292821 0.8362395 0.9715843 0.9940842 0.9835307
## [372] 0.8210148 0.8210148 0.9940842 0.9288007 0.9800503 0.9800503 0.9937170
## [379] 0.8210148 0.9887370 0.9887370 0.9940842 0.9800503 0.9800503 0.8210148
## [386] 0.9940842 0.8210148 0.8498827 0.9488023 0.9800503 0.9715843 0.8210148
## [393] 0.8210148 0.8210148 0.9800503 0.9940842 0.8287612 0.8210148 0.8210148
## [400] 0.9285689 0.8210148 0.9887370 0.9295314 0.8210148 0.8210148 0.8210148
## [407] 0.9887370 0.9940842 0.8210148 0.8210148 0.9285689 0.9285689 0.8387421
## [414] 1.0000000 0.9800503 0.8387421 0.9715843 0.9940842 0.9940842 0.9887370
## [421] 0.9292821 0.8210148 0.9715843 0.9285689 0.9940842 0.8210148 0.9436085
## [428] 0.9543261 0.8210148 0.9800503 0.9940842 0.8210148 0.8210148 0.8210148
## [435] 0.8773456 0.8210148 0.9764026 0.8210148 0.9285689 0.9800503 0.8210148
## [442] 0.8210148 0.9940842 0.9800503 0.9889643 0.8210148 0.9889643 0.8210148
## [449] 0.9889643 0.8210148 0.9940842 0.9887370 0.8210148 0.8210148 0.8210148
## [456] 0.9285689 0.8210148 0.9887370 0.9940842 0.9940842 1.0000000 0.8210148
## [463] 0.8210148 0.9940842 0.8210148 0.8574956 0.8210148 0.8210148 0.9285689
## [470] 0.8210148 0.9534217 0.9800503 0.8210148 0.9285689 0.9760851 0.8210148
## [477] 0.9940842 0.9285689 0.9887370 0.8210148 0.9672040 0.9800503 0.9436085
## [484] 0.9940842 0.9672040 0.9285689 0.9285689 0.8210148 0.9940842 0.9887370
## [491] 0.9940842 0.9800503 0.9940842 0.8210148 0.8210148 0.9285689 0.8210148
## [498] 0.8210148 0.8210148 0.8210148How many adjusted p-values are < 0.05?
table(adjusted_pvalues_kstest < 0.05)##
## FALSE
## 5004.Escapee and non escapee CpG methylation distribution plot
4.1 WGBS Meth
MeanMethLevels_df <- c(unlist(random_MeanMethLevels), assays(bsseq_obj_escapees)$Meth_WGBS %>% 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"Mean methylation levels collapsing cell types
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()
Mean methylation levels separating cell types
# Define function to compute rowMeans for specific sample types
get_rowMeans_by_type <- function(bsseq, type_label) {
sample_indices <- which(colData(bsseq)$Type == type_label)
rowMeans(assays(bsseq)$Meth_WGBS[, sample_indices, drop = FALSE], na.rm = TRUE)
}
# Apply the function to each random set for both hESC and D4hPGCLC
random_MeanMethLevels_hESC <- lapply(random_sets, get_rowMeans_by_type, type_label = "hESC")
random_MeanMethLevels_D4hPGCLC <- lapply(random_sets, get_rowMeans_by_type, type_label = "D4hPGCLC")
# Do the same for escapees
escapees_MeanMeth_hESC <- get_rowMeans_by_type(bsseq_obj_escapees, "hESC")
escapees_MeanMeth_D4hPGCLC <- get_rowMeans_by_type(bsseq_obj_escapees, "D4hPGCLC")
# Combine all into a single data frame
MeanMethLevels_df <- data.frame(
MeanMeth = c(
unlist(random_MeanMethLevels_hESC),
unlist(random_MeanMethLevels_D4hPGCLC),
escapees_MeanMeth_hESC,
escapees_MeanMeth_D4hPGCLC
),
Distribution = c(
paste0("Random_hESC_", rep(1:num_permutations, each = dim(bsseq_obj_escapees)[1])),
paste0("Random_D4hPGCLC_", rep(1:num_permutations, each = dim(bsseq_obj_escapees)[1])),
rep("Escapees_hESC", each = dim(bsseq_obj_escapees)[1]),
rep("Escapees_D4hPGCLC", each = dim(bsseq_obj_escapees)[1])
),
Category = c(
rep("Non Escapees", 2 * num_permutations * dim(bsseq_obj_escapees)[1]),
rep("Escapees", 2 * dim(bsseq_obj_escapees)[1])
),
Type = c(
rep("hESC", num_permutations * dim(bsseq_obj_escapees)[1]),
rep("D4hPGCLC", num_permutations * dim(bsseq_obj_escapees)[1]),
rep("hESC", dim(bsseq_obj_escapees)[1]),
rep("D4hPGCLC", dim(bsseq_obj_escapees)[1])
)
)
ggplot(MeanMethLevels_df, aes(x = MeanMeth, color = Category, group = Distribution)) +
geom_density(alpha = 0.5) +
labs(
title = "CpG Mean Methylation Levels by Category and Type",
x = "Mean Methylation Levels",
y = "Density",
color = "Category"
) +
facet_wrap(~ Type, scales = "free_y") +
scale_color_manual(values = c("Escapees" = "#FF7D00", "Non Escapees" = "#15616D")) +
theme_minimal()
# Create a new combined variable in the dataframe
MeanMethLevels_df$GroupLabel <- paste(MeanMethLevels_df$Category, MeanMethLevels_df$Type, sep = "_")
# Now plot using this combined label
ggplot(MeanMethLevels_df, aes(x = MeanMeth, color = GroupLabel, group = GroupLabel)) +
geom_density(alpha = 0.6, adjust = 1.5) + # smoother curves using `adjust`
labs(
title = "CpG Mean Methylation Levels by Category and Cell Type",
x = "Mean Methylation Levels",
y = "Density",
color = "Group"
) +
scale_color_manual(
values = c(
"Escapees_hESC" = "#FF7D00",
"Escapees_D4hPGCLC" = "#FFB347",
"Non Escapees_hESC" = "#15616D",
"Non Escapees_D4hPGCLC" = "#1E91A1"
)
) +
coord_cartesian(ylim = c(0, 6)) + # cut off at y = 6
theme_minimal()
4.2 bACEseq Meth
random_MeanHydroxyMethLevels <- lapply(random_sets, function(bsseq){
assays(bsseq)$Meth_hmC %>% rowMeans(na.rm = TRUE)})
MeanHydroxyMethLevels_df <- c(unlist(random_MeanHydroxyMethLevels), assays(bsseq_obj_escapees)$Meth_hmC %>% rowMeans(na.rm = TRUE) %>% as.vector()) %>% as.data.frame()
colnames(MeanHydroxyMethLevels_df) <- "MeanHydroxyMeth"
MeanHydroxyMethLevels_df$Distribution <- c(paste0("Random_", rep(1:num_permutations, each = dim(bsseq_obj_escapees)[1])), rep("Escapees_CpGs", each = dim(bsseq_obj_escapees)[1]))
MeanHydroxyMethLevels_df$Category <- "Non Escapees"
MeanHydroxyMethLevels_df$Category[(nrow(MeanHydroxyMethLevels_df)+1-nrow(bsseq_obj_escapees)):nrow(MeanHydroxyMethLevels_df)] <- "Escapees"ggplot(MeanHydroxyMethLevels_df, aes(x = MeanHydroxyMeth, color = Category, group = Distribution)) +
geom_density(alpha = 0.5) +
labs(
title = "Distribution CpG mean hydroxymethylation levels in escapees and non-escapees regions",
x = "Mean Hydroxymethylation Levels across cell types",
y = "Density",
color = "Category"
) +
scale_color_manual(values = c("Escapees" = "#FF7D00", "Non Escapees" = "#15616D")) +
coord_cartesian(xlim = c(0, 0.1)) +
theme_minimal()
# Define function to compute rowMeans for specific sample types
get_rowMeans_by_type_hydroxy <- function(bsseq, type_label) {
sample_indices <- which(colData(bsseq)$Type == type_label)
rowMeans(assays(bsseq)$Meth_hmC[, sample_indices, drop = FALSE], na.rm = TRUE)
}
# Apply the function to each random set for both hESC and D4hPGCLC
random_MeanHydroxyMethLevels_hESC <- lapply(random_sets, get_rowMeans_by_type_hydroxy, type_label = "hESC")
random_MeanHydroxyMethLevels_D4hPGCLC <- lapply(random_sets, get_rowMeans_by_type_hydroxy, type_label = "D4hPGCLC")
# Do the same for escapees
escapees_MeanHydroxyMeth_hESC <- get_rowMeans_by_type_hydroxy(bsseq_obj_escapees, "hESC")
escapees_MeanHydroxyMeth_D4hPGCLC <- get_rowMeans_by_type_hydroxy(bsseq_obj_escapees, "D4hPGCLC")
# Combine all into a single data frame
MeanHydroxyMethLevels_df <- data.frame(
MeanHydroxyMeth = c(
unlist(random_MeanHydroxyMethLevels_hESC),
unlist(random_MeanHydroxyMethLevels_D4hPGCLC),
escapees_MeanHydroxyMeth_hESC,
escapees_MeanHydroxyMeth_D4hPGCLC
),
Distribution = c(
paste0("Random_hESC_", rep(1:num_permutations, each = dim(bsseq_obj_escapees)[1])),
paste0("Random_D4hPGCLC_", rep(1:num_permutations, each = dim(bsseq_obj_escapees)[1])),
rep("Escapees_hESC", each = dim(bsseq_obj_escapees)[1]),
rep("Escapees_D4hPGCLC", each = dim(bsseq_obj_escapees)[1])
),
Category = c(
rep("Non Escapees", 2 * num_permutations * dim(bsseq_obj_escapees)[1]),
rep("Escapees", 2 * dim(bsseq_obj_escapees)[1])
),
Type = c(
rep("hESC", num_permutations * dim(bsseq_obj_escapees)[1]),
rep("D4hPGCLC", num_permutations * dim(bsseq_obj_escapees)[1]),
rep("hESC", dim(bsseq_obj_escapees)[1]),
rep("D4hPGCLC", dim(bsseq_obj_escapees)[1])
)
)
ggplot(MeanHydroxyMethLevels_df, aes(x = MeanHydroxyMeth, color = Category, group = Distribution)) +
geom_density(alpha = 0.5) +
labs(
title = "CpG Mean Hydroxy Methylation Levels by Category and Type",
x = "Mean Methylation Levels",
y = "Density",
color = "Category"
) +
facet_wrap(~ Type, scales = "free_y") +
scale_color_manual(values = c("Escapees" = "#FF7D00", "Non Escapees" = "#15616D")) +
theme_minimal()
# Create a new combined variable in the dataframe
MeanHydroxyMethLevels_df$GroupLabel <- paste(MeanHydroxyMethLevels_df$Category, MeanHydroxyMethLevels_df$Type, sep = "_")
# Now plot using this combined label
ggplot(MeanHydroxyMethLevels_df, aes(x = MeanHydroxyMeth, color = GroupLabel, group = GroupLabel)) +
geom_density(alpha = 0.6, adjust = 1.5) + # smoother curves using `adjust`
labs(
title = "Distribution CpG mean hydroxymethylation levels in escapees and non-escapees regions",
x = "Mean Hydroxymethylation Levels",
y = "Density",
color = "Group"
) +
scale_color_manual(
values = c(
"Escapees_hESC" = "#FF7D00",
"Escapees_D4hPGCLC" = "#FFB347",
"Non Escapees_hESC" = "#15616D",
"Non Escapees_D4hPGCLC" = "#1E91A1"
)
) +
coord_cartesian(xlim = c(0, 0.25), ylim = c(0, 60)) + # cut off at y = 6
theme_minimal()
5.Session Info
date()
## [1] "Wed Jul 9 13:15:31 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] xfun_0.39 crayon_1.5.2
## [23] RCurl_1.98-1.12 jsonlite_1.8.7
## [25] annotatr_1.24.0 iterators_1.0.14
## [27] glue_1.6.2 gtable_0.3.3
## [29] zlibbioc_1.44.0 XVector_0.38.0
## [31] DelayedArray_0.24.0 Rhdf5lib_1.20.0
## [33] HDF5Array_1.26.0 scales_1.2.1
## [35] DBI_1.1.3 rngtools_1.5.2
## [37] Rcpp_1.0.11 xtable_1.8-4
## [39] progress_1.2.2 bumphunter_1.40.0
## [41] bit_4.0.5 httr_1.4.6
## [43] RColorBrewer_1.1-3 ellipsis_0.3.2
## [45] pkgconfig_2.0.3 XML_3.99-0.14
## [47] R.methodsS3_1.8.2 farver_2.1.1
## [49] sass_0.4.7 dbplyr_2.3.3
## [51] locfit_1.5-9.7 utf8_1.2.3
## [53] tidyselect_1.2.0 labeling_0.4.2
## [55] rlang_1.1.1 reshape2_1.4.4
## [57] later_1.3.1 AnnotationDbi_1.60.2
## [59] munsell_0.5.0 BiocVersion_3.16.0
## [61] tools_4.2.1 cachem_1.0.8
## [63] cli_3.6.1 generics_0.1.3
## [65] RSQLite_2.3.1 evaluate_0.21
## [67] stringr_1.5.0 fastmap_1.1.1
## [69] yaml_2.3.7 outliers_0.15
## [71] knitr_1.43 bit64_4.0.5
## [73] purrr_1.0.1 KEGGREST_1.38.0
## [75] nlme_3.1-162 doRNG_1.8.6
## [77] sparseMatrixStats_1.10.0 mime_0.12
## [79] R.oo_1.25.0 xml2_1.3.5
## [81] biomaRt_2.54.1 compiler_4.2.1
## [83] rstudioapi_0.15.0 filelock_1.0.2
## [85] curl_5.0.1 png_0.1-8
## [87] interactiveDisplayBase_1.36.0 tibble_3.2.1
## [89] bslib_0.5.0 stringi_1.7.12
## [91] highr_0.10 GenomicFeatures_1.50.4
## [93] lattice_0.21-8 Matrix_1.6-0
## [95] permute_0.9-7 vctrs_0.6.3
## [97] pillar_1.9.0 lifecycle_1.0.3
## [99] rhdf5filters_1.10.1 BiocManager_1.30.20
## [101] jquerylib_0.1.4 data.table_1.14.8
## [103] bitops_1.0-7 httpuv_1.6.11
## [105] rtracklayer_1.58.0 R6_2.5.1
## [107] BiocIO_1.8.0 promises_1.2.0.1
## [109] codetools_0.2-19 gtools_3.9.4
## [111] rhdf5_2.42.1 rjson_0.2.21
## [113] withr_2.5.0 regioneR_1.30.0
## [115] GenomicAlignments_1.34.1 Rsamtools_2.14.0
## [117] GenomeInfoDbData_1.2.9 parallel_4.2.1
## [119] hms_1.1.3 grid_4.2.1
## [121] rmarkdown_2.23 DelayedMatrixStats_1.20.0
## [123] shiny_1.7.4.1 restfulr_0.0.15