diff --git a/DESCRIPTION b/DESCRIPTION
index c38a3923..000d8749 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: TwoSampleMR
 Title: Two Sample MR Functions and Interface to MRC Integrative
     Epidemiology Unit OpenGWAS Database
-Version: 0.7.6
+Version: 0.7.7
 Authors@R: c(
     person("Gibran", "Hemani", , "g.hemani@bristol.ac.uk", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-0920-1055")),
@@ -29,6 +29,7 @@ BugReports: https://github.com/MRCIEU/TwoSampleMR/issues/
 Depends:
     R (>= 4.1.0)
 Imports:
+    cli,
     cowplot,
     data.table,
     dplyr,
@@ -47,8 +48,8 @@ Imports:
     pbapply,
     psych,
     RadialMR,
-    reshape2,
-    rmarkdown
+    rmarkdown,
+    tidyr
 Suggests:
     Cairo,
     car,
@@ -58,8 +59,7 @@ Suggests:
     mr.raps,
     MRInstruments,
     randomForest,
-    testthat,
-    tidyr
+    testthat
 VignetteBuilder:
     knitr
 Remotes:
diff --git a/NEWS.md b/NEWS.md
index 9cc5f630..100dbb4c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,20 @@
+# TwoSampleMR v0.7.7
+
+(Release date 2026-06-07)
+
+* Remove dead `exposure_mat` assignment
+* Replace `reshape2::dcast()` with `tidyr::pivot_wider()`
+* Replace deprecated `ggplot2::aes_string()` with `ggplot2::aes()`
+* Add new bare names to `globalVariables()`
+* Replace deprecated `size` with `linewidth` in `geom_vline()`, `geom_errorbarh()`, `geom_errorbar()`, and `element_rect()`
+* Fix constant-in-`aes`, redundant `sapply`, and `paste`/`paste0`
+* Qualify `predict()` and `desc()` calls; remove from `globalVariables()`
+* Remove broken and unreachable `mr_mode_broken()` function and its `globalVariables()` entry
+* Replace `sapply` with `invisible(lapply())` for side-effect iteration; use `paste0()` for single leading space
+* Replace `1:n` index patterns with `seq_len()` in chunking logic
+* Hoist `paste()` out of inner loop; replace nested loops with single vectorised assignment
+* Make startup message URLs clickable via `cli::style_hyperlink()`; add cli to Imports
+
 # TwoSampleMR v0.7.6
 
 (Release date 2026-05-15)
diff --git a/R/forest_plot.R b/R/forest_plot.R
index 6860969e..2523d76d 100644
--- a/R/forest_plot.R
+++ b/R/forest_plot.R
@@ -194,7 +194,7 @@ mr_forest_plot_grouped <-
         )
       ) +
         ggplot2::geom_segment() +
-        ggplot2::geom_point(ggplot2::aes(y = space_col, x = as.numeric(eff_col), size = 4)) +
+        ggplot2::geom_point(ggplot2::aes(y = space_col, x = as.numeric(eff_col)), size = 4) +
         ggplot2::theme_bw() +
         ggplot2::theme(
           axis.text.y = ggplot2::element_blank(),
@@ -223,7 +223,7 @@ mr_forest_plot_grouped <-
       data_Fm$text_col <- data_Fm[, text_col]
 
       # A hard rule to set the width of the annotation column, which sometimes truncates very wide columns (complex disease names, numbers with 16 digits, etc)
-      text_widths <- c(-1, max(10, 0.5 * max(sapply(as.character(data_Fm[, text_col]), nchar))))
+      text_widths <- c(-1, max(10, 0.5 * max(nchar(as.character(data_Fm[, text_col])))))
 
       # GGplot rendering of the annotation column
       lefttext <- ggplot2::ggplot(
@@ -296,7 +296,7 @@ mr_forest_plot_grouped <-
       left_Grobs$relative_widths <- NULL
 
       for (i in seq_along(left_Grobs)) {
-        grob_Bag[paste('l', names(left_Grobs)[i], sep = '')] <- left_Grobs[i]
+        grob_Bag[paste0('l', names(left_Grobs)[i])] <- left_Grobs[i]
       }
 
       grob_Bag$m_forest <- ggplotGrob(forst_Pt)
@@ -305,7 +305,7 @@ mr_forest_plot_grouped <-
       right_Grobs <- right_Hs
       right_Grobs$relative_widths <- NULL
       for (i in seq_along(right_Grobs)) {
-        grob_Bag[paste('r', names(right_Grobs)[i], sep = '')] <- right_Grobs[i]
+        grob_Bag[paste0('r', names(right_Grobs)[i])] <- right_Grobs[i]
       }
 
       width_vec <- c(left_RW, 0.34, right_RW)
diff --git a/R/forest_plot2.R b/R/forest_plot2.R
index 31e81156..90200c00 100644
--- a/R/forest_plot2.R
+++ b/R/forest_plot2.R
@@ -325,13 +325,13 @@ forest_plot_basic <- function(
       ggplot2::geom_vline(
         xintercept = seq(ceiling(lo_orig), ceiling(up), by = 0.5),
         colour = "white",
-        size = 0.3
+        linewidth = 0.3
       ) +
-      ggplot2::geom_vline(xintercept = null_line, colour = "#333333", size = 0.3) +
+      ggplot2::geom_vline(xintercept = null_line, colour = "#333333", linewidth = 0.3) +
       ggplot2::geom_errorbarh(
         ggplot2::aes(xmin = lo_ci, xmax = up_ci),
         height = 0,
-        size = 0.4,
+        linewidth = 0.4,
         colour = "#aaaaaa"
       ) +
       ggplot2::geom_point(colour = "black", size = 2.2) +
@@ -358,7 +358,7 @@ forest_plot_basic <- function(
         plot.margin = ggplot2::unit(c(2, 3, 2, 0), units = "points"),
         plot.background = ggplot2::element_rect(fill = "white"),
         panel.spacing = ggplot2::unit(0, "lines"),
-        panel.background = ggplot2::element_rect(colour = "red", fill = "grey", size = 1),
+        panel.background = ggplot2::element_rect(colour = "red", fill = "grey", linewidth = 1),
         strip.text.y = ggplot2::element_blank()
         # strip.background = ggplot2::element_blank()
       ) +
@@ -376,13 +376,13 @@ forest_plot_basic <- function(
       ggplot2::geom_vline(
         xintercept = seq(ceiling(lo_orig), ceiling(up), by = 0.5),
         colour = "white",
-        size = 0.3
+        linewidth = 0.3
       ) +
-      ggplot2::geom_vline(xintercept = null_line, colour = "#333333", size = 0.3) +
+      ggplot2::geom_vline(xintercept = null_line, colour = "#333333", linewidth = 0.3) +
       ggplot2::geom_errorbar(
         ggplot2::aes(xmin = lo_ci, xmax = up_ci),
         width = 0,
-        size = 0.4,
+        linewidth = 0.4,
         colour = "#aaaaaa",
         orientation = "y"
       ) +
@@ -410,7 +410,7 @@ forest_plot_basic <- function(
         plot.margin = ggplot2::unit(c(2, 3, 2, 0), units = "points"),
         plot.background = ggplot2::element_rect(fill = "white"),
         panel.spacing = ggplot2::unit(0, "lines"),
-        panel.background = ggplot2::element_rect(colour = "red", fill = "grey", size = 1),
+        panel.background = ggplot2::element_rect(colour = "red", fill = "grey", linewidth = 1),
         strip.text.y = ggplot2::element_blank()
         # strip.background = ggplot2::element_blank()
       ) +
@@ -499,7 +499,7 @@ forest_plot_names <- function(dat, section = NULL, bottom = TRUE) {
       plot.margin = ggplot2::unit(c(2, 0, 2, 0), units = "points"),
       plot.background = ggplot2::element_rect(fill = "white"),
       panel.spacing = ggplot2::unit(0, "lines"),
-      panel.background = ggplot2::element_rect(colour = "red", fill = "grey", size = 1),
+      panel.background = ggplot2::element_rect(colour = "red", fill = "grey", linewidth = 1),
       strip.text.y = ggplot2::element_blank()
       # strip.background = ggplot2::element_blank()
     ) +
diff --git a/R/format_mr_results2.R b/R/format_mr_results2.R
index 1a1ad782..683e0457 100644
--- a/R/format_mr_results2.R
+++ b/R/format_mr_results2.R
@@ -200,14 +200,11 @@ combine_all_mrresults <- function(
   res$nsnp[is.na(res$nsnp)] <- 1
 
   for (i in unique(res$id.outcome)) {
-    Methods <- unique(res$Method[res$id.outcome == i])
-    Methods <- Methods[Methods != "Wald ratio"]
-    for (j in unique(Methods)) {
-      res$SNP[res$id.outcome == i & res$Method == j] <- paste(
-        res$SNP[res$id.outcome == i & res$Method == "Wald ratio"],
-        collapse = "; "
-      )
-    }
+    wald_snps <- paste(
+      res$SNP[res$id.outcome == i & res$Method == "Wald ratio"],
+      collapse = "; "
+    )
+    res$SNP[res$id.outcome == i & res$Method != "Wald ratio"] <- wald_snps
   }
 
   if (Exp) {
diff --git a/R/globals.R b/R/globals.R
index 7bb3efbc..33b7d889 100644
--- a/R/globals.R
+++ b/R/globals.R
@@ -1,6 +1,9 @@
 utils::globalVariables(c(
+  "Estimate",
+  "Method",
   "MOE",
   "Q",
+  "Qdash",
   "Q_df",
   "Q_pval",
   "Qj",
@@ -15,7 +18,6 @@ utils::globalVariables(c(
   "beta.outcome",
   "bind_rows",
   "category",
-  "desc",
   "eaf.exposure",
   "eaf.outcome",
   "effect",
@@ -51,9 +53,10 @@ utils::globalVariables(c(
   "melt",
   "method",
   "method2",
+  "model",
+  "model_name",
   "mr_keep",
   "mr_keep.outcome",
-  "mr_mode_internal",
   "n",
   "nsnp",
   "nsnp_removed",
@@ -66,7 +69,6 @@ utils::globalVariables(c(
   "outlier_filtered",
   "p",
   "pdf_document",
-  "predict",
   "proxy_snp.outcome",
   "pval",
   "pval.exposure",
diff --git a/R/moe.R b/R/moe.R
index 71d1746e..9b78f248 100644
--- a/R/moe.R
+++ b/R/moe.R
@@ -206,12 +206,12 @@ mr_moe_single <- function(res, rf) {
   pred <- lapply(methodlist, function(m) {
     d <- dplyr::tibble(
       method = m,
-      MOE = predict(rf[[m]], metric, type = "prob")[, 2]
+      MOE = stats::predict(rf[[m]], metric, type = "prob")[, 2]
     )
     return(d)
   }) %>%
     dplyr::bind_rows() %>%
-    dplyr::arrange(desc(MOE))
+    dplyr::arrange(dplyr::desc(MOE))
   if ("MOE" %in% names(res$estimates)) {
     message("Overwriting previous MOE estimate")
     res$estimates <- subset(res$estimates, select = -c(MOE, method2))
diff --git a/R/mr_mode.R b/R/mr_mode.R
index 0247792a..344841d0 100644
--- a/R/mr_mode.R
+++ b/R/mr_mode.R
@@ -373,19 +373,3 @@ mr_simple_mode_nome <- function(b_exp, b_out, se_exp, se_out, parameters = defau
     mode_method = "Simple mode (NOME)"
   ))
 }
-
-
-mr_mode_broken <- function(dat, parameters = default_parameters(), mode_method = "all") {
-  combos <- unique(dat[, c("id.exposure", "exposure", "id.outcome", "outcome")])
-  results <- lapply(seq_len(nrow(combos)), function(i) {
-    x <- dat[dat$id.exposure == combos$id.exposure[i] & dat$id.outcome == combos$id.outcome[i], ]
-    mr_mode_internal(x, parameters, mode_method = "all")
-  })
-  res <- data.table::rbindlist(results, fill = TRUE, use.names = TRUE)
-  data.table::setDF(res)
-  if (mode_method != "all") {
-    return(subset(res, method %in% mode_method))
-  } else {
-    return(res)
-  }
-}
diff --git a/R/multivariable_mr.R b/R/multivariable_mr.R
index cec18d5d..25f13387 100644
--- a/R/multivariable_mr.R
+++ b/R/multivariable_mr.R
@@ -417,8 +417,6 @@ mv_harmonise_data <- function(exposure_dat, outcome_dat, harmonise_strictness =
   keepsnp <- names(tab)[tab == nexp]
   exposure_dat <- subset(exposure_dat, SNP %in% keepsnp)
 
-  exposure_mat <- reshape2::dcast(exposure_dat, SNP ~ id.exposure, value.var = "beta.exposure")
-
   # Get outcome data
   dat <- harmonise_data(
     subset(exposure_dat, id.exposure == exposure_dat$id.exposure[1]),
@@ -428,15 +426,15 @@ mv_harmonise_data <- function(exposure_dat, outcome_dat, harmonise_strictness =
   dat <- subset(dat, mr_keep)
   dat$SNP <- as.character(dat$SNP)
 
-  exposure_beta <- reshape2::dcast(exposure_dat, SNP ~ id.exposure, value.var = "beta.exposure")
+  exposure_beta <- tidyr::pivot_wider(exposure_dat, id_cols = "SNP", names_from = "id.exposure", names_sort = TRUE, values_from = "beta.exposure")
   exposure_beta <- subset(exposure_beta, SNP %in% dat$SNP)
   exposure_beta$SNP <- as.character(exposure_beta$SNP)
 
-  exposure_pval <- reshape2::dcast(exposure_dat, SNP ~ id.exposure, value.var = "pval.exposure")
+  exposure_pval <- tidyr::pivot_wider(exposure_dat, id_cols = "SNP", names_from = "id.exposure", names_sort = TRUE, values_from = "pval.exposure")
   exposure_pval <- subset(exposure_pval, SNP %in% dat$SNP)
   exposure_pval$SNP <- as.character(exposure_pval$SNP)
 
-  exposure_se <- reshape2::dcast(exposure_dat, SNP ~ id.exposure, value.var = "se.exposure")
+  exposure_se <- tidyr::pivot_wider(exposure_dat, id_cols = "SNP", names_from = "id.exposure", names_sort = TRUE, values_from = "se.exposure")
   exposure_se <- subset(exposure_se, SNP %in% dat$SNP)
   exposure_se$SNP <- as.character(exposure_se$SNP)
 
diff --git a/R/query.R b/R/query.R
index 7d7e507e..2f7471a2 100644
--- a/R/query.R
+++ b/R/query.R
@@ -144,7 +144,7 @@ extract_outcome_data_internal <- function(
     n <- length(snps)
     splits <- data.frame(
       snps = snps,
-      chunk_id = rep(1:(ceiling(n / splitsize)), each = splitsize)[1:n]
+      chunk_id = rep(seq_len(ceiling(n / splitsize)), each = splitsize)[seq_len(n)]
     )
     d <- list()
     for (i in seq_along(outcomes)) {
@@ -180,7 +180,7 @@ extract_outcome_data_internal <- function(
     n <- length(outcomes)
     splits <- data.frame(
       outcomes = outcomes,
-      chunk_id = rep(1:(ceiling(n / splitsize)), each = splitsize)[1:n]
+      chunk_id = rep(seq_len(ceiling(n / splitsize)), each = splitsize)[seq_len(n)]
     )
     d <- list()
     for (i in seq_along(snps)) {
diff --git a/R/rucker.R b/R/rucker.R
index f2e510ed..d964b6a6 100644
--- a/R/rucker.R
+++ b/R/rucker.R
@@ -350,8 +350,8 @@ mr_rucker_bootstrap <- function(dat, parameters = default_parameters()) {
   res <- rbind(rucker$rucker, rucker_point, rucker_mean, rucker_median)
   rownames(res) <- NULL
 
-  p1 <- ggplot2::ggplot(bootstrap, ggplot2::aes_string(x = "Q", y = "Qdash")) +
-    ggplot2::geom_point(ggplot2::aes_string(colour = "model")) +
+  p1 <- ggplot2::ggplot(bootstrap, ggplot2::aes(x = Q, y = Qdash)) +
+    ggplot2::geom_point(ggplot2::aes(colour = model)) +
     ggplot2::geom_point(data = subset(bootstrap, i == "Full")) +
     ggplot2::scale_colour_brewer(type = "qual") +
     ggplot2::xlim(0, max(bootstrap$Q, bootstrap$Qdash)) +
@@ -375,11 +375,11 @@ mr_rucker_bootstrap <- function(dat, parameters = default_parameters()) {
   modsel$model_name <- "IVW"
   modsel$model_name[modsel$model %in% c("C", "D")] <- "Egger"
 
-  p2 <- ggplot2::ggplot(modsel, ggplot2::aes_string(x = "Estimate")) +
-    ggplot2::geom_density(ggplot2::aes_string(fill = "model_name"), alpha = 0.4) +
+  p2 <- ggplot2::ggplot(modsel, ggplot2::aes(x = Estimate)) +
+    ggplot2::geom_density(ggplot2::aes(fill = model_name), alpha = 0.4) +
     ggplot2::geom_vline(
       data = res,
-      ggplot2::aes_string(xintercept = "Estimate", colour = "Method")
+      ggplot2::aes(xintercept = Estimate, colour = Method)
     ) +
     ggplot2::scale_colour_brewer(type = "qual") +
     ggplot2::scale_fill_brewer(type = "qual") +
@@ -505,8 +505,8 @@ mr_rucker_jackknife_internal <- function(dat, parameters = default_parameters())
     res <- rbind(rucker$rucker, rucker_point, rucker_mean, rucker_median)
     rownames(res) <- NULL
 
-    p1 <- ggplot2::ggplot(bootstrap, ggplot2::aes_string(x = "Q", y = "Qdash")) +
-      ggplot2::geom_point(ggplot2::aes_string(colour = "model")) +
+    p1 <- ggplot2::ggplot(bootstrap, ggplot2::aes(x = Q, y = Qdash)) +
+      ggplot2::geom_point(ggplot2::aes(colour = model)) +
       ggplot2::geom_point(data = bootstrap[bootstrap$i == "Full", ]) +
       ggplot2::scale_colour_brewer(type = "qual") +
       ggplot2::xlim(0, max(bootstrap$Q, bootstrap$Qdash)) +
@@ -530,11 +530,11 @@ mr_rucker_jackknife_internal <- function(dat, parameters = default_parameters())
     modsel$model_name <- "IVW"
     modsel$model_name[modsel$model %in% c("C", "D")] <- "Egger"
 
-    p2 <- ggplot2::ggplot(modsel, ggplot2::aes_string(x = "Estimate")) +
-      ggplot2::geom_density(ggplot2::aes_string(fill = "model_name"), alpha = 0.4) +
+    p2 <- ggplot2::ggplot(modsel, ggplot2::aes(x = Estimate)) +
+      ggplot2::geom_density(ggplot2::aes(fill = model_name), alpha = 0.4) +
       ggplot2::geom_vline(
         data = res,
-        ggplot2::aes_string(xintercept = "Estimate", colour = "Method")
+        ggplot2::aes(xintercept = Estimate, colour = Method)
       ) +
       ggplot2::scale_colour_brewer(type = "qual") +
       ggplot2::scale_fill_brewer(type = "qual") +
diff --git a/R/zzz.R b/R/zzz.R
index aa19b6ab..11b5477d 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1,3 +1,18 @@
+# Turn any URLs in `text` into clickable terminal hyperlinks (OSC 8). cli only
+# emits the escape sequences where the console supports them and otherwise
+# returns the URL unchanged, so this degrades gracefully.
+linkify <- function(text) {
+  if (!requireNamespace("cli", quietly = TRUE)) {
+    return(text)
+  }
+  # Match URLs but stop before trailing sentence punctuation (e.g. a final ".").
+  m <- gregexpr("https?://[^[:space:]]*[^[:space:][:punct:]]", text, perl = TRUE)
+  regmatches(text, m) <- lapply(regmatches(text, m), function(urls) {
+    vapply(urls, function(u) cli::style_hyperlink(u, u), character(1))
+  })
+  text
+}
+
 .onAttach <- function(libname, pkgname) {
   packageStartupMessage(paste("TwoSampleMR version", utils::packageVersion("TwoSampleMR"), "\n"))
 
@@ -11,7 +26,7 @@
     if (length(b) > 0) {
       o <- lapply(b, function(x) {
         # packageStartupMessage(" Message date: ", x[["date"]])
-        sapply(x[["message"]], function(j) packageStartupMessage(paste(" ", j)))
+        invisible(lapply(x[["message"]], function(j) packageStartupMessage(paste0(" ", linkify(j)))))
       })
     }
   }
diff --git a/vignettes/perform_mr.Rmd b/vignettes/perform_mr.Rmd
index 80fcb658..637b0485 100644
--- a/vignettes/perform_mr.Rmd
+++ b/vignettes/perform_mr.Rmd
@@ -286,7 +286,7 @@ p4[[1]]
 
 A 1-to-many MR analysis interrogates the effect of a single exposure on multiple outcomes or multiple exposures on a single outcome. The results of this analysis can be visualised using the 1-to-many forest plot, with or without stratification on a categorical variable. From a visual point of view, the function works best for 50 or fewer results and is not really designed to handle more than a 100 results. If your number of results is much greater than 50, it may be better to split these across two separate plots. For example, if you have 100 sets of results you could divide these equally across two plots and then combine the two plots together in another program like Powerpoint. The function assumes the results are already in the right order for plotting. As such, users are advised to sort their results according to how they would like them to appear in the plot. Users can use their own code to do this or they can use the `sort_1_to_many()` function.
 
-### Step 1: generate 1-to-many MR results
+### Step 1: Generate 1-to-many MR results
 
 ```{r cache=FALSE, warning=FALSE, eval=FALSE}
 exp_dat <- extract_instruments(outcomes = c("ieu-a-2", "ieu-a-100", "ieu-a-1032", "ieu-a-104", "ieu-a-1", "ieu-a-72", "ieu-a-999"))
@@ -303,7 +303,7 @@ dat2 <- harmonise_data(
 res <- mr(dat2)
 ```
 
-### Step 2. Make the 1-to-many forest plot
+### Step 2: Make the 1-to-many forest plot
 
 #### Example 1. Effect of multiple risk factors on coronary heart disease