--- title: "Genome-wide Melting Temperature Profiling: An E. coli Case Study" author: "Junhui Li, Lihua Julie Zhu" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true vignette: > %\VignetteIndexEntry{Genome-wide Melting Temperature Profiling: An E. coli Case Study} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", # html_vignette inherits fig.retina = 2 from rmarkdown, which renders every # figure at twice the nominal resolution and then scales it down in the # browser. The pixel count, and therefore the size of the base64 blob # embedded in the self-contained HTML, is four times larger for no visible # gain in a vignette. Setting it to 1 is what keeps the installed size of # doc/ within the limit R CMD check reports on. fig.retina = 1, dpi = 72, fig.width = 7, fig.height = 6 ) ``` # 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()`. --- # Prerequisites ## R packages ```{r packages} library(TmCalculator) library(BSgenome) library(GenomicRanges) ``` # The reference genome {#build-bsgenome} `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. ```{r available-genomes, eval=FALSE} 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. ```{r forge, eval=FALSE} library(BSgenomeForge) forgeBSgenomeDataPkgFromNCBI( assembly_accession = "GCF_000005845.2", pkg_maintainer = "Junhui Li ", 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. ```{r setup-ecoli-bsgenome, message=FALSE, warning=FALSE} 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 ", 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): ```{r load-genome} 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") cat("Length: ", format(chr_length, big.mark = ","), "bp\n") ``` --- # 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. ```{r tm-profile} 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 <- as.data.frame(tm_ecoli[, c("Tm", "GC")]) summary(Tm[, c("Tm", "GC")]) ``` The result is an ordinary `GRanges`, so the coordinates, the melting temperatures and the GC percentages travel together through everything that follows. ```{r profile-head} head(tm_ecoli, 3) ``` 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: ```{r tm-calculate-seqs} tm_calculate(c("ACGTGCTAGCTAGCTAGC", "GGCCATATATGCGC"), method = "tm_nn", Na = 50) ``` --- # 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. ## Assemble the track list We overlay the Tm/GC profile with the Hasenauer *et al.* (2025) multi-omics layers from `ecoli_rep_hotspots`: | Track | Description | |-------|-------------| | MutL-AR | MutL ChIP-seq peaks marking replication-error / mismatch-repair hotspots | | Microsatellites | Tandem-repeat (mononucleotide-repeat) density per 1 kb bin | | Cruciform | Cruciform-forming sequence density per 1 kb bin | | ssDNA | Single-stranded DNA regions enriched near error hotspots | | GATC sites | Dam methylation-site (5'-GATC-3') density per 1 kb bin | ```{r track-list} # 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) ) ``` ## Circular genome map ```{r circular-full, fig.width=8, fig.height=8, fig.cap="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."} plot_genome_track( genome_name = genome_name, genome_size = chr_length, track_list = tracks, circular = TRUE, label = label ) ``` ## 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. ```{r linear-full, fig.width=10, fig.height=5, fig.cap="Linear genome view of E. coli K-12 MG1655. MutL-AR peaks are drawn inside the chromosome ideogram bar."} plot_genome_track( genome_name = genome_name, genome_size = chr_length, track_list = tracks ) ``` ## 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. ```{r zoom-single, fig.width=10, fig.height=8, fig.cap="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."} 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 ) ``` ## 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. ```{r zoom-multi-circular, fig.width=10, fig.height=9, fig.cap="Two zoomed regions, 0.1-0.5 Mb and 3.6-4.5 Mb, concatenated around the circle with a gap 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") ) ``` ## 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. ```{r circular-pan, fig.width=7, fig.height=7, fig.cap="Panned circular view showing the upper-right quadrant of the E. coli chromosome."} 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) ) ``` ## 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: ```{r per-track-highlight, fig.width=8, fig.height=8, fig.cap="Circular plot with per-track highlight bands on the Microsatellites track."} 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 ) ``` ## 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. ```{r gc-vs-tm-data} ## 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)))) cat(sprintf("GC-content formula predicts a single value: %.2f C\n", median(Cgc$Tm[Cgc$GCi == modal]))) cat(sprintf("Largest spread at any GC value: %.2f C\n", max(tapply(Cnn$Tm, Cnn$GCi, function(v) diff(range(v)))))) ``` ```{r gc-vs-tm-plot, fig.width=7, fig.height=5, fig.cap="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."} 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) 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. ## One locus, for the manuscript figure The distribution above says the effect holds across the genome. Seeing it in a single place takes a search: score every 20 kb locus by how far apart its equal-GC windows are in Tm, and take the best one. That search, and the three-panel figure it feeds (Figure 3 of the manuscript), are in `inst/scripts/make_figure3.R`, which is shipped with the package: ```{r make-figure3, eval=FALSE} script <- system.file("scripts", "make_figure3.R", package = "TmCalculator") file.edit(script) # --locus-kb, --n-groups, --pairs-in-peaks, --dpi ``` The panel it produces shows a flat GC track across each shaded pair of windows and a Tm track that is not, which is the argument for a nearest-neighbour model in one picture. Two caveats belong with it: the locus is chosen *because* it is extreme, so the box plot above is the honest summary of typical behaviour; and the pairs are extreme within the locus but need not both fall inside a peak, which `--pairs-in-peaks` enforces at the cost of a smaller effect. --- # Statistical testing We test whether Tm values inside MutL-AR peaks differ significantly from the rest of the genome. ## Build MutL-AR peak GRanges ```{r build-mutH-peaks} 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" ``` ## Annotate Tm tiles with peak membership ```{r annotate-mutH-peaks} 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) ``` ## Wilcoxon test on Tm and GC ```{r wilcoxon} res <- compare_groups( gr = tm_annot, target = c("Tm", "GC"), method = "wilcoxon", group = "in_mutH", alternative = "greater", posthoc = FALSE ) res$results res$summary ``` --- # 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'). --- # Computational Performance On a standard desktop computer (single core, R 4.4, compiled Rcpp nearest-neighbor core): | Step | Function | Time (elapsed) | |------|----------|---------------| | Tiling, sequence retrieval and Tm for 23,208 windows | `tm_calculate()` | ~`r sprintf("%.2f", runtime[["elapsed"]])` s | 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")`. --- # 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. 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. --- # Session Information ```{r session-info} sessionInfo() ```