6  Scanning the genome for regions of interest

footprintR provides a flexible scanning framework, which can be used to partition the genome (or a subset thereof) into windows and calculate a custom score for each window. Optionally, neighboring windows with similar scores can then be merged into larger signal regions. In this chapter, we show some examples of how this framework can be used to scan the genome for various types of signals of interest.

The figure below illustrates the general idea of the scanning framework.

Schematic of the footprintR genome scanning framework.
Figure 6.1

The first step is to define the windows and calculate some quantity of interest for each of them, using the quantifyWindowsInRegion function. Next, these quantities are passed to a score function, which summarizes the input data into a score for each window - this can, for example, perform a differential methylation analysis based on the counts of methylated and unmethylated nucleotides obtained from quantifyWindowsInRegion. Finally, the window scores are processed using processWindowScores. This function allows neighboring windows with similar scores to be merged into larger signal regions. The user can choose to either call the functions directly, or use the wrapper function scanForHighScoringRegions.

Thanks to the modularity, the scanning framework is flexible and can accommodate many different use cases by selecting a suitable quantification function, scoring function and aggregation of the scores. The table below lists some of the applications that can be accommodated with the functions that are provided within footprintR. However, the user can also define their own functions for specific use cases.

Table of use cases that can be addressed using the genome scanning framework in footprintR.
Figure 6.2

Below we will show several examples of using the scanning framework for different purposes.

Show/hide code
BSgenomeName <- "BSgenome.Mmusculus.GENCODE.GRCm39.gencodeM34"

library(footprintR)
library(ggplot2)
library(BSgenomeName, character.only = TRUE)
library(BiocParallel)
library(stringr)
library(Hmisc)
library(SummarizedExperiment)
library(patchwork)
library(dplyr)
library(tidyr)
library(parallel)

# Load genome
gnm <- get(BSgenomeName)
genome(gnm) <- "mm39"

# Define parallelization backend
bpParam <- MulticoreParam(workers = 24L, RNGseed = 42L)

6.1 Genome-wide DMR scanning (between samples)

The first use case involves scanning for differentially methylated regions (DMRs) between two conditions, based on 6mA data. In our case, we will compare wild-type mouse ES cells to ‘TKO’ mouse ES cells where the three DNA methyltransferases (Dnmts) have been knocked out. For time reasons, we limit the scanning to chromosome 19 - the most efficient way of performing a genome-wide analysis is typically to parallelize over chromosomes (e.g. by replacing the lapply call in the code below by a parallelized version).

We tile the genome into 30bp windows (windowSize), shifted successively by 15bp (windowStep) to achieve a high resolution and minimize the risk of splitting a region of interest in the middle (thereby potentially losing the signal). We then use the sumNmodNvalid function to count the number of modified As for each window, as well as the total number of As. The classification of nucleotides is obtained by thresholding the modification probability using the modProbThreshold value (0.5). Given the modified and total count for each sample, we next use the getDifferentiallyModifiedWindows score function to perform the differential methylation analysis for each window. This will return several statistics, including the dirNegLog10PValue, which is obtained by calculating -log10(p-value) and multiplying with the sign of the log-fold change. We will use this statistic (scoreCol) as the input for the final processing of the windows. Briefly, the statistics will be smoothed over neighboring windows (the degree of smoothing is controlled by the minperiod argument) and subsequently thresholded (thresh argument). Neighboring windows where the statistic exceeds this threshold will be fused into a larger region of interest, unless they are separated by more than 50 nucleotides (maxGap).

Show/hide code
# define bam files to include
bamfiles <- c("data/mESC_wt_6mA_rep1.bam", "data/mESC_wt_6mA_rep2.bam",
              "data/mESC_TKO_6mA_rep1.bam", "data/mESC_TKO_6mA_rep2.bam")
names(bamfiles) <- sub("\\.bam", "", basename(bamfiles))

# define sample annotation table
sampleAnnot <- data.frame(sample = names(bamfiles),
                          group = factor(str_extract(names(bamfiles), "wt|TKO"), 
                                         levels = c("wt", "TKO")))
Show/hide code
# run scanning of chromosome 19
chrlens <- seqlengths(gnm)["chr19"]
dmrL <- lapply(seq_along(chrlens), function(i) {
    scanForHighScoringRegions(
        bamfiles = bamfiles,
        sampleAnnot = sampleAnnot,
        chromosomeLengths = chrlens[i],
        quantFunction = "sumNmodNvalid",
        scoreFunction = "getDifferentiallyModifiedWindows",
        modbase = "a",
        modProbThreshold = 0.5,
        tileSize = 2e+06,
        seqinfo = seqinfo(gnm),
        sequenceContextWidth = 1,
        sequenceReference = gnm,
        sequenceContext = "A",
        windowMode = "fixed",
        windowSize = 30,
        windowStep = 15,
        scoreCol = "dirNegLog10PValue",
        scoreAction = "smoothFuse",
        thresh = 3,
        minperiod = 3,
        maxGap = 50,
        BPPARAM = bpParam,
        verbose = interactive())
})

Each element of the resulting list (in this case, there is only one) contains the inferred DMRs for one chromosome.

Show/hide code
length(dmrL[[1]])
[1] 32498
Show/hide code
head(dmrL[[1]])
GRanges object with 6 ranges and 5 metadata columns:
      seqnames          ranges strand | dirNegLog10PValueThresh
         <Rle>       <IRanges>  <Rle> |               <numeric>
  [1]    chr19 3051451-3051480      * |                 3.02470
  [2]    chr19 3051841-3051915      * |                 3.49924
  [3]    chr19 3052036-3052110      * |                 5.05907
  [4]    chr19 3052591-3052650      * |                 3.72814
  [5]    chr19 3053191-3053265      * |                 3.83881
  [6]    chr19 3053386-3053445      * |                 3.86745
      numWindowsThresh direction dirNegLog10PValue numWindows
             <integer>  <factor>         <numeric>  <integer>
  [1]                1  positive           3.02470          1
  [2]                4  positive           3.49924          4
  [3]                4  positive           5.05907          4
  [4]                3  positive           3.72814          3
  [5]                4  positive           3.83881          4
  [6]                3  positive           3.86745          3
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
Show/hide code
table(dmrL[[1]]$direction)

positive negative 
   32490        8 

We can visualize one of the top hits using plotRegion.

Show/hide code
reg <- dmrL[[1]][which.max(abs(dmrL[[1]]$dirNegLog10PValue))]
se <- readModBam(bamfiles = bamfiles, sampleAnnot = sampleAnnot, 
                 regions = reg + 500, modbase = "a", 
                 seqinfo = seqinfo(gnm), sequenceContextWidth = 1, 
                 sequenceReference = gnm, trim = TRUE, 
                 verbose = interactive())
se <- flattenReadLevelAssay(se)
plotRegion(se, region = reg + 500, sequenceContext = "A", 
           tracks = list(list(trackData = "mod_prob", trackType = "Heatmap",
                              interpolate = TRUE),
                         list(trackData = "FracMod", trackType = "Smooth",
                              colorBy = "group", highlightRegions = reg))) + 
    plot_layout(heights = c(3, 1))
Figure 6.3

6.1.1 Choosing suitable values for thresholds and smoothing parameters

In the example above, we smooth the window scores and fuse neighboring windows with smoothed scores above a certain threshold (3). The parameters for this step of the analysis were determined by inspecting the intermediate results obtained from applying the quantification and scoring functions independently to a small region of the genome. Here we illustrate how this can be done, using a 5MB region on chromosome 19. First, we call the quantifyWindowsInRegion function, with sumNmodNvalid as the quantFunction, to get the total and modified nucleotide and the fraction of modified nucleotides for each window. The windows will be automatically defined within the specified region, given the indicated window size and step.

Show/hide code
# define the small region
regSmall <- "chr19:31000000-36000000"

# define windows in the small region and quantify data in windows
se <- quantifyWindowsInRegion(bamfiles = bamfiles,
                              region = regSmall,
                              modbase = "a",
                              modProbThreshold = 0.5,
                              sampleAnnot = sampleAnnot,
                              seqinfo = seqinfo(gnm),
                              sequenceContextWidth = 1,
                              sequenceReference = gnm,
                              sequenceContext = "A",
                              windowMode = "fixed",
                              windowSize = 30,
                              windowStep = 15,
                              quantFunction = "sumNmodNvalid",
                              BPPARAM = bpParam,
                              verbose = interactive())

The output of quantifyWindowsInRegion is a SummarizedExperiment object, with assays for the Nmod, Nvalid and FracMod values for each window in each sample.

Show/hide code
se
class: RangedSummarizedExperiment 
dim: 333073 4 
metadata(1): readLevelData
assays(3): Nmod Nvalid FracMod
rownames(333073): 1 2 ... 333331 333332
rowData names(0):
colnames(4): mESC_wt_6mA_rep1 mESC_wt_6mA_rep2 mESC_TKO_6mA_rep1
  mESC_TKO_6mA_rep2
colData names(3): sample modbase group

Next, we apply the score function (getDifferentiallyModifiedWindows) to perform the differential methylation analysis for each window.

Show/hide code
gr <- getDifferentiallyModifiedWindows(se, verbose = interactive())

This function (as all other score functions compatible with the footprintR scanning framework) returns a GRanges object containing the windows, with scores as metadata columns.

Show/hide code
gr
GRanges object with 333073 ranges and 9 metadata columns:
           seqnames            ranges strand |     logFC    logCPM        LR
              <Rle>         <IRanges>  <Rle> | <numeric> <numeric> <numeric>
       [1]    chr19 31000000-31000029      * |  0.807304  0.629767   2.45710
       [2]    chr19 31000015-31000044      * |  0.742682  0.697488   2.06101
       [3]    chr19 31000030-31000059      * |  0.636205  0.516105   1.54618
       [4]    chr19 31000045-31000074      * |  0.508839  0.762750   0.95991
       [5]    chr19 31000060-31000089      * |  0.648976  0.770806   1.60876
       ...      ...               ...    ... .       ...       ...       ...
  [333069]    chr19 35999905-35999934      * | 0.1183276 0.0717503 0.0623986
  [333070]    chr19 35999920-35999949      * | 0.2911895 0.5566819 0.3401543
  [333071]    chr19 35999935-35999964      * | 0.0709167 0.5667234 0.0202314
  [333072]    chr19 35999950-35999979      * | 0.2903877 0.2983924 0.3527457
  [333073]    chr19 35999965-35999994      * | 0.7348412 0.7163173 1.8874211
              PValue       FDR dirNegLog10PValue FracMod_wt FracMod_TKO
           <numeric> <numeric>         <numeric>  <numeric>   <numeric>
       [1]  0.116994  0.402987          0.931835   0.129820    0.200790
       [2]  0.151110  0.437962          0.820708   0.145305    0.213280
       [3]  0.213700  0.494005          0.670195   0.141810    0.190575
       [4]  0.327210  0.581068          0.485174   0.136349    0.176522
       [5]  0.204666  0.486337          0.688954   0.165851    0.233447
       ...       ...       ...               ...        ...         ...
  [333069]  0.802744  0.891057         0.0954228  0.3054885    0.314199
  [333070]  0.559740  0.739766         0.2520135  0.2349981    0.259468
  [333071]  0.886893  0.939795         0.0521289  0.2306079    0.229048
  [333072]  0.552563  0.735028         0.2576182  0.1846282    0.212466
  [333073]  0.169493  0.455429         0.7708480  0.0793377    0.125558
           DeltaFracMod
              <numeric>
       [1]    0.0709699
       [2]    0.0679755
       [3]    0.0487643
       [4]    0.0401731
       [5]    0.0675964
       ...          ...
  [333069]   0.00871096
  [333070]   0.02446953
  [333071]  -0.00155955
  [333072]   0.02783778
  [333073]   0.04621999
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

At this point, we can either keep the results on the individual window level, or attempt to fuse neighboring windows with high scores, using the processWindowScores function. We will illustrate the latter using the dirNegLog10PValue as the score column. First, we plot the estimates of this score along the genomic region (just a small part of it, for visibility).

Show/hide code
ggplot(as.data.frame(gr[6700:7200]), aes(x = start, y = dirNegLog10PValue)) + 
    geom_line() + 
    theme_bw()
Figure 6.4

To determine suitable smoothing parameters, we can call the processWindow with scoreAction = "smooth".

Show/hide code
(grsmooth <- processWindowScores(gr, scoreCol = "dirNegLog10PValue", 
                                 scoreAction = "smooth", minperiod = 3))
GRanges object with 333073 ranges and 9 metadata columns:
           seqnames            ranges strand |     logFC    logCPM        LR
              <Rle>         <IRanges>  <Rle> | <numeric> <numeric> <numeric>
       [1]    chr19 31000000-31000029      * |  0.807304  0.629767   2.45710
       [2]    chr19 31000015-31000044      * |  0.742682  0.697488   2.06101
       [3]    chr19 31000030-31000059      * |  0.636205  0.516105   1.54618
       [4]    chr19 31000045-31000074      * |  0.508839  0.762750   0.95991
       [5]    chr19 31000060-31000089      * |  0.648976  0.770806   1.60876
       ...      ...               ...    ... .       ...       ...       ...
  [333069]    chr19 35999905-35999934      * | 0.1183276 0.0717503 0.0623986
  [333070]    chr19 35999920-35999949      * | 0.2911895 0.5566819 0.3401543
  [333071]    chr19 35999935-35999964      * | 0.0709167 0.5667234 0.0202314
  [333072]    chr19 35999950-35999979      * | 0.2903877 0.2983924 0.3527457
  [333073]    chr19 35999965-35999994      * | 0.7348412 0.7163173 1.8874211
              PValue       FDR dirNegLog10PValue FracMod_wt FracMod_TKO
           <numeric> <numeric>         <numeric>  <numeric>   <numeric>
       [1]  0.116994  0.402987          0.877458   0.129820    0.200790
       [2]  0.151110  0.437962          0.764158   0.145305    0.213280
       [3]  0.213700  0.494005          0.661064   0.141810    0.190575
       [4]  0.327210  0.581068          0.706547   0.136349    0.176522
       [5]  0.204666  0.486337          0.911285   0.165851    0.233447
       ...       ...       ...               ...        ...         ...
  [333069]  0.802744  0.891057         0.0687533  0.3054885    0.314199
  [333070]  0.559740  0.739766         0.1374055  0.2349981    0.259468
  [333071]  0.886893  0.939795         0.2544209  0.2306079    0.229048
  [333072]  0.552563  0.735028         0.4199374  0.1846282    0.212466
  [333073]  0.169493  0.455429         0.5510439  0.0793377    0.125558
           DeltaFracMod
              <numeric>
       [1]    0.0709699
       [2]    0.0679755
       [3]    0.0487643
       [4]    0.0401731
       [5]    0.0675964
       ...          ...
  [333069]   0.00871096
  [333070]   0.02446953
  [333071]  -0.00155955
  [333072]   0.02783778
  [333073]   0.04621999
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

This replaces the dirNegLog10PValue column with a smoothed version, which we can see by plotting it overlaid on the unsmoothed values.

Show/hide code
plotdf <- data.frame(
    position = start(gr[6700:7200]), 
    raw = gr$dirNegLog10PValue[6700:7200],
    smoothed = grsmooth$dirNegLog10PValue[6700:7200]
) |>
    pivot_longer(names_to = "type", values_to = "score", -position)
ggplot(plotdf, aes(x = position, y = score, color = type)) + 
    geom_line() +
    scale_color_manual(values = c(raw = "black", smoothed = "red")) + 
    theme_bw()
Figure 6.5

Based on these values, we can now also decide on a suitable value for the threshold parameter, which will be applied to the smoothed values to determine which windows show a strong signal.

6.2 Genome-wide estimation of CpG methylation levels (within samples)

In the example above, the scores were calculated over windows of different size. By setting the window size to 1, we can also calculate single nucleotide-level scores. Here, as an example, we illustrate how to estimate the global distribution of methylation levels across all CpGs in the genome. For time reasons, as above we only consider chromosome 19 here.

In this case, we again use sumNmodNvalid as the quantFunction, as it provides an estimate of the fraction modified bases for each window (in this case, since the window size is 1, for each modified base). We use getRangesWithAssayValues as the scoreFunction, which will move the calculated values from the FracMod and Nvalid assays to the metadata columns of the nucleotide-level GRanges object. By setting scoreAction to “pass”, we instruct the scanning framework to not process the nucleotide-level scores further, but simply return the GRanges object with the added metadata columns. Finally, by setting the sequenceContext to “NCG”, we make sure to only retain positions where the genomic sequence is CG.

Show/hide code
# Define the bam file to use and generate sample annotation table
bamfiles <- c(wt2 = "data/mESC_wt_5mCG_5hmCG_rep2.bam")
sampleAnnot <- data.frame(sample = names(bamfiles))

# Run scanning
chrlens <- seqlengths(gnm)["chr19"]
rescpgL <- lapply(seq_along(chrlens), function(i) {
    scanForHighScoringRegions(
        bamfiles = bamfiles,
        sampleAnnot = sampleAnnot,
        chromosomeLengths = chrlens[i],
        quantFunction = "sumNmodNvalid",
        quantFunctionArgs = list(),
        scoreFunction = "getRangesWithAssayValues",
        scoreFunctionArgs = list(assayName = c("FracMod", "Nvalid")),
        modbase = "m",
        modProbThreshold = 0.5,
        tileSize = 1e6,
        seqinfo = seqinfo(gnm),
        sequenceContextWidth = 3,
        sequenceReference = gnm,
        sequenceContext = "NCG",
        windowMode = "fixed",
        windowSize = 1,
        windowStep = 1,
        scoreCol = paste0("FracMod.", names(bamfiles)[1]),
        scoreAction = "pass",
        BPPARAM = bpParam,
        verbose = FALSE
    )
})

Next, we plot the distribution of modification fractions.

Show/hide code
plotdf <- as.data.frame(unname(rescpgL[[1]])) |>
    tidyr::pivot_longer(cols = -c(seqnames, start, end, width, strand), 
                        names_to = c(".value", "sample"),
                        names_pattern = '(FracMod|Nvalid)\\.(.*)') |> 
    dplyr::mutate(FracModBin = Hmisc::cut2(FracMod,cuts = seq(0, 1, by = 0.1)))
head(plotdf)
# A tibble: 6 × 9
  seqnames   start     end width strand sample FracMod Nvalid FracModBin
  <fct>      <int>   <int> <int> <fct>  <chr>    <dbl>  <dbl> <fct>     
1 chr19    3050119 3050119     1 *      wt2       1         4 [0.9,1.0] 
2 chr19    3050120 3050120     1 *      wt2       1         2 [0.9,1.0] 
3 chr19    3050316 3050316     1 *      wt2       0.75      4 [0.7,0.8) 
4 chr19    3050317 3050317     1 *      wt2       1         2 [0.9,1.0] 
5 chr19    3050603 3050603     1 *      wt2       1         6 [0.9,1.0] 
6 chr19    3050604 3050604     1 *      wt2       1         1 [0.9,1.0] 
Show/hide code
ggplot(plotdf, aes(x = FracModBin)) + 
    geom_bar() + 
    labs(title = "All CpGs", x = "FracMod (binned)") + 
    theme_bw() + 
    theme(axis.text.y = element_text(size = 12),
          axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
          axis.title = element_text(size = 14))
Figure 6.6

6.3 Annotation of predefined windows

As mentioned above, the individual functions called by scanForHighScoringRegions can also be used independently. This can be helpful in different settings - as already shown, it is helpful in order to determine suitable parameters and thresholds for the smoothing and fusing of windows. Another situation when it is useful is when the windows are predefined rather than obtained by tiling the genome into equally sized bins. To illustrate this use case, here we will use the quantifyWindowsInRegion function, together with the phasingScoreFourier quantFunction to annotate a set of CTCF sites with a score quantifying the degree of nucleosome phasing estimated from the average 6mA signal across reads in the region around the motif. This function calculates a phasing score for each window by decomposing the FracMod values into frequency components and extracting the amplitude of the component corresponding to a period consistent with the average nucleosome distance. Note that we set windowMode to "predefined", and define the desired windows via the windows argument. In addition, we use larger windows than above (760bp), to have enough data to be able to detect periodic patterns.

Show/hide code
bamfiles <- c("data/mESC_wt_6mA_rep1.bam")
names(bamfiles) <- sub("\\.bam", "", basename(bamfiles))
sampleAnnot <- data.frame(sample = names(bamfiles), 
                          group = "wt")

period <- 190
windowSize <- 4 * period
numCoef <- 5

(ctcf_sites <- readRDS("data/ctcf_sites_1000.rds"))
GRanges object with 1000 ranges and 1 metadata column:
         seqnames              ranges strand |     score
            <Rle>           <IRanges>  <Rle> | <numeric>
     [1]     chr4 106387210-106387228      - |    21.751
     [2]     chr2   29509735-29509753      + |    16.924
     [3]    chr13 110505258-110505276      - |    19.577
     [4]     chr3 121678692-121678710      + |    18.323
     [5]     chr2   29509695-29509713      + |    17.356
     ...      ...                 ...    ... .       ...
   [996]     chr5   35519412-35519430      + |    19.651
   [997]     chr8   23727858-23727876      + |    11.875
   [998]     chr6 112468049-112468067      + |    16.978
   [999]     chr4   44564745-44564763      + |    22.866
  [1000]     chr8 106101422-106101440      + |    10.054
  -------
  seqinfo: 61 sequences (1 circular) from mm39 genome
Show/hide code
ctcf_sites <- resize(ctcf_sites, width = windowSize, fix = "center")
ctcf_sites <- unstrand(ctcf_sites)

se <- readModBam(bamfiles = bamfiles, sampleAnnot = sampleAnnot, 
                 regions = ctcf_sites, modbase = "a", 
                 seqinfo = seqinfo(gnm), sequenceContextWidth = 1, 
                 sequenceReference = gnm, trim = TRUE, 
                 verbose = interactive())
se <- flattenReadLevelAssay(se)

ctcf_sites$phasing_scores <- unlist(parallel::mclapply(
    seq_along(ctcf_sites),
    function(i) {
        tmp <- phasingScoreFourier(se = subsetByOverlaps(se, ctcf_sites[i]),
                                   gr = ctcf_sites[i], numCoef = numCoef)
        if (nrow(tmp) > 0) {
            phasingscores <- assay(tmp, "phasingScoreAbs")[1, ]
        } else {
            phasingscores <- NA
        }
        # average across samples (in this case only one sample is used)
        mean(phasingscores, na.rm = TRUE)
    }, mc.cores = BiocParallel::bpnworkers(bpParam)))

6.4 Session info

Click to view session info
Show/hide code
sessioninfo::session_info(info = "packages")
═ Session info ═══════════════════════════════════════════════════════════════
─ Packages ───────────────────────────────────────────────────────────────────
 package                                      * version   date (UTC) lib source
 abind                                          1.4-8     2024-09-12 [1] CRAN (R 4.5.0)
 backports                                      1.5.0     2024-05-23 [1] CRAN (R 4.5.0)
 base64enc                                      0.1-3     2015-07-28 [1] CRAN (R 4.5.0)
 Biobase                                      * 2.68.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 BiocGenerics                                 * 0.54.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 BiocIO                                       * 1.18.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 BiocParallel                                 * 1.42.1    2025-06-01 [1] Bioconductor 3.21 (R 4.5.0)
 Biostrings                                   * 2.76.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 bitops                                         1.0-9     2024-10-03 [1] CRAN (R 4.5.0)
 BSgenome                                     * 1.76.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 BSgenome.Mmusculus.GENCODE.GRCm39.gencodeM34 * 0.1.0     2025-04-17 [1] Bioconductor
 checkmate                                      2.3.2     2024-07-29 [1] CRAN (R 4.5.0)
 cli                                            3.6.5     2025-04-23 [1] CRAN (R 4.5.0)
 cluster                                        2.1.8.1   2025-03-12 [2] CRAN (R 4.5.0)
 codetools                                      0.2-20    2024-03-31 [2] CRAN (R 4.5.0)
 colorspace                                     2.1-1     2024-07-26 [1] CRAN (R 4.5.0)
 crayon                                         1.5.3     2024-06-20 [1] CRAN (R 4.5.0)
 curl                                           6.3.0     2025-06-06 [1] CRAN (R 4.5.0)
 data.table                                     1.17.4    2025-05-26 [1] CRAN (R 4.5.0)
 DelayedArray                                   0.34.1    2025-04-17 [1] Bioconductor 3.21 (R 4.5.0)
 dichromat                                      2.0-0.1   2022-05-02 [1] CRAN (R 4.5.0)
 digest                                         0.6.37    2024-08-19 [1] CRAN (R 4.5.0)
 dplyr                                        * 1.1.4     2023-11-17 [1] CRAN (R 4.5.0)
 edgeR                                          4.6.2     2025-05-07 [1] Bioconductor 3.21 (R 4.5.0)
 evaluate                                       1.0.3     2025-01-10 [1] CRAN (R 4.5.0)
 farver                                         2.1.2     2024-05-13 [1] CRAN (R 4.5.0)
 fastmap                                        1.2.0     2024-05-15 [1] CRAN (R 4.5.0)
 footprintR                                   * 0.3.5     2025-06-04 [1] Github (fmicompbio/footprintR@612f713)
 foreign                                        0.8-90    2025-03-31 [2] CRAN (R 4.5.0)
 Formula                                        1.2-5     2023-02-24 [1] CRAN (R 4.5.0)
 generics                                     * 0.1.4     2025-05-09 [1] CRAN (R 4.5.0)
 GenomeInfoDb                                 * 1.44.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 GenomeInfoDbData                               1.2.14    2025-04-17 [1] Bioconductor
 GenomicAlignments                              1.44.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 GenomicRanges                                * 1.60.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 ggforce                                        0.4.2     2024-02-19 [1] CRAN (R 4.5.0)
 ggplot2                                      * 3.5.2     2025-04-09 [1] CRAN (R 4.5.0)
 glue                                           1.8.0     2024-09-30 [1] CRAN (R 4.5.0)
 gridExtra                                      2.3       2017-09-09 [1] CRAN (R 4.5.0)
 gtable                                         0.3.6     2024-10-25 [1] CRAN (R 4.5.0)
 Hmisc                                        * 5.2-3     2025-03-16 [1] CRAN (R 4.5.0)
 htmlTable                                      2.4.3     2024-07-21 [1] CRAN (R 4.5.0)
 htmltools                                      0.5.8.1   2024-04-04 [1] CRAN (R 4.5.0)
 htmlwidgets                                    1.6.4     2023-12-06 [1] CRAN (R 4.5.0)
 httr                                           1.4.7     2023-08-15 [1] CRAN (R 4.5.0)
 IRanges                                      * 2.42.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 jsonlite                                       2.0.0     2025-03-27 [1] CRAN (R 4.5.0)
 knitr                                          1.50      2025-03-16 [1] CRAN (R 4.5.0)
 labeling                                       0.4.3     2023-08-29 [1] CRAN (R 4.5.0)
 lattice                                        0.22-7    2025-04-02 [1] CRAN (R 4.5.0)
 lifecycle                                      1.0.4     2023-11-07 [1] CRAN (R 4.5.0)
 limma                                          3.64.1    2025-05-25 [1] Bioconductor 3.21 (R 4.5.0)
 locfit                                         1.5-9.12  2025-03-05 [1] CRAN (R 4.5.0)
 magrittr                                       2.0.3     2022-03-30 [1] CRAN (R 4.5.0)
 MASS                                           7.3-65    2025-02-28 [2] CRAN (R 4.5.0)
 Matrix                                         1.7-3     2025-03-11 [2] CRAN (R 4.5.0)
 MatrixGenerics                               * 1.20.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 matrixStats                                  * 1.5.0     2025-01-07 [1] CRAN (R 4.5.0)
 nnet                                           7.3-20    2025-01-01 [2] CRAN (R 4.5.0)
 patchwork                                    * 1.3.0     2024-09-16 [1] CRAN (R 4.5.0)
 pillar                                         1.10.2    2025-04-05 [1] CRAN (R 4.5.0)
 pkgconfig                                      2.0.3     2019-09-22 [1] CRAN (R 4.5.0)
 polyclip                                       1.10-7    2024-07-23 [1] CRAN (R 4.5.0)
 purrr                                          1.0.4     2025-02-05 [1] CRAN (R 4.5.0)
 R6                                             2.6.1     2025-02-15 [1] CRAN (R 4.5.0)
 RColorBrewer                                   1.1-3     2022-04-03 [1] CRAN (R 4.5.0)
 Rcpp                                           1.0.14    2025-01-12 [1] CRAN (R 4.5.0)
 RCurl                                          1.98-1.17 2025-03-22 [1] CRAN (R 4.5.0)
 restfulr                                       0.0.15    2022-06-16 [1] CRAN (R 4.5.0)
 rjson                                          0.2.23    2024-09-16 [1] CRAN (R 4.5.0)
 rlang                                          1.1.6     2025-04-11 [1] CRAN (R 4.5.0)
 rmarkdown                                      2.29      2024-11-04 [1] CRAN (R 4.5.0)
 rpart                                          4.1.24    2025-01-07 [2] CRAN (R 4.5.0)
 Rsamtools                                      2.24.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 rstudioapi                                     0.17.1    2024-10-22 [1] CRAN (R 4.5.0)
 rtracklayer                                  * 1.68.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 S4Arrays                                       1.8.1     2025-06-01 [1] Bioconductor 3.21 (R 4.5.0)
 S4Vectors                                    * 0.46.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 scales                                         1.4.0     2025-04-24 [1] CRAN (R 4.5.0)
 sessioninfo                                    1.2.3     2025-02-05 [1] CRAN (R 4.5.0)
 signal                                         1.8-1     2024-06-26 [1] CRAN (R 4.5.0)
 SparseArray                                    1.8.0     2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 statmod                                        1.5.0     2023-01-06 [1] CRAN (R 4.5.0)
 stringi                                        1.8.7     2025-03-27 [1] CRAN (R 4.5.0)
 stringr                                      * 1.5.1     2023-11-14 [1] CRAN (R 4.5.0)
 SummarizedExperiment                         * 1.38.1    2025-04-30 [1] Bioconductor 3.21 (R 4.5.0)
 tibble                                         3.3.0     2025-06-08 [1] CRAN (R 4.5.0)
 tidyr                                        * 1.3.1     2024-01-24 [1] CRAN (R 4.5.0)
 tidyselect                                     1.2.1     2024-03-11 [1] CRAN (R 4.5.0)
 tweenr                                         2.0.3     2024-02-26 [1] CRAN (R 4.5.0)
 UCSC.utils                                     1.4.0     2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 utf8                                           1.2.6     2025-06-08 [1] CRAN (R 4.5.0)
 vctrs                                          0.6.5     2023-12-01 [1] CRAN (R 4.5.0)
 viridisLite                                    0.4.2     2023-05-02 [1] CRAN (R 4.5.0)
 withr                                          3.0.2     2024-10-28 [1] CRAN (R 4.5.0)
 xfun                                           0.52      2025-04-02 [1] CRAN (R 4.5.0)
 XML                                            3.99-0.18 2025-01-01 [1] CRAN (R 4.5.0)
 XVector                                      * 0.48.0    2025-04-15 [1] Bioconductor 3.21 (R 4.5.0)
 yaml                                           2.3.10    2024-07-26 [1] CRAN (R 4.5.0)
 zoo                                            1.8-14    2025-04-10 [1] CRAN (R 4.5.0)

 [1] /tungstenfs/groups/gbioinfo/Appz/R/BioC/R-4.5-release-foss-2024.05_BioC-3.21-release-foss-2024.05
 [2] /tachyon/groups/gbioinfo/Appz/easybuild/software/R/4.5.0-foss-2024.05/lib64/R/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────