DMR Annotation
Sarah Stucchi
Date: August 19, 2025
1. Environment setting
library(DSS)
## Loading required package: Biobase
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: BiocParallel
## Loading required package: bsseq
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## The following object is masked from 'package:Biobase':
##
## rowMedians
## Loading required package: parallel
library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:IRanges':
##
## slice
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(ggplot2)
library(scales)
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(ChIPseeker)
##
## ChIPseeker v1.34.1 For help: https://guangchuangyu.github.io/software/ChIPseeker
##
## If you use ChIPseeker in published research, please cite:
## Qianwen Wang, Ming Li, Tianzhi Wu, Li Zhan, Lin Li, Meijun Chen, Wenqin Xie, Zijing Xie, Erqiang Hu, Shuangbin Xu, Guangchuang Yu. Exploring epigenomic datasets by ChIPseeker. Current Protocols 2022, 2(10): e585
library(org.Hs.eg.db)
## Loading required package: AnnotationDbi
##
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:plotly':
##
## select
##
library(ReactomePA)
## ReactomePA v1.42.0 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use ReactomePA in published research, please cite:
## Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Molecular BioSystems 2016, 12(2):477-479
library(clusterProfiler)
## clusterProfiler v4.6.2 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
##
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:AnnotationDbi':
##
## select
## The following object is masked from 'package:IRanges':
##
## slice
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following object is masked from 'package:stats':
##
## filter
library(dmrseq)
library(annotatr)
library(DT)
library(data.table)
##
## Attaching package: 'data.table'
## 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
GTFFile <- "/group/testa/Project/PGCLC/NewData_fromLaSitua/3.TwistBedAnn/Input/gencode.v35.annotation.gtf.gz"
txdb_v35 <- GenomicFeatures::makeTxDbFromGFF(GTFFile, format="gtf")
## Import genomic features from the file as a GRanges object ...
## OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ...
## Warning in .get_cds_IDX(mcols0$type, mcols0$phase): The "phase" metadata column contains non-NA values for features of type
## stop_codon. This information was ignored.
## OKInputFolder <- params$InputFolder
OutputFolder <- params$OutputFolder2. Loading data
2.1 BSSeq objects
Loading of bsseq object, with CpGs shared by 75% of samples (after sample selection).
bsseq_obj <- readRDS("/group/testa/Project/PGCLC/emseq_rebuttal/2.PreliminaryAnalysis/bsseq_obj_sharedby75ofall.rds") #getting bsseq object where sample selection and CpG filtering were already performed 2.2 DMRs and DMLs
The DMRs here selected are called like this:
delta=0.25 and p.threshold=0.0001
DMRs <- readRDS(paste0(InputFolder, '/', 'DMRs.rds'))
#DMLs <- readRDS(paste0(InputFolder, '/', 'DMLs.rds'))2.3 Universe of genes
Importing Twist regions that were annotated
TwistAnnotated <- readRDS(params$TwistAnnotated)
TwistAnnotated %>% head()
## seqnames start end width strand
## 1 chr1 10466 10585 120 *
## 2 chr1 10790 10909 120 *
## 3 chr1 15806 15925 120 *
## 4 chr1 18768 18887 120 *
## 5 chr1 29357 29476 120 *
## 6 chr1 36544 36663 120 *
## ann
## 1 cg14817997,cpg_inter
## 2 cg16269199,cg26928153,cpg_inter
## 3 CTCF_binding_site,cg13869341,cpg_inter,promoter_flanking_region
## 4 cg14008030,cpg_inter
## 5 cg12045430,cg20826792,cpg_islands,promoter
## 6 cg18231760,cpg_inter,promoter_flanking_region
## annotation geneChr geneStart geneEnd geneLength geneStrand
## 1 Promoter (1-2kb) 1 11869 14409 2541 1
## 2 Promoter (<=1kb) 1 11869 14409 2541 1
## 3 Promoter (1-2kb) 1 17369 17436 68 2
## 4 Promoter (1-2kb) 1 17369 17436 68 2
## 5 Promoter (<=1kb) 1 29554 31097 1544 1
## 6 Promoter (<=1kb) 1 34554 36081 1528 2
## ensembl_gene_id_version transcriptId distanceToTSS
## 1 ENSG00000223972.5 ENST00000456328.2 -1284
## 2 ENSG00000223972.5 ENST00000456328.2 -960
## 3 ENSG00000278267.1 ENST00000619216.1 1511
## 4 ENSG00000278267.1 ENST00000619216.1 -1332
## 5 ENSG00000243485.5 ENST00000473358.1 -78
## 6 ENSG00000237613.2 ENST00000417324.1 -463
## flank_txIds
## 1 ENST00000456328.2;ENST00000450305.2;ENST00000488147.1
## 2 ENST00000456328.2;ENST00000450305.2;ENST00000488147.1
## 3 ENST00000456328.2;ENST00000450305.2;ENST00000488147.1;ENST00000619216.1
## 4 ENST00000456328.2;ENST00000488147.1;ENST00000619216.1
## 5 ENST00000488147.1;ENST00000473358.1;ENST00000469289.1;ENST00000607096.1
## 6 ENST00000417324.1;ENST00000461467.1
## flank_geneIds
## 1 ENSG00000223972.5;ENSG00000223972.5;ENSG00000227232.5
## 2 ENSG00000223972.5;ENSG00000223972.5;ENSG00000227232.5
## 3 ENSG00000223972.5;ENSG00000223972.5;ENSG00000227232.5;ENSG00000278267.1
## 4 ENSG00000223972.5;ENSG00000227232.5;ENSG00000278267.1
## 5 ENSG00000227232.5;ENSG00000243485.5;ENSG00000243485.5;ENSG00000284332.1
## 6 ENSG00000237613.2;ENSG00000237613.2
## flank_gene_distances ensembl_gene_id hgnc_symbol external_gene_name
## 1 -1284;-1425;0 ENSG00000223972 DDX11L1 DDX11L1
## 2 -960;-1101;0 ENSG00000223972 DDX11L1 DDX11L1
## 3 3937;3796;0;1511 ENSG00000278267 MIR6859-1 MIR6859-1
## 4 6899;0;-1332 ENSG00000278267 MIR6859-1 MIR6859-1
## 5 0;-78;-791;-890 ENSG00000243485 MIR1302-2HG MIR1302-2HG
## 6 -463;-471 ENSG00000237613 FAM138A FAM138A
## gene_biotype
## 1 transcribed_unprocessed_pseudogene
## 2 transcribed_unprocessed_pseudogene
## 3 miRNA
## 4 miRNA
## 5 lncRNA
## 6 lncRNA
## description
## 1 DEAD/H-box helicase 11 like 1 (pseudogene) [Source:HGNC Symbol;Acc:HGNC:37102]
## 2 DEAD/H-box helicase 11 like 1 (pseudogene) [Source:HGNC Symbol;Acc:HGNC:37102]
## 3 microRNA 6859-1 [Source:HGNC Symbol;Acc:HGNC:50039]
## 4 microRNA 6859-1 [Source:HGNC Symbol;Acc:HGNC:50039]
## 5 MIR1302-2 host gene [Source:HGNC Symbol;Acc:HGNC:52482]
## 6 family with sequence similarity 138 member A [Source:HGNC Symbol;Acc:HGNC:32334]
## chromosome_name start_position end_position
## 1 1 11869 14409
## 2 1 11869 14409
## 3 1 17369 17436
## 4 1 17369 17436
## 5 1 29554 31109
## 6 1 34554 36081Genes_df <- TwistAnnotated[, c("ensembl_gene_id_version", "ensembl_gene_id", "hgnc_symbol", "gene_biotype")] %>% dplyr::distinct()%>%filter(!is.na(hgnc_symbol))%>%filter(!hgnc_symbol%in%"")
rownames(Genes_df) <- NULLThese are the genes with the same symbol but different ensembl id
Genes_df[Genes_df$hgnc_symbol%in%Genes_df[duplicated(Genes_df$hgnc_symbol), "hgnc_symbol"],]
## ensembl_gene_id_version ensembl_gene_id hgnc_symbol
## 3309 ENSG00000285053.1 ENSG00000285053 TBCE
## 3310 ENSG00000284770.2 ENSG00000284770 TBCE
## 20449 ENSG00000237940.3 ENSG00000237940 LINC01238
## 20451 ENSG00000261186.2 ENSG00000261186 LINC01238
## 21924 ENSG00000206195.11 ENSG00000206195 DUXAP8
## 21926 ENSG00000271672.1 ENSG00000271672 DUXAP8
## 24581 ENSG00000284862.3 ENSG00000284862 CCDC39
## 24582 ENSG00000145075.13 ENSG00000145075 CCDC39
## 30308 ENSG00000272655.2 ENSG00000272655 POLR2J4
## 30309 ENSG00000214783.9 ENSG00000214783 POLR2J4
## gene_biotype
## 3309 protein_coding
## 3310 protein_coding
## 20449 lncRNA
## 20451 lncRNA
## 21924 lncRNA
## 21926 transcribed_processed_pseudogene
## 24581 protein_coding
## 24582 lncRNA
## 30308 transcribed_unprocessed_pseudogene
## 30309 lncRNAGenes_df$entrez_gene_id <- as.vector(mapIds(org.Hs.eg.db, Genes_df$hgnc_symbol, "ENTREZID","SYMBOL"))
## 'select()' returned 1:many mapping between keys and columnsnrow(Genes_df[is.na(Genes_df$entrez_gene_id),]) #number of genes without an entrez id
## [1] 421GeneUniverse_symbol <- unique(Genes_df$hgnc_symbol)
GeneUniverse_entrezid <- na.omit(unique(Genes_df$entrez_gene_id))3. Transform DataFrame with DMRs info to GRanges object
Remember that position form bedGraph files are 0-based while R works
with 1-based position and therefore also the packages we will go to use,
this is the reason why I set starts.in.df.are.0based = TRUE
in creating the GRanges objects, it converts 0-based to 1-based
coordinates.
DMRs_GRanges <- DMRs
#DMLs_GRanges <- DMLs
for(i in 1:length(DMRs)){
DMRs_GRanges[[i]] <- makeGRangesFromDataFrame(DMRs_GRanges[[i]], keep.extra.columns = TRUE,
start.field = "start", end.field = "end", starts.in.df.are.0based = TRUE)
#DMLs_GRanges[[i]] <- makeGRangesFromDataFrame(DMLs_GRanges[[i]], keep.extra.columns = TRUE,
#start.field = "start", end.field = "end", starts.in.df.are.0based = TRUE)
} - For hiPSCs vs hEGCLC_p10 we have 125 DMRs.
- For iMeLCs vs hEGCLCs_p10 we have 267 DMRs.
- For hPGCLCs vs hEGCLC_p10 we have 3361 DMRs.
- For hEGCLCs_p5 vs hEGCLCs_p10 we have 770 DMRs.
4. DMRs annotation
4.1 hiPSCs vs hEGCLCs_p10
4.1.1 Genetic annotation and visualization with
ChIPseeker
DMRs are now annotated to their associated genes and genomic regions
using the annotatePeaks function of the
ChIPseekerR package, using a TxDb object. The GenomicFeatures
package uses TxDb objects to store transcript metadata.
This class maps the 5’ and 3’ untranslated regions (UTRs), protein
coding sequences (CDSs) and exons for a set of mRNA transcripts to their
associated genome. All TxDb objects are backed by a SQLite
database that manages genomic locations and the relationships between
pre-processed mRNA transcripts, exons, protein coding sequences, and
their related gene identifiers.
TxDb created by us
peakAnnoList <- lapply(list("Higher methylation in hiPSCs" = DMRs_GRanges$hEGCLCs_p10vshiPSCs[DMRs_GRanges$hEGCLCs_p10vshiPSCs$diff.Methy<0,], "Higher methylation in hEGCLCs_p10" = DMRs_GRanges$hEGCLCs_p10vshiPSCs[DMRs_GRanges$hEGCLCs_p10vshiPSCs$diff.Methy>0,]), annotatePeak, TxDb=txdb_v35, tssRegion=c(-3000, 3000), verbose=FALSE)ChIPseeker::plotAnnoBar(peakAnnoList)
ChIPseeker::plotDistToTSS(peakAnnoList)
Given a list of gene set, compareCluster function will
compute profiles of each gene cluster.
genes = lapply(peakAnnoList, function(i) unique(as.data.frame(i)$geneId))
names(genes) = sub("_", "\n", names(genes))any(duplicated(genes$`Higher methylation in hiPSCs`))
## [1] FALSEany(duplicated(genes$`Higher methylation in hEGCLCs_p10`))
## [1] FALSEThere are no duplicated ensembl ids
table(genes$`Higher methylation in hEGCLCs_p10` %in% Genes_df$ensembl_gene_id_version)
## < table of extent 0 >genes$`Higher methylation in hEGCLCs_p10` [!genes$`Higher methylation in hEGCLCs_p10` %in% Genes_df$ensembl_gene_id_version]
## NULLtable(genes$`Higher methylation in hiPSCs`%in% Genes_df$ensembl_gene_id_version)
##
## FALSE TRUE
## 16 82DMG_uphiPSCsvshEGCLCs_p10 <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hiPSCs`, "entrez_gene_id"]
DMG_uphiPSCsvshEGCLCs_p10 <- na.omit(DMG_uphiPSCsvshEGCLCs_p10[!DMG_uphiPSCsvshEGCLCs_p10 %in% ""])
DMG_downhiPSCsvshEGCLCs_p10 <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hEGCLCs_p10`, "entrez_gene_id"]
DMG_downhiPSCsvshEGCLCs_p10 <- na.omit(DMG_downhiPSCsvshEGCLCs_p10[!DMG_downhiPSCsvshEGCLCs_p10 %in% ""])any(duplicated(DMG_uphiPSCsvshEGCLCs_p10))
## [1] FALSEany(duplicated(DMG_downhiPSCsvshEGCLCs_p10))
## [1] FALSEgenes$`Higher methylation in hiPSCs` <- DMG_uphiPSCsvshEGCLCs_p10
genes$`Higher methylation in hEGCLCs_p10` <- DMG_downhiPSCsvshEGCLCs_p10compKEGG <- compareCluster(geneCluster = genes,
fun = "enrichKEGG",
organism="hsa",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
universe = GeneUniverse_entrezid)
## Warning in compareCluster(geneCluster = genes, fun = "enrichKEGG", organism =
## "hsa", : No enrichment found in any of gene cluster, please check your input...There are no significant terms (They have Padj>0.05)
#dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")compGO <- compareCluster(geneCluster = genes,
fun = "enrichGO",
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
OrgDb='org.Hs.eg.db',
universe=GeneUniverse_entrezid,
ont = "ALL")
## Warning in compareCluster(geneCluster = genes, fun = "enrichGO", keyType =
## "ENTREZID", : No enrichment found in any of gene cluster, please check your
## input...#dotplot(compGO, showCategory = 15, title = "GO Enrichment Analysis")A new DMRs_annotated object is created in which
annotation from ChIPseeker is added.
DMRs_annotated <- DMRs_GRanges
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)$ChipSeekerAnn <- NA
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)$ChipSeekerDistanceToTss <- NA
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)$ChipSeekerGeneId <- NA
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vshiPSCs, peakAnnoList[["Higher methylation in hiPSCs"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hiPSCs"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hiPSCs"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hiPSCs"]]@anno)[subjectHits(m), "geneId"]
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vshiPSCs, peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vshiPSCs)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "geneId"]Adding entrez gene id, gene symbol and gene information
DMRs_annotated$hEGCLCs_p10vshiPSCs <- as.data.frame(DMRs_annotated$hEGCLCs_p10vshiPSCs)
DMRs_annotated$hEGCLCs_p10vshiPSCs <- dplyr::left_join(DMRs_annotated$hEGCLCs_p10vshiPSCs, Genes_df, by=c('ChipSeekerGeneId' = 'ensembl_gene_id_version'))DMRs_annotated$hEGCLCs_p10vshiPSCs[!is.na(DMRs_annotated$hEGCLCs_p10vshiPSCs$hgnc_symbol),] %>% select(-c("ChipSeekerDistanceToTss", "ChipSeekerGeneId", "ensembl_gene_id", "entrez_gene_id")) %>% distinct() %>%
datatable(class = 'hover', rownames = FALSE, caption="DMRs - hEGCLCs_p10vshiPSCs", extension='Buttons', escape = FALSE,
options = list(pageLength=20, dom='Bfrtip', autoWidth=TRUE,
buttons=list(c('csv', 'excel'))))saveRDS(DMRs_annotated, paste0(OutputFolder, '/', 'DMRsAnnotated1.rds'))4.2 iMeLCs VS hEGCLC_p10
4.2.1 Genetic annotation and visualization with
ChIPseeker
DMRs are now annotated to their associated genes and genomic regions
using the annotatePeaks function of the
ChIPseekerR package, using a TxDb object. The GenomicFeatures
package uses TxDb objects to store transcript metadata.
This class maps the 5’ and 3’ untranslated regions (UTRs), protein
coding sequences (CDSs) and exons for a set of mRNA transcripts to their
associated genome. All TxDb objects are backed by a SQLite
database that manages genomic locations and the relationships between
pre-processed mRNA transcripts, exons, protein coding sequences, and
their related gene identifiers.
TxDb created by us
peakAnnoList <- lapply(list("Higher methylation in iMeLCs" = DMRs_GRanges$hEGCLCs_p10vsiMeLCs[DMRs_GRanges$hEGCLCs_p10vsiMeLCs$diff.Methy<0,], "Higher methylation in hEGCLCs_p10" = DMRs_GRanges$hEGCLCs_p10vsiMeLCs[DMRs_GRanges$hEGCLCs_p10vsiMeLCs$diff.Methy>0,]), annotatePeak, TxDb=txdb_v35, tssRegion=c(-3000, 3000), verbose=FALSE)ChIPseeker::plotAnnoBar(peakAnnoList)
ChIPseeker::plotDistToTSS(peakAnnoList)
Given a list of gene set, compareCluster function will
compute profiles of each gene cluster.
genes = lapply(peakAnnoList, function(i) unique(as.data.frame(i)$geneId))
names(genes) = sub("_", "\n", names(genes))any(duplicated(genes$`Higher methylation in hEGCLC_p10`))
## [1] FALSEany(duplicated(genes$`Higher methylation in iMeLCs`))
## [1] FALSEThere are no duplicated ensembl ids
table(genes$`Higher methylation in iMeLCs` %in% Genes_df$ensembl_gene_id_version)
##
## FALSE TRUE
## 17 73genes$`Higher methylation in iMeLCs` [!genes$`Higher methylation in iMeLCs` %in% Genes_df$ensembl_gene_id_version]
## [1] "ENSG00000251196.1" "ENSG00000259607.1" "ENSG00000249149.2"
## [4] "ENSG00000204850.4" "ENSG00000261842.1" "ENSG00000233820.2"
## [7] "ENSG00000274219.1" "ENSG00000270889.1" "ENSG00000255108.1"
## [10] "ENSG00000288571.1" "ENSG00000267713.1" "ENSG00000275201.1"
## [13] "ENSG00000224765.1" "ENSG00000269388.1" "ENSG00000238110.1"
## [16] "ENSG00000274660.1" "ENSG00000227459.1"table(genes$`Higher methylation in hEGCLC_p10`%in% Genes_df$ensembl_gene_id_version)
## < table of extent 0 >DMG_uphEGCLCs_p10vsiMeLCs <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hEGCLC_p10`, "entrez_gene_id"]
DMG_uphEGCLCs_p10vsiMeLCs <- na.omit(DMG_uphEGCLCs_p10vsiMeLCs[!DMG_uphEGCLCs_p10vsiMeLCs %in% ""])
DMG_downhEGCLCs_p10vsiMeLCs <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in iMeLCs`, "entrez_gene_id"]
DMG_downhEGCLCs_p10vsiMeLCs <- na.omit(DMG_downhEGCLCs_p10vsiMeLCs[!DMG_downhEGCLCs_p10vsiMeLCs %in% ""])any(duplicated(DMG_uphEGCLCs_p10vsiMeLCs))
## [1] FALSEany(duplicated(DMG_downhEGCLCs_p10vsiMeLCs))
## [1] FALSEgenes$`Higher methylation in hEGCLC_p10` <- DMG_uphEGCLCs_p10vsiMeLCs
genes$`Higher methylation in iMeLCs` <- DMG_downhEGCLCs_p10vsiMeLCscompKEGG <- compareCluster(geneCluster = genes,
fun = "enrichKEGG",
organism="hsa",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
universe = GeneUniverse_entrezid)
## Warning in compareCluster(geneCluster = genes, fun = "enrichKEGG", organism =
## "hsa", : No enrichment found in any of gene cluster, please check your input...There are no significant terms (They have Padj>0.05)
if (!is.null(compKEGG)) {
dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")
}compGO <- compareCluster(geneCluster = genes,
fun = "enrichGO",
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
OrgDb='org.Hs.eg.db',
universe=GeneUniverse_entrezid,
ont = "ALL")
## Warning in compareCluster(geneCluster = genes, fun = "enrichGO", keyType =
## "ENTREZID", : No enrichment found in any of gene cluster, please check your
## input...#dotplot(compGO, showCategory = 15, title = "GO Enrichment Analysis")A new DMRs_annotated object is created in which
annotation from ChIPseeker is added.
DMRs_annotated <- DMRs_GRanges
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)$ChipSeekerAnn <- NA
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)$ChipSeekerDistanceToTss <- NA
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)$ChipSeekerGeneId <- NA
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vsiMeLCs, peakAnnoList[["Higher methylation in iMeLCs"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in iMeLCs"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in iMeLCs"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in iMeLCs"]]@anno)[subjectHits(m), "geneId"]
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vsiMeLCs, peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vsiMeLCs)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "geneId"]Adding entrez gene id, gene symbol and gene information
DMRs_annotated$hEGCLCs_p10vsiMeLCs <- as.data.frame(DMRs_annotated$hEGCLCs_p10vsiMeLCs)
DMRs_annotated$hEGCLCs_p10vsiMeLCs <- dplyr::left_join(DMRs_annotated$hEGCLCs_p10vsiMeLCs, Genes_df, by=c('ChipSeekerGeneId' = 'ensembl_gene_id_version'))DMRs_annotated$hEGCLCs_p10vsiMeLCs[!is.na(DMRs_annotated$hEGCLCs_p10vsiMeLCs$hgnc_symbol),] %>% select(-c("ChipSeekerDistanceToTss", "ChipSeekerGeneId", "ensembl_gene_id", "entrez_gene_id")) %>% distinct() %>%
datatable(class = 'hover', rownames = FALSE, caption="DMRs - iMeLCsvshEGCLC_p10", extension='Buttons', escape = FALSE,
options = list(pageLength=20, dom='Bfrtip', autoWidth=TRUE,
buttons=list(c('csv', 'excel'))))saveRDS(DMRs_annotated, paste0(OutputFolder, '/', 'DMRsAnnotated2.rds'))4.3 hPGCLCs vs hEGCLCs_p10
4.3.1 Genetic annotation and visualization with
ChIPseeker
DMRs are now annotated to their associated genes and genomic regions
using the annotatePeaks function of the
ChIPseekerR package, using a TxDb object. The GenomicFeatures
package uses TxDb objects to store transcript metadata.
This class maps the 5’ and 3’ untranslated regions (UTRs), protein
coding sequences (CDSs) and exons for a set of mRNA transcripts to their
associated genome. All TxDb objects are backed by a SQLite
database that manages genomic locations and the relationships between
pre-processed mRNA transcripts, exons, protein coding sequences, and
their related gene identifiers.
TxDb created by us
peakAnnoList <- lapply(list("Higher methylation in hPGCLCs" = DMRs_GRanges$hEGCLCs_p10vshPGCLCs[DMRs_GRanges$hEGCLCs_p10vshPGCLCs$diff.Methy<0,], "Higher methylation in hEGCLC_p10" = DMRs_GRanges$hEGCLCs_p10vshPGCLCs[DMRs_GRanges$hEGCLCs_p10vshPGCLCs$diff.Methy>0,]), annotatePeak, TxDb=txdb_v35, tssRegion=c(-3000, 3000), verbose=FALSE)ChIPseeker::plotAnnoBar(peakAnnoList)
ChIPseeker::plotDistToTSS(peakAnnoList)
Given a list of gene set, compareCluster function will
compute profiles of each gene cluster.
genes = lapply(peakAnnoList, function(i) unique(as.data.frame(i)$geneId))
names(genes) = sub("_", "\n", names(genes))any(duplicated(genes$`Higher methylation in hPGCLCs`))
## [1] FALSEany(duplicated(genes$`Higher methylation in hEGCLC_p10`))
## [1] FALSEThere are no duplicated ensembl ids
table(genes$`Higher methylation in hEGCLC_p10` %in% Genes_df$ensembl_gene_id_version)
## < table of extent 0 >genes$`Higher methylation in hEGCLC_p10` [!genes$`Higher methylation in hEGCLC_p10` %in% Genes_df$ensembl_gene_id_version]
## NULLtable(genes$`Higher methylation in hPGCLCs`%in% Genes_df$ensembl_gene_id_version)
##
## FALSE TRUE
## 34 157DMG_uphPGCLCsvshEGCLC_p10 <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hPGCLCs`, "entrez_gene_id"]
DMG_uphPGCLCsvshEGCLC_p10 <- na.omit(DMG_uphPGCLCsvshEGCLC_p10[!DMG_uphPGCLCsvshEGCLC_p10 %in% ""])
DMG_downhPGCLCsvshEGCLC_p10 <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hEGCLC_p10`, "entrez_gene_id"]
DMG_downhPGCLCsvshEGCLC_p10 <- na.omit(DMG_downhPGCLCsvshEGCLC_p10[!DMG_downhPGCLCsvshEGCLC_p10 %in% ""])any(duplicated(DMG_uphPGCLCsvshEGCLC_p10))
## [1] FALSEany(duplicated(DMG_downhPGCLCsvshEGCLC_p10))
## [1] FALSEgenes$`Higher methylation in hPGCLCs` <- DMG_uphPGCLCsvshEGCLC_p10
genes$`Higher methylation in hEGCLC_p10` <- DMG_downhPGCLCsvshEGCLC_p10compKEGG <- compareCluster(geneCluster = genes,
fun = "enrichKEGG",
organism="hsa",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
universe = GeneUniverse_entrezid)
## Warning in compareCluster(geneCluster = genes, fun = "enrichKEGG", organism =
## "hsa", : No enrichment found in any of gene cluster, please check your input...There are no significant terms (They have Padj>0.05)
#dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")compGO <- compareCluster(geneCluster = genes,
fun = "enrichGO",
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
OrgDb='org.Hs.eg.db',
universe=GeneUniverse_entrezid,
ont = "ALL")dotplot(compGO, showCategory = 15, title = "GO Enrichment Analysis")
A new DMRs_annotated object is created in which
annotation from ChIPseeker is added.
DMRs_annotated <- DMRs_GRanges
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)$ChipSeekerAnn <- NA
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)$ChipSeekerDistanceToTss <- NA
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)$ChipSeekerGeneId <- NA
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vshPGCLCs, peakAnnoList[["Higher methylation in hPGCLCs"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hPGCLCs"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hPGCLCs"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hPGCLCs"]]@anno)[subjectHits(m), "geneId"]
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vshPGCLCs, peakAnnoList[["Higher methylation in hEGCLC_p10"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hEGCLC_p10"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hEGCLC_p10"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vshPGCLCs)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hEGCLC_p10"]]@anno)[subjectHits(m), "geneId"]Adding entrez gene id, gene symbol and gene information
DMRs_annotated$hEGCLCs_p10vshPGCLCs <- as.data.frame(DMRs_annotated$hEGCLCs_p10vshPGCLCs)
DMRs_annotated$hEGCLCs_p10vshPGCLCs <- dplyr::left_join(DMRs_annotated$hEGCLCs_p10vshPGCLCs, Genes_df, by=c('ChipSeekerGeneId' = 'ensembl_gene_id_version'))DMRs_annotated$hEGCLCs_p10vshPGCLCs[!is.na(DMRs_annotated$hEGCLCs_p10vshPGCLCs$hgnc_symbol),] %>% select(-c("ChipSeekerDistanceToTss", "ChipSeekerGeneId", "ensembl_gene_id", "entrez_gene_id")) %>% distinct() %>%
datatable(class = 'hover', rownames = FALSE, caption="DMRs - hEGCLCs_p10vshPGCLCs", extension='Buttons', escape = FALSE,
options = list(pageLength=20, dom='Bfrtip', autoWidth=TRUE,
buttons=list(c('csv', 'excel'))))saveRDS(DMRs_annotated, paste0(OutputFolder, '/', 'DMRsAnnotated3.rds'))4.4 hEGCLCs_p5 vs hEGCLCs_p10
4.4.1 Genetic annotation and visualization with
ChIPseeker
DMRs are now annotated to their associated genes and genomic regions
using the annotatePeaks function of the
ChIPseekerR package, using a TxDb object. The GenomicFeatures
package uses TxDb objects to store transcript metadata.
This class maps the 5’ and 3’ untranslated regions (UTRs), protein
coding sequences (CDSs) and exons for a set of mRNA transcripts to their
associated genome. All TxDb objects are backed by a SQLite
database that manages genomic locations and the relationships between
pre-processed mRNA transcripts, exons, protein coding sequences, and
their related gene identifiers.
TxDb created by us
peakAnnoList <- lapply(list("Higher methylation in hEGCLCs_p5" = DMRs_GRanges$hEGCLCs_p10vshEGCLCs_p5[DMRs_GRanges$hEGCLCs_p10vshEGCLCs_p5$diff.Methy<0,], "Higher methylation in hEGCLCs_p10" = DMRs_GRanges$hEGCLCs_p10vshEGCLCs_p5[DMRs_GRanges$hEGCLCs_p10vshEGCLCs_p5$diff.Methy>0,]), annotatePeak, TxDb=txdb_v35, tssRegion=c(-3000, 3000), verbose=FALSE)ChIPseeker::plotAnnoBar(peakAnnoList)
ChIPseeker::plotDistToTSS(peakAnnoList)
Given a list of gene set, compareCluster function will
compute profiles of each gene cluster.
genes = lapply(peakAnnoList, function(i) unique(as.data.frame(i)$geneId))
names(genes) = sub("_", "\n", names(genes))any(duplicated(genes$`Higher methylation in hEGCLCs_p5`))
## [1] FALSEany(duplicated(genes$`Higher methylation in hEGCLCs_p10`))
## [1] FALSEThere are no duplicated ensembl ids
table(genes$`Higher methylation in hEGCLCs_p10` %in% Genes_df$ensembl_gene_id_version)
## < table of extent 0 >genes$`Higher methylation in hEGCLCs_p10` [!genes$`Higher methylation in hEGCLCs_p10` %in% Genes_df$ensembl_gene_id_version]
## NULLtable(genes$`Higher methylation in hEGCLCs_p5`%in% Genes_df$ensembl_gene_id_version)
## < table of extent 0 >DMG_uphEGCLCs_p5vshEGCLCs_p10 <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hEGCLCs_p5`, "entrez_gene_id"]
DMG_uphEGCLCs_p5vshEGCLCs_p10 <- na.omit(DMG_uphEGCLCs_p5vshEGCLCs_p10[!DMG_uphEGCLCs_p5vshEGCLCs_p10 %in% ""])
DMG_downhEGCLCs_p5vshEGCLCs_p10 <- Genes_df[Genes_df$ensembl_gene_id_version %in% genes$`Higher methylation in hEGCLCs_p10`, "entrez_gene_id"]
DMG_downhEGCLCs_p5vshEGCLCs_p10 <- na.omit(DMG_downhEGCLCs_p5vshEGCLCs_p10[!DMG_downhEGCLCs_p5vshEGCLCs_p10 %in% ""])any(duplicated(DMG_uphEGCLCs_p5vshEGCLCs_p10))
## [1] FALSEany(duplicated(DMG_downhEGCLCs_p5vshEGCLCs_p10))
## [1] FALSEgenes$`Higher methylation in hEGCLCs_p5` <- DMG_uphEGCLCs_p5vshEGCLCs_p10
genes$`Higher methylation in hEGCLCs_p10` <- DMG_downhEGCLCs_p5vshEGCLCs_p10compKEGG <- compareCluster(geneCluster = genes,
fun = "enrichKEGG",
organism="hsa",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
universe = GeneUniverse_entrezid)
## Warning in compareCluster(geneCluster = genes, fun = "enrichKEGG", organism =
## "hsa", : No enrichment found in any of gene cluster, please check your input...There are no significant terms (They have Padj>0.05)
#dotplot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")compGO <- compareCluster(geneCluster = genes,
fun = "enrichGO",
keyType = "ENTREZID",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
OrgDb='org.Hs.eg.db',
universe=GeneUniverse_entrezid,
ont = "ALL")
## Warning in compareCluster(geneCluster = genes, fun = "enrichGO", keyType =
## "ENTREZID", : No enrichment found in any of gene cluster, please check your
## input...#dotplot(compGO, showCategory = 15, title = "GO Enrichment Analysis")A new DMRs_annotated object is created in which
annotation from ChIPseeker is added.
DMRs_annotated <- DMRs_GRanges
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)$ChipSeekerAnn <- NA
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)$ChipSeekerDistanceToTss <- NA
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)$ChipSeekerGeneId <- NA
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5, peakAnnoList[["Higher methylation in hEGCLCs_p5"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p5"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p5"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p5"]]@anno)[subjectHits(m), "geneId"]
m <- findOverlaps(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5, peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)[queryHits(m), "ChipSeekerAnn"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "annotation"]
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)[queryHits(m), "ChipSeekerDistanceToTss"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "distanceToTSS"]
mcols(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)[queryHits(m), "ChipSeekerGeneId"] <- mcols(peakAnnoList[["Higher methylation in hEGCLCs_p10"]]@anno)[subjectHits(m), "geneId"]Adding entrez gene id, gene symbol and gene information
DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5 <- as.data.frame(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5)
DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5 <- dplyr::left_join(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5, Genes_df, by=c('ChipSeekerGeneId' = 'ensembl_gene_id_version'))DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5[!is.na(DMRs_annotated$hEGCLCs_p10vshEGCLCs_p5$hgnc_symbol),] %>% select(-c("ChipSeekerDistanceToTss", "ChipSeekerGeneId", "ensembl_gene_id", "entrez_gene_id")) %>% distinct() %>%
datatable(class = 'hover', rownames = FALSE, caption="DMRs - hEGCLCs_p10vshEGCLCs_p5", extension='Buttons', escape = FALSE,
options = list(pageLength=20, dom='Bfrtip', autoWidth=TRUE,
buttons=list(c('csv', 'excel'))))saveRDS(DMRs_annotated, paste0(OutputFolder, '/', 'DMRsAnnotated4.rds'))5. Saving DMRs with annotation and Session Info
saveRDS(Genes_df, paste0(OutputFolder, '/', 'GeneUniverse.rds'))
saveRDS(DMRs_annotated, paste0(OutputFolder, '/', 'DMRsAnnotated_final.rds'))SessionInfo <- sessionInfo()
Date <- date()Date
## [1] "Tue Aug 19 12:31:42 2025"
SessionInfo
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] data.table_1.14.8 DT_0.28
## [3] annotatr_1.24.0 dmrseq_1.18.1
## [5] clusterProfiler_4.6.2 ReactomePA_1.42.0
## [7] org.Hs.eg.db_3.16.0 AnnotationDbi_1.60.2
## [9] ChIPseeker_1.34.1 gridExtra_2.3
## [11] scales_1.2.1 plotly_4.10.2
## [13] ggplot2_3.4.2 DSS_2.46.0
## [15] bsseq_1.34.0 SummarizedExperiment_1.28.0
## [17] MatrixGenerics_1.10.0 matrixStats_1.0.0
## [19] GenomicRanges_1.50.2 GenomeInfoDb_1.34.9
## [21] IRanges_2.32.0 S4Vectors_0.36.2
## [23] BiocParallel_1.32.6 Biobase_2.58.0
## [25] BiocGenerics_0.44.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.3
## [2] R.utils_2.12.2
## [3] tidyselect_1.2.0
## [4] RSQLite_2.3.1
## [5] htmlwidgets_1.6.2
## [6] grid_4.2.1
## [7] scatterpie_0.1.8
## [8] munsell_0.5.0
## [9] codetools_0.2-19
## [10] withr_2.5.0
## [11] colorspace_2.1-0
## [12] GOSemSim_2.24.0
## [13] filelock_1.0.2
## [14] highr_0.10
## [15] knitr_1.43
## [16] rstudioapi_0.15.0
## [17] DOSE_3.24.2
## [18] labeling_0.4.2
## [19] GenomeInfoDbData_1.2.9
## [20] polyclip_1.10-4
## [21] bit64_4.0.5
## [22] farver_2.1.1
## [23] rhdf5_2.42.1
## [24] downloader_0.4
## [25] vctrs_0.6.3
## [26] treeio_1.23.1
## [27] generics_0.1.3
## [28] gson_0.1.0
## [29] xfun_0.39
## [30] BiocFileCache_2.6.1
## [31] regioneR_1.30.0
## [32] R6_2.5.1
## [33] graphlayouts_0.8.4
## [34] locfit_1.5-9.7
## [35] bitops_1.0-7
## [36] rhdf5filters_1.10.1
## [37] cachem_1.0.8
## [38] fgsea_1.24.0
## [39] gridGraphics_0.5-1
## [40] DelayedArray_0.24.0
## [41] promises_1.2.0.1
## [42] BiocIO_1.8.0
## [43] ggraph_2.1.0
## [44] enrichplot_1.18.4
## [45] gtable_0.3.3
## [46] tidygraph_1.2.3
## [47] rlang_1.1.1
## [48] splines_4.2.1
## [49] rtracklayer_1.58.0
## [50] lazyeval_0.2.2
## [51] BiocManager_1.30.20
## [52] yaml_2.3.7
## [53] reshape2_1.4.4
## [54] crosstalk_1.2.0
## [55] GenomicFeatures_1.50.4
## [56] httpuv_1.6.11
## [57] qvalue_2.30.0
## [58] tools_4.2.1
## [59] ggplotify_0.1.0
## [60] ellipsis_0.3.2
## [61] gplots_3.1.3
## [62] jquerylib_0.1.4
## [63] RColorBrewer_1.1-3
## [64] Rcpp_1.0.11
## [65] plyr_1.8.8
## [66] sparseMatrixStats_1.10.0
## [67] progress_1.2.2
## [68] zlibbioc_1.44.0
## [69] purrr_1.0.1
## [70] RCurl_1.98-1.12
## [71] prettyunits_1.1.1
## [72] viridis_0.6.2
## [73] bumphunter_1.40.0
## [74] cowplot_1.1.1
## [75] ggrepel_0.9.3
## [76] magrittr_2.0.3
## [77] reactome.db_1.82.0
## [78] xtable_1.8-4
## [79] mime_0.12
## [80] hms_1.1.3
## [81] patchwork_1.1.2
## [82] evaluate_0.21
## [83] HDO.db_0.99.1
## [84] XML_3.99-0.14
## [85] compiler_4.2.1
## [86] biomaRt_2.54.1
## [87] tibble_3.2.1
## [88] KernSmooth_2.23-22
## [89] crayon_1.5.2
## [90] shadowtext_0.1.2
## [91] R.oo_1.25.0
## [92] htmltools_0.5.5
## [93] tzdb_0.4.0
## [94] later_1.3.1
## [95] ggfun_0.0.9
## [96] tidyr_1.3.0
## [97] aplot_0.1.10
## [98] DBI_1.1.3
## [99] tweenr_2.0.2
## [100] dbplyr_2.3.3
## [101] MASS_7.3-60
## [102] rappdirs_0.3.3
## [103] boot_1.3-28.1
## [104] readr_2.1.4
## [105] Matrix_1.6-0
## [106] permute_0.9-7
## [107] cli_3.6.1
## [108] R.methodsS3_1.8.2
## [109] igraph_1.5.0
## [110] pkgconfig_2.0.3
## [111] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [112] GenomicAlignments_1.34.1
## [113] foreach_1.5.2
## [114] xml2_1.3.5
## [115] ggtree_3.6.2
## [116] bslib_0.5.0
## [117] rngtools_1.5.2
## [118] XVector_0.38.0
## [119] doRNG_1.8.6
## [120] yulab.utils_0.0.6
## [121] stringr_1.5.0
## [122] digest_0.6.33
## [123] graph_1.76.0
## [124] Biostrings_2.66.0
## [125] rmarkdown_2.23
## [126] fastmatch_1.1-3
## [127] tidytree_0.4.2
## [128] DelayedMatrixStats_1.20.0
## [129] restfulr_0.0.15
## [130] curl_5.0.1
## [131] shiny_1.7.4.1
## [132] Rsamtools_2.14.0
## [133] gtools_3.9.4
## [134] graphite_1.44.0
## [135] rjson_0.2.21
## [136] outliers_0.15
## [137] lifecycle_1.0.3
## [138] nlme_3.1-162
## [139] jsonlite_1.8.7
## [140] Rhdf5lib_1.20.0
## [141] viridisLite_0.4.2
## [142] limma_3.54.2
## [143] BSgenome_1.66.3
## [144] fansi_1.0.4
## [145] pillar_1.9.0
## [146] lattice_0.21-8
## [147] KEGGREST_1.38.0
## [148] fastmap_1.1.1
## [149] httr_1.4.6
## [150] plotrix_3.8-2
## [151] GO.db_3.16.0
## [152] interactiveDisplayBase_1.36.0
## [153] glue_1.6.2
## [154] iterators_1.0.14
## [155] png_0.1-8
## [156] BiocVersion_3.16.0
## [157] bit_4.0.5
## [158] ggforce_0.4.1
## [159] stringi_1.7.12
## [160] sass_0.4.7
## [161] HDF5Array_1.26.0
## [162] blob_1.2.4
## [163] AnnotationHub_3.6.0
## [164] caTools_1.18.2
## [165] memoise_2.0.1
## [166] dplyr_1.1.2
## [167] ape_5.7-1