Comparison escapees vs non-escapees - list provided by Mitinori lab
Gaja Matassa
Date: May 05, 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("~/DataDir/1.PreliminaryAnalysis/Output/WithSelectedEGCLCs/bsseq_obj_sharedby75ofall.rds") #getting bsseq object where sample selection and CpG filtering were already performed
bsseq_obj_escapees <- readRDS(paste0(params$InputFolder, "bsseq_obj_escapees_Murase.rds"))
bsseq_obj_NonEscapees <- readRDS(paste0(params$InputFolder, "bsseq_obj_NonEscapees_Murase.rds"))- Our bsseq object contain 3680616 CpG loci.
- Among them 207478 CpG loci belong to escapees regions (escapee CpGs).
- While 3473138 CpG loci do not belong to escapees regions (non-escapee CpGs).
Colors for figures
cellTypes_colors <- c(hiPSCs ="#dd1c77", iMeLCs ="#377eb8", hPGCLCs ="#4daf4a", hEGCLCs = "#ff7f00")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)
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)
As expected it resembles the PCA done on the complete one, since we removed only 207478 out of a total of 3680616 CpG loci.
2.Generation of 500 random sets of 207478 non-escapee CpGs and distributions plots
Let’s generate random indices to randomly select 207478 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 iMeLC_rep1 iMeLC_rep2
## [1,] 0.8421053 0.8571429 0.8260870 0.8400000 0.7500000 0.8000000
## [2,] 1.0000000 1.0000000 1.0000000 1.0000000 0.9600000 0.9411765
## [3,] 1.0000000 0.9545455 1.0000000 0.9615385 0.9642857 1.0000000
## [4,] 1.0000000 0.9565217 0.9615385 0.9629630 0.8148148 0.9756098
## [5,] 1.0000000 0.9545455 0.9629630 0.9310345 0.9642857 0.9523810
## [6,] 1.0000000 0.9565217 0.9642857 0.9000000 0.9615385 1.0000000
## iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,] 0.8709677 0.8421053 0.6500000 0.6875 0.7931034 0.8888889
## [2,] 0.9655172 0.9444444 1.0000000 1.0000 1.0000000 0.9782609
## [3,] 0.9375000 0.8636364 0.8636364 1.0000 0.9655172 0.9074074
## [4,] 0.9375000 0.8750000 0.8571429 1.0000 0.8928571 0.9454545
## [5,] 0.9375000 0.9200000 0.8333333 0.9375 0.9259259 0.9272727
## [6,] 0.9375000 0.8518519 0.8260870 1.0000 0.8666667 0.9464286
## hEGCLC_rep2 hEGCLC_rep3
## [1,] 0.7209302 0.7297297
## [2,] 0.9729730 0.9411765
## [3,] 0.9523810 0.9565217
## [4,] 0.9347826 0.8936170
## [5,] 0.9583333 0.8750000
## [6,] 0.9607843 0.9591837Escapee CpG methylation distribution plot
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
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
getMeth(random_sets[[1]], type = "raw") %>% head()## hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 iMeLC_rep1 iMeLC_rep2
## [1,] 0.82142857 0.7321429 0.74666667 0.67500000 0.79710145 0.68595041
## [2,] 0.00000000 0.0000000 0.00000000 0.00000000 0.00000000 0.01666667
## [3,] 0.94117647 0.9787234 0.92857143 1.00000000 0.85000000 0.89473684
## [4,] 0.85714286 1.0000000 0.81818182 0.85714286 0.93333333 0.73913043
## [5,] 0.04444444 0.0000000 0.02985075 0.01492537 0.02222222 0.01923077
## [6,] 0.47916667 0.3906250 0.53424658 0.52631579 0.43939394 0.50000000
## iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,] 0.7341772 0.7319588 0.71666667 0.48387097 0.67272727 0.58823529
## [2,] 0.0000000 0.0000000 0.00000000 0.00000000 0.02857143 0.01818182
## [3,] 0.8787879 0.9534884 0.91891892 0.85714286 0.84615385 0.87301587
## [4,] 0.7500000 0.6666667 0.38095238 0.80952381 0.73076923 0.80769231
## [5,] 0.0000000 0.0000000 0.03333333 0.06896552 0.00000000 0.04950495
## [6,] 0.4200000 0.4871795 0.48076923 0.32000000 0.52000000 0.40243902
## hEGCLC_rep2 hEGCLC_rep3
## [1,] 0.66165414 0.63333333
## [2,] 0.00000000 0.00000000
## [3,] 0.79104478 0.89743590
## [4,] 0.90625000 0.83870968
## [5,] 0.04651163 0.03921569
## [6,] 0.50961538 0.53703704First set of randomly selected non-escapee CpG methylation distribution plot
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
getMeth(random_sets[[2]], type = "raw") %>% head()## hiPSC_rep1 hiPSC_rep2 hiPSC_rep3 hiPSC_rep4 iMeLC_rep1 iMeLC_rep2
## [1,] 0.01492537 0.00000000 0.02564103 0.009009009 0.01538462 0.01342282
## [2,] 0.58064516 0.73684211 0.75000000 0.783783784 0.68000000 0.62790698
## [3,] 0.97727273 0.95312500 0.98437500 0.970588235 0.97435897 0.93396226
## [4,] 0.00000000 0.01298701 0.00000000 0.000000000 0.00000000 0.00000000
## [5,] 0.89655172 0.87500000 0.83333333 0.878787879 0.86363636 0.77049180
## [6,] 0.85714286 0.93478261 0.85714286 0.936363636 0.92424242 0.87500000
## iMeLC_rep3 iMeLC_rep4 hPGCLC_rep1 hPGCLC_rep2 hPGCLC_rep3 hEGCLC_rep1
## [1,] 0.02061856 0.04545455 0.0000000 0.02500000 0.01612903 0.0000000
## [2,] 0.63333333 0.68965517 0.5312500 0.61111111 0.56756757 0.3400000
## [3,] 0.90909091 1.00000000 0.9189189 1.00000000 0.97872340 0.9411765
## [4,] 0.00000000 0.00000000 0.0000000 0.04878049 0.00000000 0.0000000
## [5,] 0.82758621 0.77419355 0.8250000 0.80000000 0.70731707 0.8979592
## [6,] 0.86111111 0.84158416 0.5614035 0.68571429 0.68421053 0.3629032
## hEGCLC_rep2 hEGCLC_rep3
## [1,] 0.007246377 0.013698630
## [2,] 0.266666667 0.294117647
## [3,] 0.956043956 0.976744186
## [4,] 0.000000000 0.008264463
## [5,] 0.793650794 0.898305085
## [6,] 0.672413793 0.669421488Second set of randomly selected non-escapee CpG methylation distribution plot
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 hEGCLC_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_rep1"], getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_rep1"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(bsseq_obj_escapees, type = "raw")[, "hEGCLC_rep1"] and getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_rep1"]
## D = 0.46225, p-value < 2.2e-16
## alternative hypothesis: two-sidedExample of KS test result comparing distributions in case of hEGCLC_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_rep1"], getMeth(random_sets[[2]], type = "raw")[, "hEGCLC_rep1"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(random_sets[[1]], type = "raw")[, "hEGCLC_rep1"] and getMeth(random_sets[[2]], type = "raw")[, "hEGCLC_rep1"]
## D = 0.0018763, p-value = 0.8585
## alternative hypothesis: two-sidedExample of KS test result comparing distributions: escapee CpG methylation levels (hiPSC_rep2 sample) vs first randomly selected non-escapee set (hPGCLC_rep2 sample).
ks_test_result <- ks.test(getMeth(random_sets[[1]], type = "raw")[, "hiPSC_rep2"], getMeth(random_sets[[1]], type = "raw")[, "hPGCLC_rep2"])
print(ks_test_result)##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: getMeth(random_sets[[1]], type = "raw")[, "hiPSC_rep2"] and getMeth(random_sets[[1]], type = "raw")[, "hPGCLC_rep2"]
## D = 0.12154, 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(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.49739, 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.0025256, p-value = 0.5224
## 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 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1How 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] "Mon May 5 11:50:32 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