1 Introduction

Spontaneous replication errors are not distributed uniformly along the Escherichia coli chromosome. Hasenauer et al. (2025) located them genome-wide by tagging the mismatch-repair protein MutL in rapidly proliferating mutH-deficient cells, in which a mismatch is recognised but cannot be corrected, so that MutL accumulates where errors arise. The resulting ChIP-seq peaks, termed MutL-associated regions (MutL-AR), are enriched for mononucleotide repeats, cruciform-forming DNA, single-stranded DNA and sequences of low thermal stability, and are depleted for GATC methylation sites.

The last of those associations is the one pursued here. Low thermal stability was inferred in that study from sequence composition; melting temperature (Tm) measures it directly, as the temperature at which half the duplex is dissociated, computed from nearest-neighbour stacking thermodynamics rather than from base content alone. Asking whether MutL-AR windows melt at a lower temperature than the rest of the chromosome is therefore a test of the proposed mechanism on its own terms, and one that needs no new experiment: the genome sequence and the published peak coordinates are sufficient.

This vignette computes Tm with TmCalculator across the E. coli K-12 MG1655 chromosome (NCBI assembly GCF_000005845.2 / ASM584v2) for every non-overlapping 200 bp window, the bin size used by Hasenauer et al. themselves, and overlays their annotation tracks, which ship with this package as ecoli_rep_hotspots. Thermodynamic parameters are those of Breslauer et al. (1986) at 50 mM Na+, again matching the published Methods, so that any difference between the two analyses is attributable to what is being compared rather than to how Tm was obtained.

Any window-based profile depends on the bandwidth chosen, and 200 bp is a choice inherited rather than derived. The primary analysis is therefore accompanied by a sensitivity analysis over 50, 100, 200 and 500 bp windows, reported in full in vignette("window_size_sensitivity"). Absolute Tm scales with window length, as duplex thermodynamics requires, but the spatial landscape does not: cross-scale correlations are at least 0.975, and the MutL-AR association, about 1 degree C lower median Tm and about 2.7 percentage points lower GC inside peaks, is reproduced at every window size tested (all p < 1e-12). The conclusions below do not rest on the 200 bp bin.

The analysis proceeds in four steps:

1. Obtain a reference genome as a BSgenome package: use the one Bioconductor ships for the assembly, or forge one when none exists.
2. Compute the Tm profile for every window with tm_calculate(), which tiles the genome, retrieves the sequence and computes Tm in one call.
3. Visualise the profile with plot_genome_track(): circular and linear plots, zoomed views, and multi-omics overlays.
4. Test whether Tm inside MutL-AR peaks differs from the rest of the chromosome, with compare_groups().

2 Prerequisites

library(TmCalculator)
library(BSgenome)
library(GenomicRanges)

3 The reference genome

tm_calculate() retrieves sequence from a BSgenome package when you give it one, so the first question is whether one already exists for your assembly.

If Bioconductor ships your assembly, use it. BSgenome::available.genomes() lists them, and the common references are all there: BSgenome.Hsapiens.UCSC.hg38, BSgenome.Mmusculus.UCSC.mm39, BSgenome.Scerevisiae.UCSC.sacCer3 and so on. Install it and go straight to the next section; nothing in this vignette needs the genome to have been built by hand.

BSgenome::available.genomes()          # is your assembly already packaged?
BiocManager::install("BSgenome.Hsapiens.UCSC.hg38")

Otherwise, forge one. That is the case here. Bioconductor does ship an E. coli package, BSgenome.Ecoli.NCBI.20080805, but it is an older multi-strain snapshot and does not contain the RefSeq assembly this analysis needs (GCF_000005845.2, ASM584v2, sequence U00096.3). The MutL-AR peaks and the microsatellite, cruciform, ssDNA and GATC tracks are all defined against that assembly, so using a different one would silently misplace every feature. BSgenomeForge builds a matching package from an NCBI accession. This is a once-per-assembly step, not a once-per-analysis one.

library(BSgenomeForge)

forgeBSgenomeDataPkgFromNCBI(
  assembly_accession = "GCF_000005845.2",
  pkg_maintainer     = "Junhui Li <ljh.biostat@gmail.com>",
  destdir            = "."
)

install.packages(
  "./BSgenome.Ecoli.NCBI.ASM584v2",
  repos = NULL,
  type  = "source"
)

During TmCalculator package development and vignette building (e.g. R CMD check), we cannot call forgeBSgenomeDataPkgFromNCBI() in the vignette environment, so the chunks below install the pre-forged BSgenome.Ecoli.NCBI.ASM584v2 package from GitHub instead.

ecoli_pkg   <- "BSgenome.Ecoli.NCBI.ASM584v2"
genome_obj  <- "Ecoli"   # BSgenomeObjname in DESCRIPTION; not the package name

.ecoli_genome_ready <- function() {
  if (!requireNamespace(ecoli_pkg, quietly = TRUE)) return(FALSE)
  exists(genome_obj, envir = asNamespace(ecoli_pkg), inherits = FALSE)
}

if (!requireNamespace(ecoli_pkg, quietly = TRUE)) {
  if (!requireNamespace("remotes", quietly = TRUE)) {
    utils::install.packages("remotes", repos = "https://cloud.r-project.org")
  }
  remotes::install_github(
    "JunhuiLi1017/BSgenome.Ecoli.NCBI.ASM584v2",
    upgrade = "never",
    quiet = TRUE
  )
}

if (!.ecoli_genome_ready()) {
  if (!requireNamespace("BiocManager", quietly = TRUE)) {
    utils::install.packages("BiocManager", repos = "https://cloud.r-project.org")
  }
  if (!requireNamespace("BSgenomeForge", quietly = TRUE)) {
    BiocManager::install("BSgenomeForge", ask = FALSE, update = FALSE)
  }
  pkgdir <- BSgenomeForge::forgeBSgenomeDataPkgFromNCBI(
    assembly_accession = "GCF_000005845.2",
    pkg_maintainer     = "Junhui Li <ljh.biostat@gmail.com>",
    destdir            = tempdir()
  )
  utils::install.packages(pkgdir, repos = NULL, type = "source", quiet = TRUE)
  if (ecoli_pkg %in% loadedNamespaces()) {
    unloadNamespace(ecoli_pkg)
  }
}

if (!.ecoli_genome_ready()) {
  stop(
    "Could not load genome object '", genome_obj, "' from package '", ecoli_pkg, "'.\n",
    "Run the manual 'forge' chunk above, or forgeBSgenomeDataPkgFromNCBI() locally.",
    call. = FALSE
  )
}

Load the genome (object Ecoli; package BSgenome.Ecoli.NCBI.ASM584v2 for coordinates):

ecoli_pkg  <- "BSgenome.Ecoli.NCBI.ASM584v2"
genome_obj <- "Ecoli"

suppressPackageStartupMessages(library(ecoli_pkg, character.only = TRUE))
genome      <- base::get(genome_obj, envir = asNamespace(ecoli_pkg))
genome_name <- ecoli_pkg
chr_name    <- "U00096.3"
chr_length  <- length(genome[[chr_name]])

cat("Chromosome:", chr_name, "\n")

## Chromosome: U00096.3

cat("Length:    ", format(chr_length, big.mark = ","), "bp\n")

## Length:     4,641,652 bp

4 Compute the genome-wide Tm profile

tm_calculate() does the whole of it in one call: it tiles the requested regions, retrieves the sequence for each window, and computes a melting temperature for every window. Its first argument is the source the sequence comes from, either the name of a BSgenome package as here or the path to a FASTA file. It takes a package name rather than a loaded object because when the work is spread over several processes each one opens the source itself.

We use nearest-neighbour thermodynamics with the Breslauer et al. (1986) parameter set at 50 mM Na+, matching the Methods of Hasenauer et al. (2025), who computed Tm with this package and the same parameters.

runtime <- system.time({
  tm_ecoli <- tm_calculate(
    input_seq = genome_name,   # a BSgenome package name, a FASTA path, or sequences
    regions  = chr_name,     # the whole chromosome; "U00096.3:1-100000" would
    unit     = "region",     #   tile a sub-region instead
    window   = 200L,
    slide    = 200L,         # slide == window gives a non-overlapping tiling
    method   = "tm_nn",
    nn_table = "DNA_NN_Breslauer_1986",
    Na       = 50,           # mM; standard PCR-like conditions
    verbose  = FALSE
  )$gr                       # $gr is the profile; the object also carries $options
})

cat(sprintf(
  "Tm profile: %.2f s (elapsed) for %s windows\n",
  runtime[["elapsed"]], format(length(tm_ecoli), big.mark = ",")
))

## Tm profile: 1.41 s (elapsed) for 23,208 windows

Tm <- as.data.frame(tm_ecoli[, c("Tm", "GC")])
summary(Tm[, c("Tm", "GC")])

##        Tm               GC       
##  Min.   : 81.05   Min.   :20.50  
##  1st Qu.: 97.07   1st Qu.:47.50  
##  Median : 99.61   Median :52.00  
##  Mean   : 98.91   Mean   :50.79  
##  3rd Qu.:101.49   3rd Qu.:55.00  
##  Max.   :111.36   Max.   :75.00

The result is an ordinary GRanges, so the coordinates, the melting temperatures and the GC percentages travel together through everything that follows.

head(tm_ecoli, 3)

## GRanges object with 3 ranges and 4 metadata columns:
##           seqnames    ranges strand |   region_id             genome_pkg
##              <Rle> <IRanges>  <Rle> | <character>            <character>
##   region1 U00096.3     1-200      + |     region1 BSgenome.Ecoli.NCBI...
##   region2 U00096.3   201-400      + |     region2 BSgenome.Ecoli.NCBI...
##   region3 U00096.3   401-600      + |     region3 BSgenome.Ecoli.NCBI...
##                  GC        Tm
##           <numeric> <numeric>
##   region1      37.0   89.6470
##   region2      51.5   99.8653
##   region3      58.0  103.0685
##   -------
##   seqinfo: 1 sequence from an unspecified genome

A chromosome of 4.6 Mb is small enough to profile in one process. On a mammalian genome the same call takes a BPPARAM and spreads the tiles over workers; vignette("hg38_performance_parallel") covers that, together with the forms regions accepts and how to choose a worker count.

The same function takes sequences directly, which is the form to use when they are already in hand rather than in a genome:

tm_calculate(c("ACGTGCTAGCTAGCTAGC", "GGCCATATATGCGC"), method = "tm_nn", Na = 50)

## GRanges object with 2 ranges and 4 metadata columns:
##     seqnames    ranges strand |           sequence         complement        GC
##        <Rle> <IRanges>  <Rle> |        <character>        <character> <numeric>
##   1        1      1-18      * | ACGTGCTAGCTAGCTAGC TGCACGATCGATCGATCG   55.5556
##   2        2      1-14      * |     GGCCATATATGCGC     CCGGTATATACGCG   57.1429
##            Tm
##     <numeric>
##   1   43.4936
##   2   33.6973
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

5 Visualise with plot_genome_track()

plot_genome_track() is the unified plotting function in TmCalculator. It supports both linear (karyoploteR-based) and circular (base R graphics) layouts from the same track list, with features including ideogram tracks, per-track highlights, proportional track heights, multi-region zoom, and customisable legends.

# Reference labels: replication origin (ori) and terminus (dif)
label <- data.frame(
  seqnames = genome_name,
  start    = c(3925804, 1590777),
  end      = c(3925804, 1590777),
  label    = c("ori", "dif")
)
data(ecoli_rep_hotspots)
tracks <- list(
  # Ideogram: MutL-AR peaks shown inside the chromosome bar
  list(type = "rect", data = ecoli_rep_hotspots$all_peaks_IP_mutH,
       col = "#2C3E50", bg.col = "grey", name = "MutL-AR",
       legend_font_col = "#2C3E50", ideogram = TRUE, height = 0.5),

  # Sequence thermodynamics
  list(type = "line", data = Tm, value_col = "GC",
       name = "GC content", col = "#4A90E2",
       legend_font_col = "#4A90E2"),
  list(type = "line", data = Tm, value_col = "Tm",
       name = "Melting temp", col = "#E06666",
       legend_font_col = "#E06666", height = 2),

  # Repeat / structural features
  list(type = "line", data = ecoli_rep_hotspots$bins_rep,
       value_col = "count", name = "Microsatellites", col = "#2ECC71",
       legend_font_col = "#2ECC71"),
  list(type = "line", data = ecoli_rep_hotspots$bins_cru,
       value_col = "count", name = "Cruciform", col = "#3B3E6B",
       legend_font_col = "#3B3E6B"),

  # ssDNA regions
  list(data = ecoli_rep_hotspots$ssdna, name = "ssDNA",
       col = "#8E44AD", legend_font_col = "#8E44AD"),

  # GATC methylation sites
  list(type = "line", data = ecoli_rep_hotspots$bins_gatc,
       value_col = "count", name = "GATC sites", col = "#D35400",
       legend_font_col = "#D35400"),

  # Global highlight: translucent bands at MutL-AR peaks across all tracks
  list(type = "highlight", data = ecoli_rep_hotspots$all_peaks_IP_mutH,
       col = "#F1C40F", alpha = 0.18)
)

5.2 Circular genome map

plot_genome_track(
  genome_name = genome_name,
  genome_size = chr_length,
  track_list  = tracks,
  circular    = TRUE,
  label       = label
)

Circular genome map of E. coli K-12 MG1655. Concentric rings from outside in: MutL-AR peaks (grey ideogram), GC content, melting temperature, microsatellite density, cruciform sequences, ssDNA regions, and GATC site density. Yellow highlight bands mark MutL-AR peak regions.

5.3 Linear genome map

The same tracks list works for a linear karyoploteR layout — simply omit circular = TRUE. The ideogram track is drawn inside the chromosome bar; all other tracks are stacked as horizontal panels.

plot_genome_track(
  genome_name = genome_name,
  genome_size = chr_length,
  track_list  = tracks
)

Linear genome view of E. coli K-12 MG1655. MutL-AR peaks are drawn inside the chromosome ideogram bar.

5.4 Zoom — single region

The zoom parameter accepts a character string (e.g. "chr:start-end") or a GRanges object. In linear mode, the karyoploteR view is restricted to that region; in circular mode, only data overlapping the region is drawn.

plot_genome_track(
  genome_name = genome_name,
  genome_size = chr_length,
  track_list  = tracks,
  zoom        = "U00096.3:100000-500000",
  ## Seven tracks in a linear panel leave each one little vertical room. The
  ## gap is relative to the panel, so it has to grow with the track count.
  track.gap   = 0.03,
  axis.cex    = 0.55
)

Zoomed linear view of the 0.1-0.5 Mb region. Seven tracks share one panel, so track.gap separates them and axis.cex shrinks the tick labels; without both the y-axis labels of adjacent tracks overlap.

5.5 Zoom — multiple regions

Pass a character vector to zoom to view several disjoint regions at once. In linear mode each region is drawn as a separate stacked panel; in circular mode the regions are concatenated around the circle with small gaps between them.

plot_genome_track(
  genome_name = genome_name,
  genome_size = chr_length,
  track_list  = tracks,
  circular    = TRUE,
  zoom        = c("U00096.3:100000-500000",
                  "U00096.3:3600000-4500000")
)

Two zoomed regions, 0.1-0.5 Mb and 3.6-4.5 Mb, concatenated around the circle with a gap between them.

5.6 Circular canvas panning

Use canvas.xlim and canvas.ylim to pan and magnify a portion of the circular plot. The circle.margin parameter controls whitespace around the plot.

plot_genome_track(
  genome_name   = genome_name,
  genome_size   = chr_length,
  track_list    = tracks,
  circular      = TRUE,
  canvas.xlim   = c(0.5, 1),
  canvas.ylim   = c(0,   1),
  circle.margin = c(0.05, 0.05)
)

Panned circular view showing the upper-right quadrant of the E. coli chromosome.

5.7 Per-track highlights

Individual tracks can carry a highlight field to draw coloured bands within that track only. This is useful for marking specific regions of interest on a per-track basis:

tracks_hl <- tracks
tracks_hl[[4]]$highlight <- list(
  data  = ecoli_rep_hotspots$bins_rep[1100:1200, ],
  col   = "black",
  alpha = 0.12
)

plot_genome_track(
  genome_name = genome_name,
  genome_size = chr_length,
  track_list  = tracks_hl,
  circular    = TRUE
)

Circular plot with per-track highlight bands on the Microsatellites track.

5.8 Same GC content, different Tm

The tracks above show GC content and Tm as two separate rings, which raises the obvious question of whether the second adds anything to the first. It does, and the genome-wide profile is enough to measure how much.

Group the windows by GC content. Within a group every window has the same length, the same salt, and the same GC percentage, so an empirical formula whose only sequence input is GC percentage assigns all of them a single temperature. A nearest-neighbour model reads the dinucleotide stacks and does not.

## The GC-content prediction for the same windows, for comparison.
tm_gc_ecoli <- tm_calculate(genome_name, regions = chr_name, unit = "region",
                            window = 200L, slide = 200L,
                            method = "tm_gc", variant = "Schildkraut1965",
                            Na = 50, verbose = FALSE)$gr
GCmod <- as.data.frame(tm_gc_ecoli[, c("Tm", "GC")])

## Keep only windows whose GC percentage is an exact integer. A 200 bp
## window can only take GC values in steps of 0.5%, so binning 49.5% with
## 50.0% would put a real GC difference inside a group and some of the
## spread below could be attributed to GC rather than to composition.
## The same test also removes windows that contained an ambiguous base,
## whose GC denominator is smaller than the window.
prep <- function(d) {
  d <- d[!is.na(d$Tm) & !is.na(d$GC) & d$width == 200L, ]
  d <- d[d$GC == round(d$GC), ]
  d$GCi <- as.integer(d$GC)
  d
}
Cnn <- prep(as.data.frame(tm_ecoli[, c("Tm", "GC")]))
Cgc <- prep(GCmod)

cnt <- table(Cnn$GCi)
lv  <- sort(as.integer(names(cnt[cnt >= 30])))   # enough windows to be stable
Cnn <- Cnn[Cnn$GCi %in% lv, ]
Cgc <- Cgc[Cgc$GCi %in% lv, ]

modal   <- lv[which.max(cnt[as.character(lv)])]
at_mode <- Cnn$Tm[Cnn$GCi == modal]

cat(sprintf("At %d%% GC: %s windows, Tm from %.2f to %.2f (spread %.2f C)\n",
            modal, format(length(at_mode), big.mark = ","),
            min(at_mode), max(at_mode), diff(range(at_mode))))

## At 53% GC: 927 windows, Tm from 97.92 to 102.52 (spread 4.60 C)

cat(sprintf("GC-content formula predicts a single value: %.2f C\n",
            median(Cgc$Tm[Cgc$GCi == modal])))

## GC-content formula predicts a single value: 78.26 C

cat(sprintf("Largest spread at any GC value: %.2f C\n",
            max(tapply(Cnn$Tm, Cnn$GCi, function(v) diff(range(v))))))

## Largest spread at any GC value: 6.07 C

gc_line <- tapply(Cgc$Tm, Cgc$GCi, median)[as.character(lv)]

op <- par(mar = c(4.6, 4.6, 1.2, 1.2), las = 1, mgp = c(2.9, 0.7, 0))
boxplot(Tm ~ factor(GCi, levels = lv), data = Cnn, at = seq_along(lv),
        outline = FALSE, border = "#34495E", col = "#D6E4F0", lwd = 0.7,
        xaxt = "n", bty = "n", xlab = "GC content (%)",
        ylab = expression(paste("Melting temperature (", degree, "C)")))
sel <- seq(1, length(lv), by = 5)          # one label per box is unreadable
axis(1, at = seq_along(lv)[sel], labels = lv[sel])
lines(seq_along(lv), gc_line, col = "#C0392B", lwd = 2.2)
legend("topleft", bty = "n", cex = 0.85,
       legend = c("Nearest-neighbour (Breslauer 1986)",
                  "GC-content formula (Schildkraut 1965)"),
       col = c("#34495E", "#C0392B"), lwd = c(0.7, 2.2), seg.len = 1.4)

Melting temperature against GC content for 200 bp windows of the E. coli chromosome. Each box holds windows of identical length, salt and GC content. The red line is the GC-content formula, which by construction has no width within a group; the boxes do, and that width is the composition effect the formula cannot represent.

par(op)

Two things to read off this, and one not to.

The boxes have real width. Windows that agree on GC content to the base still differ in Tm, because AA/TT, AT/TA and TA/AT stacks are not interchangeable and neither are the ten GC-containing ones. That difference is the whole content of a nearest-neighbour model, and it is invisible to any formula parameterised on GC percentage alone.

The width is not uniform. It is largest in the middle of the GC range, where the number of distinct sequences compatible with a given GC content is largest, and narrows at both extremes, where composition is increasingly constrained. A GC-content formula is therefore least reliable exactly where most of the genome sits.

What should not be read off it is the vertical offset between the line and the boxes. The two are different parameterisations calibrated against different reference conditions, and comparing their absolute values is a separate question from the one this panel asks. The claim here concerns the spread at fixed GC, which is a property of the sequences and does not depend on where either scale is anchored.

script <- system.file("scripts", "make_figure3.R", package = "TmCalculator")
file.edit(script)      # --locus-kb, --n-groups, --pairs-in-peaks, --dpi

6 Statistical testing

We test whether Tm values inside MutL-AR peaks differ significantly from the rest of the genome.

mutH_peaks <- GRanges(
  seqnames = ecoli_rep_hotspots$all_peaks_IP_mutH$chr,
  ranges   = IRanges(start = ecoli_rep_hotspots$all_peaks_IP_mutH$start,
                     end   = ecoli_rep_hotspots$all_peaks_IP_mutH$end)
)
seqlevels(mutH_peaks) <- "U00096.3"

mutH_peaks$peak_id <- paste0("mutH_", seq_along(mutH_peaks))

tm_annot <- integrate_granges(
  gr_tm          = tm_ecoli,
  gr_features    = mutH_peaks,
  strategy       = "overlap",
  feature_cols   = "peak_id",
  keep_unmatched = TRUE
)

tm_annot$in_mutH <- ifelse(is.na(tm_annot$peak_id), "non_peak", "peak")
table(tm_annot$in_mutH)

## 
## non_peak     peak 
##    22402      806

res <- compare_groups(
  gr          = tm_annot,
  target      = c("Tm", "GC"),
  method      = "wilcoxon",
  group       = "in_mutH",
  alternative = "greater",
  posthoc     = FALSE
)
res$results

##   target   group   method              test statistic df      p.value n_groups
## 1     Tm in_mutH wilcoxon Wilcoxon rank-sum  11650602 NA 4.824689e-45        2
## 2     GC in_mutH wilcoxon Wilcoxon rank-sum  10806568 NA 8.549859e-22        2
##   n_total    group_levels                  group_n
## 1   23208 non_peak | peak non_peak=22402, peak=806
## 2   23208 non_peak | peak non_peak=22402, peak=806

res$summary

##      group     n     mean       sd   median target
## 1 non_peak 22402 98.97778 3.791850 99.65325     Tm
## 2     peak   806 96.96551 4.206166 97.88956     Tm
## 3 non_peak 22402 50.88340 6.346952 52.00000     GC
## 4     peak   806 48.22333 7.394710 50.00000     GC

7 Interpreting the Results

MutL-AR-associated windows have lower Tm and lower GC content than the rest of the genome (p < 0.001 for both; median Tm about 1 degree C lower and mean GC 0.48 vs 0.51 inside peaks), consistent with the enrichment of low-thermal-stability, AT-rich sequence reported by Hasenauer et al. MutL-AR peaks cluster near the replication terminus, where replication forks converge and mismatch density is highest. This spatial coincidence with locally elevated microsatellite density and cruciform-forming sequences is consistent with the replication stress model of MMR recruitment.

GATC methylation sites are distributed genome-wide but show local density fluctuations that partially anti-correlate with GC content, reflecting the sequence context requirements for Dam methyltransferase (5’-GATC-3’).

8 Computational Performance

On a standard desktop computer (single core, R 4.4, compiled Rcpp nearest-neighbor core):

The same call on the human genome, where there are 14.7 million windows rather than 23,208, takes about ten minutes in one process and about three minutes over five workers; see vignette("hg38_performance_parallel").

9 References

Breslauer KJ, Frank R, Blocker H, Marky LA (1986). Predicting DNA duplex stability from the base sequence. Proceedings of the National Academy of Sciences USA 83(11): 3746-3750. doi:10.1073/pnas.83.11.3746

Hasenauer FC, Barreto HC, Lotton C, Matic I (2025). Genome-wide mapping of spontaneous DNA replication error-hotspots using mismatch repair proteins in rapidly proliferating Escherichia coli. Nucleic Acids Research 53(2): gkae1196. doi:10.1093/nar/gkae1196

10 Session Information

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sonoma 14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] BSgenome.Ecoli.NCBI.ASM584v2_1.0.0 BSgenome_1.72.0                   
##  [3] rtracklayer_1.64.0                 BiocIO_1.14.0                     
##  [5] Biostrings_2.72.1                  XVector_0.44.0                    
##  [7] GenomicRanges_1.56.2               GenomeInfoDb_1.40.1               
##  [9] IRanges_2.38.1                     S4Vectors_0.42.1                  
## [11] BiocGenerics_0.50.0                TmCalculator_1.1.1                
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.3                   bitops_1.0-9               
##   [3] gridExtra_2.3               rlang_1.1.7                
##   [5] magrittr_2.0.4              biovizBase_1.52.0          
##   [7] otel_0.2.0                  matrixStats_1.5.0          
##   [9] compiler_4.4.1              RSQLite_2.4.0              
##  [11] GenomicFeatures_1.56.0      png_0.1-8                  
##  [13] vctrs_0.7.1                 ProtGenerics_1.36.0        
##  [15] stringr_1.6.0               pkgconfig_2.0.3            
##  [17] crayon_1.5.3                fastmap_1.2.0              
##  [19] backports_1.5.0             Rsamtools_2.20.0           
##  [21] rmarkdown_2.30              UCSC.utils_1.0.0           
##  [23] bit_4.6.0                   xfun_0.58                  
##  [25] zlibbioc_1.50.0             cachem_1.1.0               
##  [27] jsonlite_2.0.0              blob_1.2.4                 
##  [29] DelayedArray_0.30.1         BiocParallel_1.38.0        
##  [31] parallel_4.4.1              cluster_2.1.6              
##  [33] R6_2.6.1                    VariantAnnotation_1.50.0   
##  [35] stringi_1.8.7               bslib_0.10.0               
##  [37] RColorBrewer_1.1-3          bezier_1.1.2               
##  [39] rpart_4.1.23                jquerylib_0.1.4            
##  [41] Rcpp_1.1.2                  SummarizedExperiment_1.34.0
##  [43] knitr_1.51                  base64enc_0.1-6            
##  [45] Matrix_1.7-0                nnet_7.3-19                
##  [47] tidyselect_1.2.1            rstudioapi_0.18.0          
##  [49] dichromat_2.0-0.1           abind_1.4-8                
##  [51] yaml_2.3.12                 codetools_0.2-20           
##  [53] curl_6.2.3                  lattice_0.22-6             
##  [55] tibble_3.2.1                regioneR_1.36.0            
##  [57] Biobase_2.64.0              KEGGREST_1.44.1            
##  [59] evaluate_1.0.5              foreign_0.8-87             
##  [61] karyoploteR_1.30.0          pillar_1.10.2              
##  [63] MatrixGenerics_1.16.0       checkmate_2.3.2            
##  [65] generics_0.1.4              RCurl_1.98-1.17            
##  [67] ensembldb_2.28.1            ggplot2_3.5.2              
##  [69] scales_1.4.0                glue_1.8.0                 
##  [71] lazyeval_0.2.2              Hmisc_5.2-3                
##  [73] tools_4.4.1                 data.table_1.17.4          
##  [75] GenomicAlignments_1.40.0    XML_3.99-0.18              
##  [77] grid_4.4.1                  colorspace_2.1-1           
##  [79] AnnotationDbi_1.66.0        GenomeInfoDbData_1.2.12    
##  [81] htmlTable_2.4.3             restfulr_0.0.15            
##  [83] Formula_1.2-5               cli_3.6.5                  
##  [85] S4Arrays_1.4.1              dplyr_1.1.4                
##  [87] AnnotationFilter_1.28.0     gtable_0.3.6               
##  [89] sass_0.4.10                 digest_0.6.39              
##  [91] SparseArray_1.4.8           rjson_0.2.23               
##  [93] htmlwidgets_1.6.4           farver_2.1.2               
##  [95] memoise_2.0.1               htmltools_0.5.9            
##  [97] lifecycle_1.0.5             httr_1.4.7                 
##  [99] bit64_4.6.0-1               bamsignals_1.36.0