Generation of bsseq object with escapee CpGs provided by MitinoriLab (Murase et al. paper)
Gaja Matassa
Date: May 05, 2025
1. Environment setting
library(ggplot2)
library(bsseq)
## 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
## 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: '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(scales)
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:bsseq':
##
## combine
## The following object is masked from 'package:Biobase':
##
## combine
## The following object is masked from 'package:BiocGenerics':
##
## combine
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following object is masked from 'package:bsseq':
##
## combine
## The following object is masked from 'package:Biobase':
##
## combine
## The following object is masked from 'package:matrixStats':
##
## count
## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
##
## intersect
## The following objects are masked from 'package:IRanges':
##
## collapse, desc, intersect, setdiff, slice, union
## The following objects are masked from 'package:S4Vectors':
##
## first, intersect, rename, setdiff, setequal, union
## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(dmrseq)
library(annotatr)
library(sechm)
library(rtracklayer)
source("../Methylation_helper.R")
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
##
## expandOutputFolder <- params$OutputFolder2. Loading data
2.1 BSSeq object
Loading of bsseq object with CpGs shared by 75% of the samples.
bsseq_obj <- readRDS("~/DataDir/1.PreliminaryAnalysis/Output/WithSelectedEGCLCs/bsseq_obj_sharedby75ofall.rds") #getting bsseq object where sample selection and CpG filtering were already performed 2.3 Escapees regions
These escapees regions were provided by Mitinori lab and are both from the Tang et al. 2015 and Murase et al. papers.
library(data.table)
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:SummarizedExperiment':
##
## shift
## The following object is masked from 'package:GenomicRanges':
##
## shift
## The following object is masked from 'package:IRanges':
##
## shift
## The following objects are masked from 'package:S4Vectors':
##
## first, second
EscapeesRegions_F <- fread(
"~/DataDir/7.EscapeesExploration/Escapees_Murase/hypermrF5_F1cF1csF2cF2aghGCF_hypermr_merge_sum_homer_label.bed",
sep = "\t",
header = FALSE,
fill = TRUE # fills in missing columns with NA
)
EscapeesRegions_M <- fread(
"~/DataDir/7.EscapeesExploration/Escapees_Murase/hypermrM4_M1_M1s_M2_hGC_hypermr_merge_sum_homer_label.bed",
sep = "\t",
header = FALSE,
fill = TRUE # fills in missing columns with NA
)head(EscapeesRegions_M)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14
## 1: chr1 10468 10620 M1 M1_single - - 1 chr1 10469 10620 + 0 NA
## 2: chr1 10688 10903 - - M2 - 2 chr1 10689 10903 + 0 NA
## 3: chr1 12756 14521 M1 - M2 - 3 chr1 12757 14521 + 0 NA
## 4: chr1 14652 16141 M1 M1_single M2 hGC_M 4 chr1 14653 16141 + 0 NA
## 5: chr1 73528 77244 M1 M1_single M2 - 5 chr1 73529 77244 + 0 NA
## 6: chr1 180998 185231 M1 M1_single M2 hGC_M 6 chr1 180999 185231 + 0 NA
## V15 V16 V17 V18
## 1: Intergenic TAR1|Satellite|telo Satellite|telo TAR1
## 2: Intergenic TAR1|Satellite|telo Satellite|telo TAR1
## 3: TTS (NR_024540) TTS | TTS
## 4: TTS (NR_046018) TTS | TTS
## 5: Intergenic L1PA2|LINE|L1 LINE|L1 L1PA2
## 6: intron (NR_148357, intron 1 of 2) intron | intron
## V19 V20 V21 V22
## 1: Satellite telo -1330 NR_046018
## 2: Satellite telo -1078 NR_046018
## 3: 1765 NR_046018
## 4: 2039 NR_107062
## 5: LINE L1 6295 NM_001005484
## 6: 727 NR_148357
EscapeesRegions_M <- EscapeesRegions_M[, c(1,2,3,15,16,18,19,20)]
colnames(EscapeesRegions_M) <- c("Chr","Start","End","GenomicType","HomerAnno","Gene","Repetitive","Repetitive2")Let’s now focus on escapees found in the male line. We have 192309 escapees regions.
These regions are in GRCh38.p12 assembly and with 0-based coordinates
EscapeesRegions_GRanges <- makeGRangesFromDataFrame(EscapeesRegions_M, seqnames.field = "Chr", start.field = "Start", end.field = "End", keep.extra.columns = FALSE, na.rm=TRUE, starts.in.df.are.0based=TRUE)3. Generation of BSSeq oject with escapee CpGs and non-escapee CpGs
bsseq_ranges_df <- granges(bsseq_obj) %>% as.data.frame() #granges has 0-based position they seem 1-based. Look here for explanation: https://www.biostars.org/p/84686/
bsseq_ranges_df$end <- bsseq_ranges_df$start + 1 #to make them 0-based
bsseq_granges <- makeGRangesFromDataFrame(bsseq_ranges_df, starts.in.df.are.0based = TRUE)
bsseq_obj_1based <- bsseq_obj
bsseq_obj_1based@rowRanges <- bsseq_granges
granges(bsseq_obj_1based)
## GRanges object with 3680616 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 10469 *
## [2] chr1 10471 *
## [3] chr1 10484 *
## [4] chr1 10489 *
## [5] chr1 10493 *
## ... ... ... ...
## [3680612] chrM 11647 *
## [3680613] chrM 11689 *
## [3680614] chrM 11692 *
## [3680615] chrM 11710 *
## [3680616] chrM 11716 *
## -------
## seqinfo: 25 sequences from an unspecified genome; no seqlengthsThis to obtain a bsseq object with only CpGs contained in the escapees regions
m <- findOverlaps(granges(bsseq_obj_1based), EscapeesRegions_GRanges)
bsseq_obj_escapees <- bsseq_obj_1based[queryHits(m),]With our targeted kit we are able to cover 207478 CpGs belonging to the escapees regions.
Let’s create also a bsseq object with non-escapee CpGs (all the other CpGs that are not in escapees regions)
bsseq_obj_NonEscapees <- bsseq_obj_1based[-queryHits(m),]4. Violin plot showing methylation levels of escapee CpGs
assays(bsseq_obj_escapees)$Meth <- getMeth(bsseq_obj_escapees, type = "raw")
assays(bsseq_obj_NonEscapees)$Meth <- getMeth(bsseq_obj_NonEscapees, type = "raw")
sample_colors <- c("hiPSC_rep1" = "#dd1c77", "hiPSC_rep2" = "#dd1c77", "hiPSC_rep3" = "#dd1c77", "hiPSC_rep4" = "#dd1c77",
"iMeLC_rep1" = "#377eb8", "iMeLC_rep2" = "#377eb8", "iMeLC_rep3" = "#377eb8", "iMeLC_rep4" = "#377eb8",
"hPGCLC_rep1" = "#4daf4a", "hPGCLC_rep2" = "#4daf4a", "hPGCLC_rep3" = "#4daf4a",
"hEGCLC_rep1" = "#ff7f00", "hEGCLC_rep2" = "#ff7f00", "hEGCLC_rep3" = "#ff7f00")
violin_plot(assays(bsseq_obj_escapees)$Meth, stat="median") + ggplot2::scale_fill_manual(values = sample_colors)
## No id variables; using all as measure variables
5. Saving and Session Info
Saving:
DataFrame of Escapee regions
bsseq object subset with only CpGs in escapees regions
bsseq object subset with non-escapee CpGs (all the other CpGs that are not in escapees regions)
openxlsx::write.xlsx(data.frame(
seqnames = as.character(seqnames(EscapeesRegions_GRanges)),
start = start(EscapeesRegions_GRanges),
end = end(EscapeesRegions_GRanges),
strand = as.character(strand(EscapeesRegions_GRanges)),
mcols(EscapeesRegions_GRanges)
), file = paste0(OutputFolder, "EscapeeRegions_Murase_hg38.xlsx"))
saveRDS(bsseq_obj_escapees, paste0(OutputFolder, "bsseq_obj_escapees_Murase.rds"))
saveRDS(bsseq_obj_NonEscapees, paste0(OutputFolder, "bsseq_obj_NonEscapees_Murase.rds"))date()
## [1] "Mon May 5 10:58:34 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] data.table_1.14.8 tidyr_1.3.0
## [3] rtracklayer_1.58.0 sechm_1.6.0
## [5] annotatr_1.24.0 dmrseq_1.18.1
## [7] dplyr_1.1.2 gridExtra_2.3
## [9] scales_1.2.1 bsseq_1.34.0
## [11] SummarizedExperiment_1.28.0 Biobase_2.58.0
## [13] MatrixGenerics_1.10.0 matrixStats_1.0.0
## [15] GenomicRanges_1.50.2 GenomeInfoDb_1.34.9
## [17] IRanges_2.32.0 S4Vectors_0.36.2
## [19] BiocGenerics_0.44.0 ggplot2_3.4.2
##
## loaded via a namespace (and not attached):
## [1] circlize_0.4.15 AnnotationHub_3.6.0
## [3] BiocFileCache_2.6.1 plyr_1.8.8
## [5] splines_4.2.1 BiocParallel_1.32.6
## [7] digest_0.6.33 ca_0.71.1
## [9] foreach_1.5.2 htmltools_0.5.5
## [11] fansi_1.0.4 magrittr_2.0.3
## [13] memoise_2.0.1 BSgenome_1.66.3
## [15] cluster_2.1.4 doParallel_1.0.17
## [17] openxlsx_4.2.5.2 tzdb_0.4.0
## [19] limma_3.54.2 ComplexHeatmap_2.14.0
## [21] Biostrings_2.66.0 readr_2.1.4
## [23] R.utils_2.12.2 prettyunits_1.1.1
## [25] colorspace_2.1-0 blob_1.2.4
## [27] rappdirs_0.3.3 xfun_0.39
## [29] crayon_1.5.2 RCurl_1.98-1.12
## [31] jsonlite_1.8.7 iterators_1.0.14
## [33] glue_1.6.2 registry_0.5-1
## [35] gtable_0.3.3 zlibbioc_1.44.0
## [37] XVector_0.38.0 V8_4.3.0
## [39] GetoptLong_1.0.5 DelayedArray_0.24.0
## [41] Rhdf5lib_1.20.0 shape_1.4.6
## [43] HDF5Array_1.26.0 DBI_1.1.3
## [45] rngtools_1.5.2 randomcoloR_1.1.0.1
## [47] Rcpp_1.0.11 xtable_1.8-4
## [49] progress_1.2.2 clue_0.3-64
## [51] bumphunter_1.40.0 bit_4.0.5
## [53] httr_1.4.6 RColorBrewer_1.1-3
## [55] ellipsis_0.3.2 farver_2.1.1
## [57] pkgconfig_2.0.3 XML_3.99-0.14
## [59] R.methodsS3_1.8.2 sass_0.4.7
## [61] dbplyr_2.3.3 locfit_1.5-9.7
## [63] utf8_1.2.3 labeling_0.4.2
## [65] tidyselect_1.2.0 rlang_1.1.1
## [67] reshape2_1.4.4 later_1.3.1
## [69] AnnotationDbi_1.60.2 munsell_0.5.0
## [71] BiocVersion_3.16.0 tools_4.2.1
## [73] cachem_1.0.8 cli_3.6.1
## [75] generics_0.1.3 RSQLite_2.3.1
## [77] evaluate_0.21 stringr_1.5.0
## [79] fastmap_1.1.1 yaml_2.3.7
## [81] outliers_0.15 knitr_1.43
## [83] bit64_4.0.5 zip_2.3.0
## [85] purrr_1.0.1 KEGGREST_1.38.0
## [87] nlme_3.1-162 doRNG_1.8.6
## [89] sparseMatrixStats_1.10.0 mime_0.12
## [91] R.oo_1.25.0 xml2_1.3.5
## [93] biomaRt_2.54.1 compiler_4.2.1
## [95] rstudioapi_0.15.0 filelock_1.0.2
## [97] curl_5.0.1 png_0.1-8
## [99] interactiveDisplayBase_1.36.0 tibble_3.2.1
## [101] bslib_0.5.0 stringi_1.7.12
## [103] highr_0.10 GenomicFeatures_1.50.4
## [105] lattice_0.21-8 Matrix_1.6-0
## [107] permute_0.9-7 vctrs_0.6.3
## [109] pillar_1.9.0 lifecycle_1.0.3
## [111] rhdf5filters_1.10.1 BiocManager_1.30.20
## [113] GlobalOptions_0.1.2 jquerylib_0.1.4
## [115] bitops_1.0-7 seriation_1.4.2
## [117] httpuv_1.6.11 R6_2.5.1
## [119] BiocIO_1.8.0 TSP_1.2-4
## [121] promises_1.2.0.1 codetools_0.2-19
## [123] gtools_3.9.4 rhdf5_2.42.1
## [125] rjson_0.2.21 withr_2.5.0
## [127] regioneR_1.30.0 GenomicAlignments_1.34.1
## [129] Rsamtools_2.14.0 GenomeInfoDbData_1.2.9
## [131] parallel_4.2.1 hms_1.1.3
## [133] grid_4.2.1 rmarkdown_2.23
## [135] DelayedMatrixStats_1.20.0 Rtsne_0.16
## [137] shiny_1.7.4.1 restfulr_0.0.15