Sensitivity analysis
For each window size we retile the chromosome, recompute Tm and GC,
and re-run the MutL-AR comparison, keeping every other setting fixed.
Nothing downstream needs to change: the annotation layers keep their own
native resolutions (1 kb bins for the microsatellite, cruciform and GATC
density tracks; variable-width intervals for MutL-AR peaks and ssDNA
regions), and both plot_genome_track() and the
overlap-based statistics work across mixed resolutions.
window_sizes <- c(50L, 100L, 200L, 500L)
sens <- lapply(window_sizes, function(w) {
bins_w <- make_genomiccoord(
bsgenome = genome_name, chromosomes = chr_name,
window = w, slide = w, start = 1, end = chr_length,
strand = "+", verbose = FALSE
)
gr_w <- to_genomic_ranges_fast(list(pkg_name = genome_name, seq = bins_w))
tm_w <- tm_calculate(gr_w, method = "tm_nn",
nn_table = "DNA_NN_Breslauer_1986", Na = 50)$gr
ann <- integrate_granges(gr_tm = tm_w, gr_features = mutH_peaks,
strategy = "overlap", feature_cols = "peak_id",
keep_unmatched = TRUE)
ann$in_mutH <- ifelse(is.na(ann$peak_id), "non_peak", "peak")
cg <- compare_groups(gr = ann, target = c("Tm", "GC"),
method = "wilcoxon", group = "in_mutH",
alternative = "greater", posthoc = FALSE)
list(gr = tm_w, ann = ann, test = cg)
})
names(sens) <- paste0("w", window_sizes)Distributional stability
dist_tbl <- do.call(rbind, lapply(seq_along(window_sizes), function(i) {
g <- sens[[i]]$gr
data.frame(
window_bp = window_sizes[i],
n_windows = length(g),
Tm_mean = mean(g$Tm, na.rm = TRUE),
Tm_sd = stats::sd(g$Tm, na.rm = TRUE),
Tm_IQR = stats::IQR(g$Tm, na.rm = TRUE),
GC_mean = mean(g$GC, na.rm = TRUE)
)
}))
knitr::kable(dist_tbl, digits = 3,
caption = "Tm/GC distribution by window size.")| window_bp | n_windows | Tm_mean | Tm_sd | Tm_IQR | GC_mean |
|---|---|---|---|---|---|
| 50 | 92833 | 83.891 | 5.389 | 6.877 | 50.791 |
| 100 | 46416 | 93.818 | 4.423 | 5.320 | 50.791 |
| 200 | 23208 | 98.908 | 3.825 | 4.422 | 50.791 |
| 500 | 9283 | 101.998 | 3.251 | 3.628 | 50.791 |
The same two effects are easier to read as distributions than as summary statistics: the curve moves to the right as the window lengthens, and it narrows at the same time.
## Palettes are defined here because they are used both by this panel and by
## the multi-scale track figures further down.
## darkest = finest window (50 bp), lightest = coarsest (500 bp)
scale_cols_gc <- c("#1A5276", "#2E86C1", "#5DADE2", "#AED6F1")
scale_cols_tm <- c("#7B241C", "#CB4335", "#EC7063", "#F5B7B1")
## Densities are normalised, so the ~93,000 windows at 50 bp and the ~9,300
## at 500 bp can be compared directly on one pair of axes.
tm_by_size <- lapply(sens, function(x) x$gr$Tm[is.finite(x$gr$Tm)])
dens <- lapply(tm_by_size, stats::density)
xlim <- range(vapply(dens, function(z) range(z$x), numeric(2)))
ylim <- c(0, max(vapply(dens, function(z) max(z$y), numeric(1))))
op <- par(mar = c(4.2, 4.4, 0.8, 0.8), las = 1)
plot(NA, xlim = xlim, ylim = ylim, bty = "n",
xlab = expression(italic(T)[m] ~ "(" * degree * "C)"),
ylab = "Density")
for (i in seq_along(dens)) {
lines(dens[[i]], col = scale_cols_tm[i], lwd = 2)
abline(v = dist_tbl$Tm_mean[i], col = scale_cols_tm[i], lwd = 1, lty = 3)
}
legend("topleft", bty = "n", lwd = 2, seg.len = 1.6,
col = scale_cols_tm[seq_along(window_sizes)],
legend = sprintf("%d bp: mean %.1f, SD %.2f",
dist_tbl$window_bp, dist_tbl$Tm_mean, dist_tbl$Tm_sd))
Tm distribution at each window size. The distribution shifts to higher Tm as the window lengthens and narrows at the same time; dotted vertical lines mark the group means. Colours match the multi-scale tracks below (darkest = 50 bp, lightest = 500 bp).
par(op)
## Spread of the finest window relative to the coarsest. This is the ratio
## quoted in the manuscript, so both come from the same computation.
sd_ratio <- dist_tbl$Tm_sd[1] / dist_tbl$Tm_sd[nrow(dist_tbl)]Two systematic effects are visible, both expected. First, mean Tm rises with window length (83.9 to 102.0 degrees C from 50 to 500 bp): duplex stability increases with length in the nearest-neighbor model, so absolute Tm values are only comparable within one window size. Second, shorter windows widen the distribution (SD 5.39 at 50 bp vs 3.25 at 500 bp, a ratio of 1.66) – window size acts as a smoothing bandwidth. The mean GC is invariant (50.791 throughout), confirming that the composition landscape itself does not depend on the tiling.
Cross-scale agreement
Every profile is aggregated onto one common grid (mean Tm per bin) and the resulting series are correlated.
The grid is 1 kb rather than 500 bp, and the choice matters. Correlating against the 500 bp profile directly makes 200 bp the odd one out: 50 and 100 divide 500 exactly, into ten and five windows per bin, but 200 does not, so each bin receives two or three windows that straddle its boundaries and the aggregation is itself inexact. That alone drove the 200 bp correlation down to 0.976 while the two finer profiles reached 0.997 and 0.999, an ordering that would invite the reader to conclude that the closest resolution agrees worst. All four window sizes divide 1 kb (into 20, 10, 5 and 2), so on that grid every profile is aggregated exactly and the comparison is between profiles rather than between remainders. The 500 bp profile is aggregated too, and serves as the reference.
grid_bp <- 1000L
stopifnot(all(grid_bp %% window_sizes == 0L)) # exact aggregation for each
gr_ref <- sens[["w500"]]$gr
ref_agg <- tapply(gr_ref$Tm, (start(gr_ref) - 1L) %/% grid_bp,
mean, na.rm = TRUE)
cor_tbl <- vapply(c("w50", "w100", "w200"), function(k) {
g <- sens[[k]]$gr
agg <- tapply(g$Tm, (start(g) - 1L) %/% grid_bp, mean, na.rm = TRUE)
common <- intersect(names(agg), names(ref_agg))
stats::cor(agg[common], ref_agg[common], use = "complete.obs")
}, numeric(1))
round(cor_tbl, 3)## w50 w100 w200
## 0.998 0.999 1.000Aggregated onto the common 1 kb grid, the profiles are nearly interchangeable (r = 0.998, 0.999, 1.000 for 50, 100 and 200 bp respectively). Agreement now rises monotonically as the window approaches the grid, which is what averaging fewer windows per bin should do, and is the behaviour the 500 bp grid obscured. Window size rescales and smooths the profile but preserves the spatial landscape.
Robustness of the MutL-AR association
Finally, the vignette’s biological conclusion – that Tm and GC differ between MutL-AR peak windows and the genomic background – is re-tested at every window size:
## Per-window-size test statistics
test_results <- do.call(rbind, lapply(seq_along(window_sizes), function(i) {
data.frame(window_bp = window_sizes[i], sens[[i]]$test$results)
}))
test_results## window_bp target group method test statistic df p.value
## 1 50 Tm in_mutH wilcoxon Wilcoxon rank-sum 168746210 NA 2.432173e-90
## 2 50 GC in_mutH wilcoxon Wilcoxon rank-sum 160420084 NA 7.452204e-48
## 3 100 Tm in_mutH wilcoxon Wilcoxon rank-sum 44130350 NA 1.530886e-65
## 4 100 GC in_mutH wilcoxon Wilcoxon rank-sum 41400570 NA 1.192782e-32
## 5 200 Tm in_mutH wilcoxon Wilcoxon rank-sum 11650602 NA 4.824689e-45
## 6 200 GC in_mutH wilcoxon Wilcoxon rank-sum 10806568 NA 8.549859e-22
## 7 500 Tm in_mutH wilcoxon Wilcoxon rank-sum 2033341 NA 4.905312e-25
## 8 500 GC in_mutH wilcoxon Wilcoxon rank-sum 1873428 NA 1.417503e-12
## n_groups n_total group_levels group_n
## 1 2 92833 non_peak | peak non_peak=89730, peak=3103
## 2 2 92833 non_peak | peak non_peak=89730, peak=3103
## 3 2 46416 non_peak | peak non_peak=44845, peak=1571
## 4 2 46416 non_peak | peak non_peak=44845, peak=1571
## 5 2 23208 non_peak | peak non_peak=22402, peak=806
## 6 2 23208 non_peak | peak non_peak=22402, peak=806
## 7 2 9283 non_peak | peak non_peak=8940, peak=343
## 8 2 9283 non_peak | peak non_peak=8940, peak=343## Per-window-size group summaries, plus the median Tm effect size
test_summary <- do.call(rbind, lapply(seq_along(window_sizes), function(i) {
ann <- sens[[i]]$ann
med <- tapply(ann$Tm, ann$in_mutH, stats::median, na.rm = TRUE)
data.frame(window_bp = window_sizes[i],
n_peak = sum(ann$in_mutH == "peak"),
Tm_med_diff = unname(med["peak"] - med["non_peak"]),
sens[[i]]$test$summary)
}))
test_summary## window_bp n_peak Tm_med_diff group n mean sd median
## 1 50 3103 -2.003979 non_peak 89730 83.95988 5.365263 84.44240
## 2 50 3103 -2.003979 peak 3103 81.88670 5.694383 82.43842
## 3 50 3103 -2.003979 non_peak 89730 50.87955 8.832840 52.00000
## 4 50 3103 -2.003979 peak 3103 48.22172 9.789787 50.00000
## 5 100 1571 -1.663404 non_peak 44845 93.88741 4.393684 94.51009
## 6 100 1571 -1.663404 peak 1571 91.85059 4.780258 92.84668
## 7 100 1571 -1.663404 non_peak 44845 50.88135 7.326739 52.00000
## 8 100 1571 -1.663404 peak 1571 48.21260 8.384205 50.00000
## 9 200 806 -1.763693 non_peak 22402 98.97778 3.791850 99.65325
## 10 200 806 -1.763693 peak 806 96.96551 4.206166 97.88956
## 11 200 806 -1.763693 non_peak 22402 50.88340 6.346952 52.00000
## 12 200 806 -1.763693 peak 806 48.22333 7.394710 50.00000
## 13 500 343 -1.703604 non_peak 8940 102.06820 3.213329 102.71352
## 14 500 343 -1.703604 peak 343 100.16319 3.679005 101.00992
## 15 500 343 -1.703604 non_peak 8940 50.88456 5.386122 51.80000
## 16 500 343 -1.703604 peak 343 48.35219 6.463921 50.20000
## target
## 1 Tm
## 2 Tm
## 3 GC
## 4 GC
## 5 Tm
## 6 Tm
## 7 GC
## 8 GC
## 9 Tm
## 10 Tm
## 11 GC
## 12 GC
## 13 Tm
## 14 Tm
## 15 GC
## 16 GCOne table for the whole analysis
Both questions are answered by the same four rows, so they are assembled into a single table: the first four columns say what window size does to the landscape, and the last three say what it does to the conclusion drawn from it.
sens_tbl <- do.call(rbind, lapply(seq_along(window_sizes), function(i) {
g <- sens[[i]]$gr
ann <- sens[[i]]$ann
med <- tapply(ann$Tm, ann$in_mutH, stats::median, na.rm = TRUE)
res <- sens[[i]]$test$results
w <- window_sizes[i]
data.frame(
`Window (bp)` = w,
Windows = length(g),
`Tm mean (C)` = mean(g$Tm, na.rm = TRUE),
`Tm SD` = stats::sd(g$Tm, na.rm = TRUE),
`GC mean (%)` = mean(g$GC, na.rm = TRUE),
## The reference profile correlates with itself by construction.
`r vs 1 kb grid` = if (w == max(window_sizes)) NA_real_
else unname(cor_tbl[[paste0("w", w)]]),
`MutL-AR windows` = sum(ann$in_mutH == "peak"),
`Tm peak - bg (C)` = unname(med["peak"] - med["non_peak"]),
## Formatted as character: these p values span sixty-five orders of
## magnitude, and rounding them to a fixed number of decimals prints
## every one of them as zero.
`p (Tm)` = format(res$p.value[res$target == "Tm"],
digits = 2, scientific = TRUE),
check.names = FALSE, stringsAsFactors = FALSE)
}))
knitr::kable(sens_tbl, digits = c(0, 0, 2, 2, 2, 3, 0, 3, 0),
caption = paste("Window-size sensitivity. GC mean is invariant;",
"the Tm distribution shifts and narrows; the",
"MutL-AR effect size varies within a few tenths",
"of a degree while the p value tracks the number",
"of windows."))| Window (bp) | Windows | Tm mean (C) | Tm SD | GC mean (%) | r vs 1 kb grid | MutL-AR windows | Tm peak - bg (C) | p (Tm) |
|---|---|---|---|---|---|---|---|---|
| 50 | 92833 | 83.89 | 5.39 | 50.79 | 0.998 | 3103 | -2.004 | 2.4e-90 |
| 100 | 46416 | 93.82 | 4.42 | 50.79 | 0.999 | 1571 | -1.663 | 1.5e-65 |
| 200 | 23208 | 98.91 | 3.82 | 50.79 | 1.000 | 806 | -1.764 | 4.8e-45 |
| 500 | 9283 | 102.00 | 3.25 | 50.79 | NA | 343 | -1.704 | 4.9e-25 |
The last two columns are the point of the section. The p value falls by sixty-five orders of magnitude between the coarsest and the finest tiling, which is a statement about how many windows were tested rather than about how far apart the two groups are. The difference in medians, which is that statement, stays between -1.7 and -2.0 degrees C. The finest tiling gives the largest separation, as expected: a 500 bp window overlapping a MutL-AR peak also contains flanking sequence that does not, and averaging over it dilutes the contrast, whereas a 50 bp window resolves the peak more sharply. The three coarser tilings agree within 0.1 degrees C.
Multi-scale tracks alongside the multi-omics layers
Because every track keeps its native resolution, the four Tm/GC profiles can be drawn as concentric layers of one integrated figure, together with the annotation data:
## scale_cols_gc / scale_cols_tm are defined in the density-panel chunk above
sens_dfs <- lapply(sens, function(x) as.data.frame(x$gr[, c("Tm", "GC")]))
tracks_scale <- c(
list(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)),
lapply(seq_along(window_sizes), function(i)
list(type = "line", data = sens_dfs[[i]], value_col = "GC",
name = paste0("GC ", window_sizes[i], " bp"),
col = scale_cols_gc[i], legend_font_col = scale_cols_gc[i])),
lapply(seq_along(window_sizes), function(i)
list(type = "line", data = sens_dfs[[i]], value_col = "Tm",
name = paste0("Tm ", window_sizes[i], " bp"),
col = scale_cols_tm[i], legend_font_col = scale_cols_tm[i])),
list(list(type = "line", data = ecoli_rep_hotspots$bins_rep,
value_col = "count", name = "Microsatellites", col = "#2ECC71",
legend_font_col = "#2ECC71"),
list(type = "highlight",
data = ecoli_rep_hotspots$all_peaks_IP_mutH,
col = "#F1C40F", alpha = 0.18))
)
plot_genome_track(
genome_name = genome_name,
genome_size = chr_length,
track_list = tracks_scale,
circular = TRUE,
label = label
)
Multi-scale integration. Concentric rings from outside in: MutL-AR peaks (ideogram), GC content at 50/100/200/500 bp (blues), Tm at 50/100/200/500 bp (reds), and microsatellite density (green). Yellow bands mark MutL-AR peaks.
At whole-genome scale the four bandwidths look nearly identical – the landscape is invariant. The difference window size makes is local smoothing, best seen in a zoomed view of the same interval used by the main case study:
## Only the four Tm profiles are drawn. Ten line tracks in one linear panel
## leave each too little vertical space for its own axis labels, which then
## collide, and everything except Tm is answering a different question: the
## GC layers duplicate what the summary table already gives numerically, and
## the microsatellite density is the subject of the main case study's figure,
## not of this one. Removing both leaves five tracks with room to be read.
## The full ten-track version is the circular figure above.
##
## The interval is the one used in the main case study's zoomed panel, so the
## two figures can be read against each other rather than against different
## parts of the chromosome.
tracks_zoom <- c(
list(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.6)),
lapply(seq_along(window_sizes), function(i)
list(type = "line", data = sens_dfs[[i]], value_col = "Tm",
name = paste0("Tm ", window_sizes[i], " bp"),
col = scale_cols_tm[i], legend_font_col = scale_cols_tm[i],
height = 1.2)),
list(list(type = "highlight",
data = ecoli_rep_hotspots$all_peaks_IP_mutH,
col = "#F1C40F", alpha = 0.18))
)
plot_genome_track(
genome_name = genome_name,
genome_size = chr_length,
track_list = tracks_zoom,
zoom = "U00096.3:100000-300000",
track.gap = 0.03,
axis.cex = 0.55
)
The 0.1-0.3 Mb region, within the interval shown in the genome-wide figure of the main case study, with Tm computed at 50, 100, 200 and 500 bp (dark to light). The 50 bp profile resolves single-window fluctuations that the 500 bp profile averages away, but all four trace the same underlying landscape. Yellow bands mark MutL-AR peaks.
The association is reproduced at every scale. MutL-AR windows are lower in median Tm than the background at all four window sizes (-2.00, -1.66, -1.76, -1.70 degrees C at 50, 100, 200, 500 bp) and lower in GC (48.252 in peaks vs 50.882 in the background), and every Wilcoxon test is significant (largest p = 1.4e-12). As expected, p values shrink with the number of windows while the effect size stays essentially constant, so the biological conclusion does not depend on the window choice. The primary analysis uses 200 bp to match Hasenauer et al. (2025), whose Methods compute GC content and melting temperature in 200 bp bins (using TmCalculator) and bin the ChIP-seq coverage at the same 200 bp resolution.