# About this analysis

This file runs the usual TF activity analysis from Decoupler analysis.Rmd, then keeps only the transcription factors that regulate MYC, meaning MYC is one of their targets in CollecTRI. It lists them ranked by activity score.

  • Input and methods (unchanged): the DESeq2 Wald stat from DESeq2_mraz.tsv (set with deseq_file in the YAML params), where positive = higher in tFL. Every TF is scored on the full CollecTRI network with decoupleR ULM and MLM (minsize = 5), and on bcellViper with msVIPER with shadow for support.
  • Filter after scoring: only TFs with a CollecTRI edge to MYC are kept, whether they activate or repress it. Scoring the full network first means the MLM and msVIPER scores are identical to the full analysis, and FDR is computed across all TFs.
  • Direction relative to MYC: MYC mRNA is higher in tFL. A TF’s change is consistent with MYC going up if it activates MYC and is more active in tFL, or represses MYC and is less active in tFL.
  • Dependence on MYC itself: each TF is also scored with MYC removed from the signature, to check whether its call rests on its MYC edge.
  • Per-sample link to MYC: each TF’s activity is scored in every sample with the MYC gene left out, then correlated with MYC mRNA across samples (levels) and within patients (changes). Samples with a MYC count of 0 are left out by default (exclude_myc_dropouts), and so are their patients for the within-patient changes, because a zero in a shallow library is more likely dropout than absent MYC.
  • Without MYC target genes: every TF is also scored without two sets of MYC target genes (Hallmark MYC_TARGETS_V1 + V2, and MYC’s CollecTRI targets), to see which TFs are still active beyond MYC’s own program. A summary ranks TFs by how many methods and gene sets agree.
  • Significance: nominal p < 0.05. You can change this in the YAML params.

1 Libraries and paths

suppressPackageStartupMessages({
  library(decoupleR)
  library(viper)
  library(bcellViper)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(patchwork)
  library(openxlsx)
})

proj_dir       <- "C:/Users/User/Desktop/MM Work/Main Work/fl_tFL analysis R"
deseq_file     <- if (grepl("^([A-Za-z]:)?[/\\\\]", params$deseq_file)) params$deseq_file else
                    file.path(proj_dir, params$deseq_file)
collectri_file <- "C:/Users/User/Desktop/MM Work/collectri_human.csv"
out_xlsx       <- if (grepl("^([A-Za-z]:)?[/\\\\]", params$out_xlsx)) params$out_xlsx else
                    file.path(proj_dir, params$out_xlsx)
if (!file.exists(deseq_file))
  stop(deseq_file, " not found. Knit 'Deseq analysis.Rmd' first, or set deseq_file: \"DESeq2_mraz.tsv\".")

# plot tokens
col_act     <- "#2a78d6"   # activates MYC
col_rep     <- "#eb6834"   # represses MYC
col_tfl     <- "#1baf7a"   # tFL samples
col_fl      <- "#898781"   # FL samples
col_context <- "#c3c2b7"
col_grid    <- "#e1e0d9"
col_surface <- "#fcfcfb"
ink_primary <- "#0b0b0b"
ink_second  <- "#52514e"
ink_muted   <- "#898781"
edge_cols   <- c("activates MYC" = col_act, "represses MYC" = col_rep)
# outcome colours for the cross-method tiles (the same slots as in Decoupler proliferation sensitivity.Rmd)
col_survive <- "#2a78d6"
col_atten   <- "#eb6834"
col_lost    <- "#1baf7a"

theme_sens <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(plot.background       = element_rect(fill = col_surface, colour = NA),
          panel.grid.major      = element_line(colour = col_grid, linewidth = 0.3),
          panel.grid.minor      = element_blank(),
          axis.text             = element_text(colour = ink_muted),
          axis.title            = element_text(colour = ink_second),
          plot.title            = element_text(colour = ink_primary, face = "bold"),
          plot.subtitle         = element_text(colour = ink_second),
          plot.caption          = element_text(colour = ink_muted, hjust = 0),
          plot.title.position   = "plot",
          plot.caption.position = "plot",
          strip.text            = element_text(colour = ink_primary, face = "bold", hjust = 0),
          legend.position       = "top",
          legend.justification  = "left",
          legend.title          = element_blank(),
          legend.text           = element_text(colour = ink_second))
}

show_table <- function(df, caption = NULL) {
  num <- which(vapply(df, is.double, logical(1)))
  DT::datatable(df, rownames = FALSE, caption = caption, filter = "top",
                options = list(pageLength = 15, scrollX = TRUE)) |>
    DT::formatSignif(columns = num, digits = 3)
}

is_sig <- function(p) {
  q <- if (params$p_adjust == "BH") stats::p.adjust(p, "BH") else p
  !is.na(q) & q < params$p_cut
}

2 Input signature

Data prep is identical to sections 1 and 3 of Decoupler analysis.Rmd.

deseq_unfiltered <- read.table(deseq_file, sep = "\t", header = TRUE)
deseq <- deseq_unfiltered |> dplyr::filter(baseMean > 4.99)

# decoupleR signature: stat as a genes x 1 matrix (first row per gene_name).
# Labelled "FL_vs_tFL" in Decoupler analysis.Rmd; the label does not affect any result.
deseq_stats <- deseq[!duplicated(deseq$gene_name), c("gene_name", "stat")]
deseq_stats <- deseq_stats[!is.na(deseq_stats$stat), ]
mat_stat <- matrix(deseq_stats$stat, ncol = 1,
                   dimnames = list(deseq_stats$gene_name, "tFL_vs_FL"))

# msVIPER signature: baseMean > 10, one row per gene (largest |stat|)
deg <- deseq_unfiltered |>
  dplyr::filter(!is.na(gene_name) & !is.na(stat) & baseMean > 10) |>
  dplyr::arrange(dplyr::desc(abs(stat))) |>
  dplyr::distinct(gene_name, .keep_all = TRUE)
sign_C <- setNames(deg$stat, deg$gene_name)

# each TF's own DESeq2 results, for context
deseq_lookup <- deseq_unfiltered |>
  dplyr::filter(!duplicated(gene_name)) |>
  dplyr::select(gene_name, deseq_baseMean = baseMean,
                deseq_log2FoldChange = log2FoldChange, deseq_stat = stat,
                deseq_pvalue = pvalue, deseq_padj = padj)

deseq_lookup |> dplyr::filter(gene_name == "MYC") |>
  dplyr::mutate(dplyr::across(where(is.numeric), ~ signif(.x, 3))) |>
  knitr::kable(caption = "MYC itself (positive stat = higher in tFL)")
MYC itself (positive stat = higher in tFL)
gene_name deseq_baseMean deseq_log2FoldChange deseq_stat deseq_pvalue deseq_padj
MYC 54.3 1.28 3.67 0.000243 0.00633

3 TFs that regulate MYC

net_collectri <- read.csv(collectri_file)
net_collectri <- net_collectri |> dplyr::select(source, target, mor = weight)

targets_in_sig <- net_collectri |>
  dplyr::filter(target %in% rownames(mat_stat)) |>
  dplyr::distinct(source, target) |>
  dplyr::count(source, name = "n_targets")

myc_targets <- net_collectri |>
  dplyr::filter(source == "MYC", target %in% rownames(mat_stat)) |>
  dplyr::pull(target)

shared_with_myc <- net_collectri |>
  dplyr::filter(target %in% rownames(mat_stat)) |>
  dplyr::distinct(source, target) |>
  dplyr::group_by(source) |>
  dplyr::summarise(shared_with_MYC_targets = mean(target %in% myc_targets), .groups = "drop")

myc_regulators <- net_collectri |>
  dplyr::filter(target == "MYC") |>
  dplyr::distinct(source, .keep_all = TRUE) |>
  dplyr::transmute(source, mor_to_MYC = mor,
                   edge_to_MYC = ifelse(mor > 0, "activates MYC", "represses MYC")) |>
  dplyr::left_join(targets_in_sig, by = "source") |>
  dplyr::mutate(n_targets = dplyr::coalesce(n_targets, 0L),
                scored    = n_targets >= params$minsize)

myc_regulators |>
  dplyr::count(edge_to_MYC, scored) |>
  tidyr::pivot_wider(names_from = scored, values_from = n, names_prefix = "scored_") |>
  knitr::kable(caption = sprintf("%d CollecTRI TFs have MYC as a target; scored = at least %d targets in the signature",
                                 nrow(myc_regulators), params$minsize))
155 CollecTRI TFs have MYC as a target; scored = at least 5 targets in the signature
edge_to_MYC scored_FALSE scored_TRUE
activates MYC 6 87
represses MYC 2 60 
cat("Not scored (fewer than", params$minsize, "targets in the signature):",
    paste(myc_regulators$source[!myc_regulators$scored], collapse = ", "), "\n")
## Not scored (fewer than 5 targets in the signature): NME2, AHRR, ZBTB32, IKZF3, L3MBTL1, TBL1X, NDN, ZMIZ1
cat("MYC is in the list because CollecTRI has a MYC -> MYC autoregulation edge.\n")
## MYC is in the list because CollecTRI has a MYC -> MYC autoregulation edge.

4 Scoring: the normal analysis, and without MYC

The decouple() call is identical to section 2a of Decoupler analysis.Rmd. It runs on the full signature, and again with the MYC gene removed, which drops the MYC edge from every regulon.

run_collectri <- function(mat) {
  decoupleR::decouple(
    mat = mat, network = net_collectri,
    .source = "source", .target = "target",
    statistics = c("mlm", "ulm"),
    args = list(mlm = list(.mor = "mor"),
                ulm = list(.mor = "mor")),
    consensus_score = FALSE,
    minsize = params$minsize
  )
}

# one row per TF: <method>_score, <method>_p, <method>_fdr (FDR across ALL scored TFs)
to_wide <- function(res, suffix = "") {
  res |>
    dplyr::select(source, statistic, score, p_value) |>
    dplyr::group_by(statistic) |>
    dplyr::mutate(fdr = stats::p.adjust(p_value, "BH")) |>
    dplyr::ungroup() |>
    tidyr::pivot_wider(names_from = statistic, values_from = c(score, p_value, fdr),
                       names_glue = paste0("{statistic}_{.value}", suffix)) |>
    dplyr::rename_with(~ sub("_p_value", "_p", .x))
}

res_all   <- run_collectri(mat_stat)
res_nomyc <- run_collectri(mat_stat[rownames(mat_stat) != "MYC", , drop = FALSE])

scores_all   <- to_wide(res_all)
scores_nomyc <- to_wide(res_nomyc, "_without_MYC") |>
  dplyr::select(source, dplyr::matches("_(score|p)_without_MYC$"))

tibble::tibble(run = c("all genes", "without MYC"),
               TFs_scored = c(dplyr::n_distinct(res_all$source), dplyr::n_distinct(res_nomyc$source))) |>
  knitr::kable()
run TFs_scored
all genes 700
without MYC 698 
data(bcellViper, package = "bcellViper", envir = environment())
regulon_filtered <- regulon[intersect(names(regulon), deg$gene_name)]
mrs <- viper::msviper(sign_C, regulon_filtered, pleiotropy = TRUE, verbose = FALSE)

viper_tab <- tibble::tibble(source      = names(mrs$es$nes),
                            msviper_nes = as.numeric(mrs$es$nes),
                            msviper_p   = as.numeric(mrs$es$p.value)) |>
  dplyr::mutate(msviper_fdr = stats::p.adjust(msviper_p, "BH"),
                MYC_in_bcellviper_regulon = source %in%
                  names(regulon)[vapply(regulon, function(r) "MYC" %in% names(r$tfmode), logical(1))])

6 MYC regulators ranked by activity

myc_table <- myc_regulators |>
  dplyr::filter(scored) |>
  dplyr::left_join(scores_all,      by = "source") |>
  dplyr::left_join(scores_nomyc,    by = "source") |>
  dplyr::left_join(shared_with_myc, by = "source") |>
  dplyr::left_join(viper_tab,       by = "source") |>
  dplyr::left_join(link,            by = "source") |>
  dplyr::left_join(deseq_lookup,    by = c("source" = "gene_name")) |>
  dplyr::mutate(
    ulm_sig = is_sig(ulm_p), mlm_sig = is_sig(mlm_p),
    consistent_with_MYC_up = sign(ulm_score) * mor_to_MYC > 0,
    per_sample_support = dplyr::case_when(
      source == "MYC"                                      ~ "not applicable (autoregulation)",
      is.na(rho_MYC_within_patients)                       ~ NA_character_,
      sign(rho_MYC_within_patients) == mor_to_MYC &
        p_rho_within < 0.05                                ~ "tracks MYC as expected",
      sign(rho_MYC_within_patients) == mor_to_MYC          ~ "expected direction, not significant",
      TRUE                                                 ~ "opposite direction")) |>
  dplyr::arrange(dplyr::desc(ulm_score)) |>
  dplyr::mutate(rank = dplyr::row_number()) |>
  dplyr::select(rank, TF = source, edge_to_MYC, consistent_with_MYC_up, n_targets, shared_with_MYC_targets,
                ulm_score, ulm_p, ulm_fdr, ulm_sig, mlm_score, mlm_p, mlm_fdr, mlm_sig,
                ulm_score_without_MYC, ulm_p_without_MYC, mlm_score_without_MYC, mlm_p_without_MYC,
                msviper_nes, msviper_p, msviper_fdr, MYC_in_bcellviper_regulon,
                rho_MYC_samples, p_rho_samples, rho_MYC_tFL, p_rho_tFL,
                rho_MYC_within_patients, p_rho_within, per_sample_support,
                rho_MYC_within_all_patients,
                change_tFL_vs_FL, p_change, change_adj_for_MYC, p_change_adj, p_MYC_term,
                mor_to_MYC, dplyr::starts_with("deseq_"))

myc_table |>
  dplyr::select(rank, TF, edge_to_MYC, consistent_with_MYC_up, n_targets, ulm_score, ulm_p,
                mlm_score, mlm_p, ulm_score_without_MYC, msviper_nes,
                rho_MYC_within_patients, per_sample_support) |>
  show_table("All scored MYC regulators, sorted by ULM activity score (positive = more active in tFL)")
sig_tab   <- myc_table |> dplyr::filter(ulm_sig)
both_tab  <- myc_table |> dplyr::filter(ulm_sig, mlm_sig)
push_up   <- sig_tab |> dplyr::filter(consistent_with_MYC_up, TF != "MYC")
tracks    <- sig_tab |> dplyr::filter(per_sample_support == "tracks MYC as expected")
max_shift <- max(abs(myc_table$ulm_score - myc_table$ulm_score_without_MYC), na.rm = TRUE)
fmt <- function(d) if (nrow(d)) paste0(d$TF, ifelse(d$ulm_score > 0, " (up)", " (down)"), collapse = ", ") else "none"

At a glance

  • Scored MYC regulators: 147; significant in ULM: 37; in both ULM and MLM: 11.
  • Significant and consistent with MYC going up in tFL (activator more active, or repressor less active): 28: E2F1 (up), TFDP1 (up), FOXM1 (up), MYCN (up), TBP (up), SP1 (up), TFDP2 (up), STAT3 (up), KLF5 (up), JUN (up), EWSR1 (up), ESR1 (up), AR (up), MYBL2 (up), KAT5 (up), SP3 (up), CREBBP (up), AHR (up), FOSB (up), MYB (up), NFKB (up), PATZ1 (up), CTCFL (up), SPI1 (down), RB1 (down), RUNX1 (down), ETV3 (down), IKZF1 (down).
  • Significant and tracking MYC mRNA within patients in the expected direction (9 pairs): 0: none.
  • Removing MYC from the signature changes any regulator’s ULM score by at most 1.05, so no call rests on the MYC edge alone.
plot_df <- myc_table |>
  dplyr::mutate(msviper_sig = is_sig(msviper_p)) |>
  dplyr::filter(ulm_sig | mlm_sig) |>
  dplyr::mutate(TF = factor(TF, levels = TF[order(ulm_score)]))

# one method's score as a lollipop; filled = significant in that method
score_lollipop <- function(score, sig, title, xlab, first = FALSE) {
  d <- plot_df |>
    dplyr::filter(!is.na(.data[[score]])) |>
    dplyr::mutate(x = .data[[score]], s = .data[[sig]])
  p <- ggplot(d, aes(y = TF)) +
    geom_vline(xintercept = 0, colour = col_context, linewidth = 0.3) +
    geom_segment(aes(x = 0, xend = x, yend = TF), colour = col_context, linewidth = 0.4) +
    geom_point(data = ~ dplyr::filter(.x, s), aes(x = x, fill = edge_to_MYC),
               shape = 21, size = 2.6, colour = col_surface, stroke = 0.5) +
    geom_point(data = ~ dplyr::filter(.x, !s), aes(x = x, colour = edge_to_MYC),
               shape = 21, size = 2.4, fill = col_surface, stroke = 1) +
    scale_y_discrete(limits = levels(plot_df$TF)) +   # the same rows in every panel
    scale_fill_manual(values = edge_cols, guide = if (first) "legend" else "none") +
    scale_colour_manual(values = edge_cols, guide = "none") +
    labs(title = title, x = xlab, y = NULL) +
    theme_sens()
  if (first) p else p + theme(axis.text.y = element_blank())
}

p_ulm <- score_lollipop("ulm_score",   "ulm_sig",     "ULM",     "ULM score", first = TRUE)
p_mlm <- score_lollipop("mlm_score",   "mlm_sig",     "MLM",     "MLM score")
p_vip <- score_lollipop("msviper_nes", "msviper_sig", "msVIPER", "msVIPER NES (bcellViper)")

# a Spearman rho as a lollipop; filled = p < 0.05
rho_lollipop <- function(rho, p, title, xlab) {
  d <- plot_df |>
    dplyr::filter(!is.na(.data[[rho]])) |>
    dplyr::mutate(x = .data[[rho]], s = .data[[p]] < 0.05)
  ggplot(d, aes(y = TF)) +
    geom_vline(xintercept = 0, colour = col_context, linewidth = 0.3) +
    geom_segment(aes(x = 0, xend = x, yend = TF), colour = col_context, linewidth = 0.4) +
    geom_point(data = ~ dplyr::filter(.x, s), aes(x = x, fill = edge_to_MYC),
               shape = 21, size = 2.6, colour = col_surface, stroke = 0.5) +
    geom_point(data = ~ dplyr::filter(.x, !s), aes(x = x, colour = edge_to_MYC),
               shape = 21, size = 2.4, fill = col_surface, stroke = 1) +
    scale_y_discrete(limits = levels(plot_df$TF)) +
    scale_fill_manual(values = edge_cols, guide = "none") +
    scale_colour_manual(values = edge_cols, guide = "none") +
    scale_x_continuous(limits = c(-1, 1)) +
    labs(title = title, x = xlab, y = NULL) +
    theme_sens() +
    theme(axis.text.y = element_blank())
}

p_chg <- rho_lollipop("rho_MYC_within_patients", "p_rho_within", "Tracks MYC change",
                      sprintf("Spearman rho, tFL - FL change (%d pairs)", n_pairs_used))
p_lvl <- rho_lollipop(if (lvl_tfl) "rho_MYC_tFL" else "rho_MYC_samples",
                      if (lvl_tfl) "p_rho_tFL" else "p_rho_samples",
                      "Tracks MYC level",
                      sprintf("Spearman rho, %s samples (n = %d)", if (lvl_tfl) "tFL" else "all", n_lvl))

p_ulm + p_mlm + p_vip + p_chg + p_lvl +
  patchwork::plot_layout(widths = c(1.3, 1, 1, 1, 1)) +
  patchwork::plot_annotation(
    title = "MYC regulators significant in ULM or MLM",
    subtitle = paste0("Activity, tFL vs FL (positive = more active in tFL); its within-patient link to the change in MYC mRNA; ",
                      "and its link to the MYC mRNA level across ", if (lvl_tfl) "tFL" else "all", " samples"),
    caption = paste0("Filled points: significant in that method (rho: p < 0.05); hollow: not. ",
                     "No msVIPER point: no bcellViper regulon. Rows sorted by ULM score.\n",
                     "Activators of MYC should sit right of zero and repressors left of zero if their change is consistent with MYC going up in tFL. ",
                     sprintf("|rho| above about %.2f is significant with %d pairs (change) and %.2f with %d samples (level).\n",
                             rho_crit, n_pairs_used, crit_rho(n_lvl), n_lvl),
                     "Samples with 0 MYC reads are left out of the level correlation. MYC's own row reflects autoregulation."),
    theme = theme_sens())

Per-sample activity vs MYC mRNA. One figure per TF. Activity is the ULM t-value in each sample, with the MYC gene left out, and lines join each patient’s FL and tFL samples:

  • Top left: paired activity, FL to tFL, and in how many patients it is higher in tFL.
  • Top right: the same paired activity against MYC mRNA.
  • Bottom left: tFL samples only, against MYC mRNA, with a straight-line trend. The title gives the Spearman correlation of activity with MYC mRNA across tFL samples: whether tFL samples with more MYC have a more active TF, rather than how the two change from FL to tFL. tFL samples with 0 MYC reads are shown hollow and left out of it.
  • Bottom right: the TF’s own mRNA, log2(normalized count + 1), against MYC mRNA. A line’s rise is that patient’s log2 fold change in TF mRNA. The subtitle gives the TF’s DESeq2 log2 fold change from the input file (already shrunken) and its padj.
  • Activity and mRNA are in different units, so they have separate axes. CollecTRI complexes such as NFKB and AP1 have no single gene, so their mRNA panel is empty.

Points at 0 MYC are dropouts (patients 17 and 4), left out of the within-patient correlation. Choose the TFs with per_sample_tfs: in the YAML header: a number such as 15 (the top N regulators above by |ULM score|), "all" with the quotes, or a list of names such as ["GATA3", "RUNX1"].

# which TFs to plot: a number = top N regulators above by |ULM score|, "all" = every one of them,
# or TF names (any TF scored per sample)
sel    <- params$per_sample_tfs
ranked <- plot_df |> dplyr::arrange(dplyr::desc(abs(ulm_score))) |> dplyr::pull(TF) |> as.character()
top_tf <- if (is.numeric(sel)) head(ranked, sel) else if (identical(sel, "all")) ranked else
            intersect(sel, unique(act$source))

ps_all <- act |>
  dplyr::select(source, sample = condition, activity = score) |>
  dplyr::left_join(samples, by = "sample")

# points per sample, lines per patient, on MYC mRNA
patient_panel <- function(df, y, title) {
  ggplot(df, aes(MYC_expr, .data[[y]])) +
    geom_line(aes(group = patient), colour = col_context, linewidth = 0.3) +
    geom_point(aes(fill = group), shape = 21, size = 2.4, colour = col_surface, stroke = 0.4) +
    scale_fill_manual(values = c(FL = col_fl, tFL = col_tfl)) +
    labs(title = title, x = "MYC mRNA, log2(normalized count + 1)", y = NULL) +
    theme_sens()
}

plot_per_sample <- function(tf) {
  d    <- ps_all |> dplyr::filter(source == tf)
  rho  <- pair_rho(d[!d$patient %in% drop_patients, ])
  edge <- myc_regulators$edge_to_MYC[match(tf, myc_regulators$source)]
  row  <- match(tf, deseq_unfiltered$gene_name)   # first row per gene, as in deseq_lookup
  de   <- deseq_lookup[match(tf, deseq_lookup$gene_name), ]

  # paired activity, FL to tFL
  chg    <- tidyr::pivot_wider(d, id_cols = patient, names_from = group, values_from = activity)
  p_pair <- ggplot(d, aes(group, activity)) +
    geom_line(aes(group = patient), colour = col_context, linewidth = 0.3) +
    geom_point(aes(fill = group), shape = 21, size = 2.4, colour = col_surface, stroke = 0.4) +
    scale_fill_manual(values = c(FL = col_fl, tFL = col_tfl)) +
    scale_x_discrete(expand = expansion(add = 0.4)) +
    labs(title = "Activity, FL to tFL (ULM t-value)",
         subtitle = sprintf("Higher in tFL in %d of %d patients", sum(chg$tFL > chg$FL, na.rm = TRUE),
                            sum(!is.na(chg$tFL - chg$FL))),
         x = NULL, y = NULL) +
    theme_sens()

  p_act  <- patient_panel(d, "activity", "Activity vs MYC mRNA") + guides(fill = "none")

  # across tFL samples: activity vs MYC mRNA (a level-to-level link, not a change)
  tfl     <- d |>
    dplyr::filter(group == "tFL") |>
    dplyr::mutate(used = !(isTRUE(params$exclude_myc_dropouts) & MYC_expr == 0))
  rho_tfl <- suppressWarnings(stats::cor.test(tfl$activity[tfl$used], tfl$MYC_expr[tfl$used],
                                              method = "spearman", exact = FALSE))
  p_tfl <- ggplot(tfl, aes(MYC_expr, activity)) +
    geom_smooth(data = ~ dplyr::filter(.x, used), method = "lm", formula = y ~ x, se = FALSE,
                colour = col_context, linewidth = 0.5) +
    geom_point(data = ~ dplyr::filter(.x, used), shape = 21, size = 2.4,
               fill = col_tfl, colour = col_surface, stroke = 0.4) +
    geom_point(data = ~ dplyr::filter(.x, !used), shape = 21, size = 2.2,
               fill = col_surface, colour = col_tfl, stroke = 0.8) +
    labs(title = sprintf("tFL only: rho %+.2f, p = %.2g (n = %d)",
                         unname(rho_tfl$estimate), rho_tfl$p.value, sum(tfl$used)),
         x = "MYC mRNA, log2(normalized count + 1)", y = NULL) +
    theme_sens()

  p_mrna <- if (is.na(row)) {
    ggplot() +
      annotate("text", x = 0, y = 0, label = paste(tf, "is a CollecTRI complex: no single gene"),
               colour = ink_muted, size = 3.2) +
      theme_void() +
      theme(plot.background = element_rect(fill = col_surface, colour = NA))
  } else {
    tibble::tibble(sample = nc, tf_mrna = log2(as.numeric(unlist(deseq_unfiltered[row, nc])) + 1)) |>
      dplyr::left_join(samples, by = "sample") |>
      patient_panel("tf_mrna", paste(tf, "mRNA, log2(normalized count + 1)")) +
      guides(fill = "none")
  }

  mrna_txt <- if (is.na(row) || is.na(de$deseq_log2FoldChange)) "n/a" else
    sprintf("log2FC %+.2f (DESeq2, shrunken), padj %s", de$deseq_log2FoldChange,
            if (is.na(de$deseq_padj)) "n/a" else sprintf("%.2g", de$deseq_padj))

  (p_pair | p_act) / (p_tfl | p_mrna) +
    patchwork::plot_annotation(
      title    = if (is.na(edge)) tf else sprintf("%s (%s)", tf, edge),
      subtitle = sprintf("Activity: ULM %+.1f; within-patient rho with MYC %+.2f, p = %.2g.   mRNA: %s",
                         scores_all$ulm_score[match(tf, scores_all$source)],
                         unname(rho$estimate), rho$p.value, mrna_txt),
      theme    = theme_sens())
}

for (tf in top_tf) print(plot_per_sample(tf))

7 Without MYC target genes

Many TFs share targets with MYC’s own program. When that program goes up in tFL, the shared targets go up with it, and a TF can look active only because of them. This section asks which TFs are still active once MYC’s target genes are taken out. It covers all TFs, with MYC’s regulators highlighted, and uses two definitions of MYC target genes:

  • Hallmark: the MSigDB sets MYC_TARGETS_V1 and MYC_TARGETS_V2, MYC’s core program (ribosome biogenesis, nucleolus, translation, replication).
  • CollecTRI: every target of MYC in CollecTRI, the network the TFs are scored on, whether MYC activates or represses it. This set is much larger.

Each set is tested the same way:

  • Exclusion: the MYC target genes are removed from the signature, so they count as neither TF targets nor background genes. Every TF is rescored with ULM, MLM and msVIPER. This is the design of Decoupler proliferation sensitivity.Rmd, with the MYC target genes in place of the E2F/G2M genes.
  • Random-target null (ULM): a TF’s t-value falls when it loses targets, even ordinary ones. For a TF with m MYC target genes among its k targets, m random targets are removed instead, 1000 times. If the real score falls below this null, the TF leaned on its MYC target genes more than on its average target.
  • Still active means significant (nominal p < 0.05) with the same sign as on all genes. In ULM this includes Attenuated TFs, which stay significant but fall below the null.
  • Summary across methods: each TF gets 6 cells: ULM, MLM and msVIPER, for each gene set. A cell counts if the TF is still active there and points in the TF’s overall direction, taken from ULM if the TF is significant there, otherwise MLM, otherwise msVIPER. TFs are ranked by how many cells count.

7.1 MYC target gene sets

The Hallmark lists below are the MSigDB gene sets MYC_TARGETS_V1 (200 genes) and MYC_TARGETS_V2 (58 genes), human, copied from gsea-msigdb.org on 2026-09-29. Gene names are current HGNC symbols; renamed genes are also matched on the original symbol MSigDB gives for them. The CollecTRI set is every target of MYC in collectri_human.csv. MYC itself is in both sets, so it is removed too. The E2F_TARGETS and G2M_CHECKPOINT lists are the ones in Decoupler proliferation sensitivity.Rmd, used here only to report the overlap.

words <- function(x) strsplit(trimws(x), "\\s+")[[1]]

hallmark_myc_v1 <- words("
  ABCE1 ACP1 AIMP2 AP3S1 APEX1 BUB3 C1QBP CAD CANX CBX3 CCNA2 CCT2 CCT3 CCT4 CCT5
  CCT7 CDC20 CDC45 CDK2 CDK4 CLNS1A CNBP COPS5 COX5A CSTF2 CTPS1 CUL1 CYC1 DDX18 DDX21
  DEK DHX15 DUT EEF1B2 EIF1AX EIF2S1 EIF2S2 EIF3B EIF3D EIF3J EIF4A1 EIF4E EIF4G2 EIF4H EPRS1
  ERH ETF1 EXOSC7 FAM120A FBL G3BP1 GLO1 GNL3 GOT2 GSPT1 H2AZ1 HDAC2 HDDC2 HDGF HNRNPA1
  HNRNPA2B1 HNRNPA3 HNRNPC HNRNPD HNRNPR HNRNPU HPRT1 HSP90AB1 HSPD1 HSPE1 IARS1 IFRD1 ILF2 IMPDH2 KARS1
  KPNA2 KPNB1 LDHA LSM2 LSM7 MAD2L1 MCM2 MCM4 MCM5 MCM6 MCM7 MRPL23 MRPL9 MRPS18B MYC
  NAP1L1 NCBP1 NCBP2 NDUFAB1 NHP2 NME1 NOLC1 NOP16 NOP56 NPM1 ODC1 ORC2 PA2G4 PABPC1 PABPC4
  PCBP1 PCNA PGK1 PHB1 PHB2 POLD2 POLE3 PPIA PPM1G PRDX3 PRDX4 PRPF31 PRPS2 PSMA1 PSMA2
  PSMA4 PSMA6 PSMA7 PSMB2 PSMB3 PSMC4 PSMC6 PSMD1 PSMD14 PSMD3 PSMD7 PSMD8 PTGES3 PWP1 RACK1
  RAD23B RAN RANBP1 RFC4 RNPS1 RPL14 RPL18 RPL22 RPL34 RPL6 RPLP0 RPS10 RPS2 RPS3 RPS5
  RPS6 RRM1 RRP9 RSL1D1 RUVBL2 SERBP1 SET SF3A1 SF3B3 SLC25A3 SMARCC1 SNRPA SNRPA1 SNRPB2 SNRPD1
  SNRPD2 SNRPD3 SNRPG SRM SRPK1 SRSF1 SRSF2 SRSF3 SRSF7 SSB SSBP1 STARD7 SYNCRIP TARDBP TCP1
  TFDP1 TOMM70 TRA2B TRIM28 TUFM TXNL4A TYMS U2AF1 UBA2 UBE2E1 UBE2L3 USP1 VBP1 VDAC1 VDAC3
  XPO1 XPOT XRCC6 YWHAE YWHAQ")

hallmark_myc_v2 <- words("
  AIMP2 BYSL CBX3 CDK4 DCTPP1 DDX18 DUSP2 EXOSC5 FARSA GNL3 GRWD1 HK2 HSPD1 HSPE1 IMP4
  IPO4 LAS1L MAP3K6 MCM4 MCM5 MPHOSPH10 MRTO4 MYBBP1A MYC NDUFAF4 NIP7 NOC4L NOLC1 NOP16 NOP2
  NOP56 NPM1 PA2G4 PES1 PHB1 PLK1 PLK4 PPAN PPRC1 PRMT3 PUS1 RABEPK RCL1 RRP12 RRP9
  SLC19A1 SLC29A2 SORD SRM SUPV3L1 TBRG4 TCOF1 TFB2M TMEM97 UNG UTP20 WDR43 WDR74")

hallmark_e2f <- words("
  AK2 ANP32E ASF1A ASF1B ATAD2 AURKA AURKB BARD1 BIRC5 BRCA1 BRCA2 BRMS1L BUB1B CBX5 CCNB2
  CCNE1 CCP110 CDC20 CDC25A CDC25B CDCA3 CDCA8 CDK1 CDK4 CDKN1A CDKN1B CDKN2A CDKN2C CDKN3
  CENPE CENPM CHEK1 CHEK2 CIT CKS1B CKS2 CNOT9 CSE1L CTCF CTPS1 DCK DCLRE1B DCTPP1 DDX39A
  DEK DEPDC1 DIAPH3 DLGAP5 DNMT1 DONSON DSCC1 DUT E2F8 EED EIF2S1 ESPL1 EXOSC8 EZH2 GINS1
  GINS3 GINS4 GSPT1 H2AX H2AZ1 HELLS HMGA1 HMGB2 HMGB3 HMMR HNRNPD HUS1 ILF3 ING3 IPO7
  JPT1 KIF18B KIF22 KIF2C KIF4A KPNA2 LBR LIG1 LMNB1 LUC7L3 LYAR MAD2L1 MCM2 MCM3 MCM4
  MCM5 MCM6 MCM7 MELK MKI67 MLH1 MMS22L MRE11 MSH2 MTHFD2 MXD3 MYBL2 MYC NAA38 NAP1L1 NASP
  NBN NCAPD2 NME1 NOLC1 NOP56 NUDT21 NUP107 NUP153 NUP205 ORC2 ORC6 PA2G4 PAICS PAN2 PCNA
  PDS5B PHF5A PLK1 PLK4 PMS2 PNN POLA2 POLD1 POLD2 POLD3 POLE POLE4 POP7 PPM1D PPP1R8
  PRDX4 PRIM2 PRKDC PRPS1 PSIP1 PSMC3IP PTTG1 RACGAP1 RAD1 RAD21 RAD50 RAD51AP1 RAD51C RAN
  RANBP1 RBBP7 RFC1 RFC2 RFC3 RNASEH2A RPA1 RPA2 RPA3 RRM2 SHMT1 SLBP SMC1A SMC3 SMC4 SMC6
  SNRPB SPAG5 SPC24 SPC25 SRSF1 SRSF2 SSRP1 STAG1 STMN1 SUV39H1 SYNCRIP TACC3 TBRG4 TCF19
  TFRC TIMELESS TIPIN TK1 TMPO TOP2A TP53 TRA2B TRIP13 TUBB TUBG1 UBE2S UBE2T UBR7 UNG
  USP1 WDR90 WEE1 XPO1 XRCC6 ZW10")

hallmark_g2m <- words("
  ABL1 AMD1 ARID4A ATF5 ATRX AURKA AURKB BARD1 BCL3 BIRC5 BRCA2 BUB1 BUB3 CASP8AP2 CBX1
  CCNA2 CCNB2 CCND1 CCNF CCNT1 CDC20 CDC25A CDC25B CDC27 CDC45 CDC6 CDC7 CDK1 CDK4 CDKN1B
  CDKN2C CDKN3 CENPA CENPE CENPF CHAF1A CHEK1 CHMP1A CKS1B CKS2 CTCF CUL1 CUL3 CUL4A CUL5
  DBF4 DDX39A DKC1 DMD DR1 DTYMK E2F1 E2F2 E2F3 E2F4 EFNA5 EGF ESPL1 EWSR1 EXO1 EZH2 FANCC
  FBXO5 FOXN3 G3BP1 GINS2 GSPT1 H2AX H2AZ1 H2AZ2 H2BC12 HIF1A HIRA HMGA1 HMGB3 HMGN2 HMMR
  HNRNPD HNRNPU HOXC10 HSPA8 HUS1 ILF3 INCENP JPT1 KATNA1 KIF11 KIF15 KIF20B KIF22 KIF23
  KIF2C KIF4A KIF5B KMT5A KNL1 KPNA2 KPNB1 LBR LIG3 LMNB1 MAD2L1 MAP3K20 MAPK14 MARCKS
  MCM2 MCM3 MCM5 MCM6 MEIS1 MEIS2 MKI67 MNAT1 MT2A MTF2 MYBL2 MYC NASP NCL NDC80 NEK2
  NOLC1 NOTCH2 NSD2 NUMA1 NUP50 NUP98 NUSAP1 ODC1 ODF2 ORC5 ORC6 PAFAH1B1 PBK PDS5B PLK1
  PLK4 PML POLA2 POLE POLQ PRC1 PRIM2 PRMT5 PRP4K PTTG1 PTTG3P PURA RACGAP1 RAD21 RAD23B
  RAD54L RASAL2 RBL1 RBM14 RPA2 RPS6KA5 SAP30 SFPQ SLC12A2 SLC38A1 SLC7A1 SLC7A5 SMAD3
  SMARCC1 SMC1A SMC2 SMC4 SNRPD1 SQLE SRSF1 SRSF10 SRSF2 SS18 STAG1 STIL STMN1 SUV39H1
  SYNCRIP TACC3 TENT4A TFDP1 TGFB1 TLE3 TMPO TNPO2 TOP1 TOP2A TPX2 TRA2B TRAIP TROAP TTK
  UBE2C UBE2S UCK2 UPF1 WRN XPO1 YTHDC1")

stopifnot(length(unique(hallmark_myc_v1)) == 200, length(unique(hallmark_myc_v2)) == 58,
          length(unique(hallmark_e2f)) == 200, length(unique(hallmark_g2m)) == 200)

# current HGNC symbol -> original MSigDB source symbol, for renamed genes
msigdb_alias <- c(EPRS1 = "EPRS",   IARS1 = "IARS",  KARS1 = "KARS",       PHB1 = "PHB",
                  CTPS1 = "CTPS",   H2AX = "H2AFX",  H2AZ1 = "H2AFZ",      H2AZ2 = "H2AFV",
                  H2BC12 = "HIST1H2BK", JPT1 = "HN1", MRE11 = "MRE11A",    CNOT9 = "RQCD1",
                  KNL1 = "CASC5",   TENT4A = "PAPD7", PRP4K = "PRPF4B",    KMT5A = "SETD8",
                  NSD2 = "WHSC1",   MAP3K20 = "ZAK")

to_data_symbol <- function(g, pool) {
  hit <- ifelse(g %in% pool, g, NA_character_)
  alt <- unname(msigdb_alias[g])
  ifelse(is.na(hit) & !is.na(alt) & alt %in% pool, alt, hit)
}

e2f_g2m      <- union(hallmark_e2f, hallmark_g2m)
prolif_genes <- sort(unique(stats::na.omit(to_data_symbol(e2f_g2m, deseq_unfiltered$gene_name))))
n_prolif_sig <- sum(prolif_genes %in% rownames(mat_stat))

myc_set_tbl <- tibble::tibble(msigdb_symbol = sort(union(hallmark_myc_v1, hallmark_myc_v2))) |>
  dplyr::mutate(in_MYC_TARGETS_V1 = msigdb_symbol %in% hallmark_myc_v1,
                in_MYC_TARGETS_V2 = msigdb_symbol %in% hallmark_myc_v2,
                in_E2F_TARGETS    = msigdb_symbol %in% hallmark_e2f,
                in_G2M_CHECKPOINT = msigdb_symbol %in% hallmark_g2m,
                data_symbol       = to_data_symbol(msigdb_symbol, deseq_unfiltered$gene_name))
myc_set_genes <- sort(unique(stats::na.omit(myc_set_tbl$data_symbol)))

# CollecTRI: every target of MYC, activated or repressed (the MYC -> MYC edge included)
collectri_myc <- net_collectri |>
  dplyr::filter(source == "MYC") |>
  dplyr::distinct(target, .keep_all = TRUE) |>
  dplyr::select(gene = target, collectri_MYC_mor = mor)

myc_sets   <- list(hallmark = myc_set_genes, collectri = collectri_myc$gene)
set_labels <- c(hallmark = "Hallmark MYC_TARGETS_V1 + V2", collectri = "CollecTRI MYC targets")
set_keys   <- setNames(names(set_labels), set_labels)

# one row per MYC target gene of either set, for the workbook
e2f_data <- to_data_symbol(hallmark_e2f, deseq_unfiltered$gene_name)
g2m_data <- to_data_symbol(hallmark_g2m, deseq_unfiltered$gene_name)
myc_gene_tbl <- myc_set_tbl |>
  dplyr::transmute(gene = dplyr::coalesce(data_symbol, msigdb_symbol), msigdb_symbol,
                   in_MYC_TARGETS_V1, in_MYC_TARGETS_V2) |>
  dplyr::full_join(collectri_myc, by = "gene") |>
  dplyr::mutate(in_MYC_TARGETS_V1      = dplyr::coalesce(in_MYC_TARGETS_V1, FALSE),
                in_MYC_TARGETS_V2      = dplyr::coalesce(in_MYC_TARGETS_V2, FALSE),
                collectri_MYC_target   = !is.na(collectri_MYC_mor),
                in_E2F_TARGETS         = gene %in% e2f_data,
                in_G2M_CHECKPOINT      = gene %in% g2m_data,
                in_decoupler_signature = gene %in% rownames(mat_stat),
                in_msviper_signature   = gene %in% names(sign_C)) |>
  dplyr::left_join(deseq_lookup |> dplyr::select(gene_name, deseq_baseMean, deseq_stat),
                   by = c("gene" = "gene_name")) |>
  dplyr::arrange(gene)

overlap_row <- function(set_name, s) {
  tibble::tibble(MYC_gene_set        = set_name,
                 genes               = length(s),
                 also_E2F_TARGETS    = sum(s %in% hallmark_e2f),
                 also_G2M_CHECKPOINT = sum(s %in% hallmark_g2m),
                 also_E2F_or_G2M     = sum(s %in% e2f_g2m),
                 pct_E2F_or_G2M      = 100 * mean(s %in% e2f_g2m))
}
dplyr::bind_rows(overlap_row("MYC_TARGETS_V1", hallmark_myc_v1),
                 overlap_row("MYC_TARGETS_V2", hallmark_myc_v2),
                 overlap_row("V1 + V2 (used here)", union(hallmark_myc_v1, hallmark_myc_v2))) |>
  knitr::kable(digits = 1, caption = "Overlap of the Hallmark MYC sets with E2F_TARGETS and G2M_CHECKPOINT (all MSigDB genes)")
Overlap of the Hallmark MYC sets with E2F_TARGETS and G2M_CHECKPOINT (all MSigDB genes)
MYC_gene_set genes also_E2F_TARGETS also_G2M_CHECKPOINT also_E2F_or_G2M pct_E2F_or_G2M
MYC_TARGETS_V1 200 35 29 47 23.5
MYC_TARGETS_V2 58 12 6 12 20.7
V1 + V2 (used here) 240 40 31 52 21.7 
shared_genes <- myc_set_tbl |> dplyr::filter(in_E2F_TARGETS | in_G2M_CHECKPOINT)
cat("Hallmark MYC target genes that are also E2F/G2M genes (", nrow(shared_genes), "): ",
    paste(shared_genes$msigdb_symbol, collapse = ", "), "\n", sep = "")
## Hallmark MYC target genes that are also E2F/G2M genes (52): BUB3, CCNA2, CDC20, CDC45, CDK4, CTPS1, CUL1, DCTPP1, DEK, DUT, EIF2S1, G3BP1, GSPT1, H2AZ1, HNRNPD, HNRNPU, KPNA2, KPNB1, MAD2L1, MCM2, MCM4, MCM5, MCM6, MCM7, MYC, NAP1L1, NME1, NOLC1, NOP56, ODC1, ORC2, PA2G4, PCNA, PLK1, PLK4, POLD2, PRDX4, RAD23B, RAN, RANBP1, SMARCC1, SNRPD1, SRSF1, SRSF2, SYNCRIP, TBRG4, TFDP1, TRA2B, UNG, USP1, XPO1, XRCC6
renamed <- myc_set_tbl |> dplyr::filter(!is.na(data_symbol), data_symbol != msigdb_symbol)
cat("Hallmark genes matched on the older symbol:",
    paste(renamed$msigdb_symbol, renamed$data_symbol, sep = " -> ", collapse = ", "), "\n")
## Hallmark genes matched on the older symbol: EPRS1 -> EPRS, H2AZ1 -> H2AFZ, IARS1 -> IARS, KARS1 -> KARS, PHB1 -> PHB
cat("Hallmark genes not in the DESeq2 file:",
    paste(myc_set_tbl$msigdb_symbol[is.na(myc_set_tbl$data_symbol)], collapse = ", "), "\n")
## Hallmark genes not in the DESeq2 file: EIF4A1, HNRNPA1, HSPE1, NME1, PPAN, PSMA2, PSMA6, RPS10, TARDBP
sets_tbl <- dplyr::bind_rows(lapply(names(myc_sets), function(key) {
  in_sig <- intersect(myc_sets[[key]], rownames(mat_stat))
  other  <- myc_sets[[setdiff(names(myc_sets), key)]]
  tibble::tibble(MYC_gene_set           = set_labels[[key]],
                 in_decoupleR_signature = length(in_sig),
                 in_msVIPER_signature   = sum(names(sign_C) %in% myc_sets[[key]]),
                 also_E2F_G2M           = sum(in_sig %in% prolif_genes),
                 pct_of_E2F_G2M_removed = 100 * sum(in_sig %in% prolif_genes) / n_prolif_sig,
                 also_in_the_other_set  = sum(in_sig %in% other))
}))
knitr::kable(sets_tbl, digits = 1,
             caption = sprintf("The two MYC target gene sets, genes in the signatures (the decoupleR signature has %d E2F/G2M genes)",
                               n_prolif_sig))
The two MYC target gene sets, genes in the signatures (the decoupleR signature has 319 E2F/G2M genes)
MYC_gene_set in_decoupleR_signature in_msVIPER_signature also_E2F_G2M pct_of_E2F_G2M_removed also_in_the_other_set
Hallmark MYC_TARGETS_V1 + V2 229 228 50 15.7 56
CollecTRI MYC targets 755 720 80 25.1 56

How the gene groups move between FL and tFL:

group_levels <- c("MYC target gene, also E2F/G2M", "MYC target gene, not E2F/G2M",
                  "E2F/G2M gene, not MYC target", "All other genes")

shift_tbl <- function(key) {
  s <- myc_sets[[key]]
  tibble::tibble(gene = rownames(mat_stat), stat = mat_stat[, 1]) |>
    dplyr::mutate(group = dplyr::case_when(
                    gene %in% s & gene %in% prolif_genes ~ group_levels[1],
                    gene %in% s                          ~ group_levels[2],
                    gene %in% prolif_genes               ~ group_levels[3],
                    TRUE                                 ~ group_levels[4]),
                  group = factor(group, levels = group_levels)) |>
    dplyr::group_by(group) |>
    dplyr::summarise(genes = dplyr::n(), mean_stat = mean(stat), median_stat = stats::median(stat),
                     pct_higher_in_tFL = 100 * mean(stat > 0), .groups = "drop") |>
    dplyr::mutate(MYC_gene_set = set_labels[[key]], .before = 1)
}

dplyr::bind_rows(lapply(names(myc_sets), shift_tbl)) |>
  knitr::kable(digits = 2, caption = "DESeq2 stat by gene group, decoupleR signature (positive = higher in tFL)")
DESeq2 stat by gene group, decoupleR signature (positive = higher in tFL)
MYC_gene_set group genes mean_stat median_stat pct_higher_in_tFL
Hallmark MYC_TARGETS_V1 + V2 MYC target gene, also E2F/G2M 50 3.43 3.44 98.00
Hallmark MYC_TARGETS_V1 + V2 MYC target gene, not E2F/G2M 179 2.60 2.63 96.65
Hallmark MYC_TARGETS_V1 + V2 E2F/G2M gene, not MYC target 269 2.73 2.83 88.48
Hallmark MYC_TARGETS_V1 + V2 All other genes 19129 0.00 -0.07 48.09
CollecTRI MYC targets MYC target gene, also E2F/G2M 80 2.96 3.08 86.25
CollecTRI MYC targets MYC target gene, not E2F/G2M 675 0.71 0.67 63.56
CollecTRI MYC targets E2F/G2M gene, not MYC target 239 2.80 2.86 91.21
CollecTRI MYC targets All other genes 18633 0.00 -0.07 48.00

7.2 Scoring without MYC target genes

The same decouple() and msviper() calls as above, run once for each gene set on the signatures without it.

The random-target null uses a closed-form ULM, copied from Decoupler proliferation sensitivity.Rmd. It computes the same t-value as decoupleR::run_ulm, which is checked on all genes and after each exclusion. MYC target genes that are not targets of the TF are removed in every permutation, as in the real exclusion, so the only difference is which targets are removed.

y     <- mat_stat[, 1]
genes <- rownames(mat_stat)

net_u <- net_collectri |>
  dplyr::filter(target %in% genes) |>
  dplyr::distinct(source, target, .keep_all = TRUE)
tf_keep <- net_u |> dplyr::count(source) |> dplyr::filter(n >= params$minsize) |> dplyr::pull(source)
net_u   <- net_u |> dplyr::filter(source %in% tf_keep)
tfs     <- sort(unique(net_u$source))
W <- Matrix::sparseMatrix(i = match(net_u$target, genes), j = match(net_u$source, tfs),
                          x = net_u$mor, dims = c(length(genes), length(tfs)),
                          dimnames = list(genes, tfs))

# ULM t-value of the slope in y ~ w (with intercept), from sufficient statistics
ulm_t <- function(n, sy, syy, sw, sww, swy) {
  sxx <- sww - sw^2 / n
  sxy <- swy - sw * sy / n
  syc <- syy - sy^2 / n
  r   <- sxy / sqrt(sxx * syc)
  r * sqrt((n - 2) / (1 - r^2))
}
ulm_closed <- function(y, W) {
  ulm_t(length(y), sum(y), sum(y^2), Matrix::colSums(W), Matrix::colSums(W^2),
        as.numeric(Matrix::crossprod(W, y)))
}

ulm_before  <- res_all |> dplyr::filter(statistic == "ulm")
t_cf_before <- setNames(ulm_closed(y, W), tfs)
if (!setequal(ulm_before$source, tfs) ||
    max(abs(ulm_before$score - t_cf_before[ulm_before$source])) > 1e-6)
  stop("Closed-form ULM does not reproduce decoupleR; the random-target null would not be valid.")

# size-matched random-target null for every TF; is_x marks the excluded genes
ulm_null <- function(is_x) {
  set.seed(params$seed)
  n_all <- length(y);  sy_all <- sum(y);        syy_all <- sum(y^2)
  n_x   <- sum(is_x);  sy_x   <- sum(y[is_x]);  syy_x   <- sum(y[is_x]^2)
  null_cols <- c("k", "m", "t_before", "t_after", "retention",
                 "null_lo", "null_med", "null_hi", "p_null")

  null_one <- function(j) {
    out <- setNames(rep(NA_real_, length(null_cols)), null_cols)
    rng <- if (W@p[j + 1] > W@p[j]) (W@p[j] + 1):W@p[j + 1] else integer(0)
    idx <- W@i[rng] + 1L
    w   <- W@x[rng]
    yj  <- y[idx]
    xj  <- is_x[idx]
    k   <- length(idx); m <- sum(xj)
    sw  <- sum(w); sww <- sum(w^2); swy <- sum(w * yj)
    out[c("k", "m")] <- c(k, m)
    out["t_before"]  <- ulm_t(n_all, sy_all, syy_all, sw, sww, swy)
    if (k - m < params$minsize) return(out)

    # universe with every excluded gene that is NOT a target of this TF removed
    nA   <- n_all   - (n_x - m)
    syA  <- sy_all  - (sy_x  - sum(yj[xj]))
    syyA <- syy_all - (syy_x - sum(yj[xj]^2))
    t_drop <- function(s_y, s_yy, s_w, s_ww, s_wy)
      ulm_t(nA - m, syA - s_y, syyA - s_yy, sw - s_w, sww - s_ww, swy - s_wy)

    out["t_after"]   <- t_drop(sum(yj[xj]), sum(yj[xj]^2), sum(w[xj]), sum(w[xj]^2), sum(w[xj] * yj[xj]))
    out["retention"] <- out[["t_after"]] / out[["t_before"]]
    if (m == 0) return(out)

    R  <- matrix(replicate(params$n_perm, sample.int(k, m)), nrow = m)   # m x n_perm
    cs <- function(v) colSums(matrix(v[R], nrow = m))
    ret_null <- t_drop(cs(yj), cs(yj^2), cs(w), cs(w^2), cs(w * yj)) / out[["t_before"]]

    out[c("null_lo", "null_med", "null_hi")] <- stats::quantile(ret_null, c(0.025, 0.5, 0.975),
                                                                names = FALSE)
    out["p_null"] <- (1 + sum(ret_null <= out[["retention"]])) / (params$n_perm + 1)
    out
  }

  do.call(rbind, lapply(seq_along(tfs), null_one)) |>
    tibble::as_tibble() |>
    dplyr::mutate(source = tfs, .before = 1)
}

classify <- function(score_b, sig_b, score_a, sig_a) {
  testable_a <- !is.na(score_a)
  dplyr::case_when(
    is.na(score_b)                                 ~ NA_character_,   # TF not in this network
    sig_b & !testable_a                            ~ "Untestable after exclusion",
    sig_b & sig_a & sign(score_a) == sign(score_b) ~ "Survives",
    sig_b & sig_a                                  ~ "Lost (direction reversed)",
    sig_b                                          ~ "Lost",
    testable_a & sig_a                             ~ "Emerged after exclusion",
    TRUE                                           ~ "Not significant")
}
run_without <- function(key) {
  excl  <- myc_sets[[key]]
  is_x  <- genes %in% excl
  res   <- run_collectri(mat_stat[!is_x, , drop = FALSE])
  mrs_x <- viper::msviper(sign_C[!names(sign_C) %in% excl], regulon_filtered, pleiotropy = TRUE, verbose = FALSE)
  nul   <- ulm_null(is_x)

  # the closed-form exclusion inside the null must equal decoupleR's run without these genes
  chk <- nul |>
    dplyr::inner_join(res |> dplyr::filter(statistic == "ulm") |> dplyr::select(source, score), by = "source")
  if (nrow(chk) != sum(!is.na(nul$t_after)) || max(abs(chk$t_after - chk$score)) > 1e-6)
    stop("Closed-form ULM does not reproduce decoupleR without the ", set_labels[[key]], ".")

  list(scores = to_wide(res),
       viper  = tibble::tibble(source      = names(mrs_x$es$nes),
                               msviper_nes = as.numeric(mrs_x$es$nes),
                               msviper_p   = as.numeric(mrs_x$es$p.value),
                               msviper_fdr = stats::p.adjust(msviper_p, "BH")),
       null   = nul,
       size   =tibble::tibble(MYC_gene_set            = set_labels[[key]],
                               decoupleR_signature     = sum(!is_x),
                               msVIPER_signature       = sum(!names(sign_C) %in% excl),
                               TFs_scored_CollecTRI    = dplyr::n_distinct(res$source),
                               regulons_scored_msVIPER = length(mrs_x$es$nes)))
}

runs     <- lapply(setNames(names(myc_sets), names(myc_sets)), run_without)
n_untest <- vapply(runs, function(r) sum(is.na(r$null$t_after)), integer(1))

dplyr::bind_rows(
  tibble::tibble(MYC_gene_set = "none (all genes)", decoupleR_signature = nrow(mat_stat),
                 msVIPER_signature = length(sign_C), TFs_scored_CollecTRI = dplyr::n_distinct(res_all$source),
                 regulons_scored_msVIPER = length(mrs$es$nes)),
  dplyr::bind_rows(lapply(runs, `[[`, "size"))) |>
  knitr::kable(caption = "Signatures and TFs scored with and without each MYC target gene set")
Signatures and TFs scored with and without each MYC target gene set
MYC_gene_set decoupleR_signature msVIPER_signature TFs_scored_CollecTRI regulons_scored_msVIPER
none (all genes) 19627 17093 700 529
Hallmark MYC_TARGETS_V1 + V2 19398 16865 694 529
CollecTRI MYC targets 18872 16373 574 529

7.3 Which TFs are still active: summary across methods and gene sets

# one row per TF on all genes, CollecTRI and bcellViper, with its overall direction
tf_base <- scores_all |>
  dplyr::full_join(viper_tab |> dplyr::select(source, msviper_nes, msviper_p, msviper_fdr), by = "source") |>
  dplyr::left_join(targets_in_sig, by = "source") |>
  dplyr::left_join(myc_regulators |> dplyr::select(source, edge_to_MYC), by = "source") |>
  dplyr::left_join(myc_table |> dplyr::select(source = TF, consistent_with_MYC_up), by = "source") |>
  dplyr::mutate(MYC_regulator = source %in% myc_regulators$source,
                ulm_sig       = is_sig(ulm_p),
                mlm_sig       = is_sig(mlm_p),
                msviper_sig   = is_sig(msviper_p),
                direction     = dplyr::case_when(ulm_sig     ~ sign(ulm_score),
                                                 mlm_sig     ~ sign(mlm_score),
                                                 msviper_sig ~ sign(msviper_nes)),
                direction     = dplyr::case_when(direction > 0 ~ "up in tFL", direction < 0 ~ "down in tFL")) |>
  dplyr::rename(TF = source)

verdicts <- function(key) {
  r <- runs[[key]]
  tf_base |>
    dplyr::select(TF, ulm_score, ulm_sig, mlm_score, mlm_sig, msviper_nes, msviper_sig) |>
    dplyr::left_join(r$scores |> dplyr::select(TF = source, ulm_score_after = ulm_score, ulm_p_after = ulm_p,
                                               ulm_fdr_after = ulm_fdr, mlm_score_after = mlm_score,
                                               mlm_p_after = mlm_p, mlm_fdr_after = mlm_fdr), by = "TF") |>
    dplyr::left_join(r$null |> dplyr::select(TF = source, n_targets_removed = m, retention,
                                             null_lo, null_med, null_hi, p_null), by = "TF") |>
    dplyr::left_join(r$viper |> dplyr::select(TF = source, msviper_nes_after = msviper_nes,
                                              msviper_p_after = msviper_p, msviper_fdr_after = msviper_fdr),
                     by = "TF") |>
    dplyr::mutate(
      MYC_gene_set  = set_labels[[key]],
      ulm_call      = classify(ulm_score, ulm_sig, ulm_score_after, is_sig(ulm_p_after)),
      ulm_verdict   = dplyr::case_when(
        ulm_call == "Survives" & n_targets_removed == 0 ~ "Survives (no MYC target genes among its targets)",
        ulm_call == "Survives" & retention >= null_lo   ~ "Survives: robust",
        ulm_call == "Survives"                          ~ "Attenuated: still significant, MYC-target-weighted",
        ulm_call == "Lost (direction reversed)"         ~ "Reversed: significant in the opposite direction",
        ulm_call == "Lost" & n_targets_removed == 0     ~ "Lost: background shift only (no MYC target genes among its targets)",
        ulm_call == "Lost" & retention < null_lo        ~ "Lost: MYC-target-driven",
        ulm_call == "Lost"                              ~ "Lost: no worse than random target loss",
        TRUE                                            ~ ulm_call),
      still_active  = dplyr::if_else(ulm_sig, ulm_call == "Survives", NA),   # ULM, Attenuated included
      mlm_call      = classify(mlm_score, mlm_sig, mlm_score_after, is_sig(mlm_p_after)),
      msviper_call  = classify(msviper_nes, msviper_sig, msviper_nes_after, is_sig(msviper_p_after))) |>
    dplyr::select(TF, MYC_gene_set, n_targets_removed, ulm_score_after, ulm_p_after, ulm_fdr_after,
                  retention, null_lo, null_med, null_hi, p_null, ulm_verdict, still_active,
                  mlm_score_after, mlm_p_after, mlm_fdr_after, mlm_call,
                  msviper_nes_after, msviper_p_after, msviper_fdr_after, msviper_call)
}

# one row per TF and gene set
ex_long <- dplyr::bind_rows(lapply(names(myc_sets), verdicts)) |>
  dplyr::left_join(tf_base, by = "TF") |>
  dplyr::relocate(TF, MYC_gene_set, MYC_regulator, edge_to_MYC, consistent_with_MYC_up, direction, n_targets,
                  ulm_score, ulm_p, ulm_fdr, ulm_sig, mlm_score, mlm_p, mlm_fdr, mlm_sig,
                  msviper_nes, msviper_p, msviper_fdr, msviper_sig)

# one row per TF, gene set and method; a cell counts if still active in the TF's overall direction
method_levels <- c("ULM", "MLM", "msVIPER")
cells <- dplyr::bind_rows(
  ex_long |> dplyr::transmute(TF, MYC_gene_set, method = "ULM",      outcome = ulm_verdict,
                              sig_before = ulm_sig,     sign_before = sign(ulm_score)),
  ex_long |> dplyr::transmute(TF, MYC_gene_set, method = "MLM",      outcome = mlm_call,
                              sig_before = mlm_sig,     sign_before = sign(mlm_score)),
  ex_long |> dplyr::transmute(TF, MYC_gene_set, method = "msVIPER",  outcome = msviper_call,
                              sig_before = msviper_sig, sign_before = sign(msviper_nes))) |>
  dplyr::left_join(tf_base |> dplyr::select(TF, direction), by = "TF") |>
  dplyr::mutate(
    dir_sign = dplyr::case_when(direction == "up in tFL" ~ 1, direction == "down in tFL" ~ -1),
    agrees   = dplyr::coalesce(sign_before == dir_sign, FALSE),
    code     = dplyr::case_when(
      is.na(outcome)                          ~ "Not in network",
      sig_before & !agrees                    ~ "Opposite direction",
      grepl("^Survives", outcome)             ~ "Active",
      grepl("^Attenuated", outcome)           ~ "Active, MYC-target-weighted",
      grepl("^Lost|^Reversed", outcome)       ~ "Lost",
      outcome == "Untestable after exclusion" ~ "Untestable",
      grepl("^Emerged", outcome)              ~ "Emerged",
      TRUE                                    ~ "Not significant on all genes"),
    active   = code %in% c("Active", "Active, MYC-target-weighted"),
    eligible = sig_before & agrees,
    cell     = paste0(set_keys[MYC_gene_set], "_", method))

cell_cols <- as.vector(outer(method_levels, names(myc_sets),
                             function(m, s) paste0(s, "_", m)))

summary_all <- cells |>
  dplyr::group_by(TF) |>
  dplyr::summarise(n_active = sum(active), n_eligible = sum(eligible),
                   direction_conflict = any(code == "Opposite direction"), .groups = "drop") |>
  dplyr::left_join(cells |>
                     dplyr::select(TF, cell, code) |>
                     tidyr::pivot_wider(names_from = cell, values_from = code), by = "TF") |>
  dplyr::left_join(tf_base, by = "TF") |>
  dplyr::arrange(dplyr::desc(n_active), ulm_p, msviper_p)

# TFs significant on all genes in at least one method, most active cells first (ties: ULM p, then msVIPER p)
summary_tf <- summary_all |>
  dplyr::filter(ulm_sig | mlm_sig | msviper_sig) |>
  dplyr::mutate(rank = dplyr::row_number()) |>
  dplyr::select(rank, TF, MYC_regulator, edge_to_MYC, consistent_with_MYC_up, direction,
                n_active, n_eligible, direction_conflict, dplyr::all_of(cell_cols),
                n_targets, ulm_score, ulm_p, mlm_score, mlm_p, msviper_nes, msviper_p)

cells |>
  dplyr::filter(!is.na(outcome)) |>
  dplyr::mutate(method = factor(method, levels = method_levels)) |>
  dplyr::group_by(MYC_gene_set, method) |>
  dplyr::summarise(significant_on_all_genes = sum(sig_before),
                   still_active             = sum(sig_before & grepl("^Survives|^Attenuated", outcome)),
                   of_which_MYC_target_weighted = sum(grepl("^Attenuated", outcome)),
                   lost                     = sum(sig_before & grepl("^Lost|^Reversed", outcome)),
                   too_few_targets_left     = sum(outcome == "Untestable after exclusion"),
                   emerged                  = sum(grepl("^Emerged", outcome)), .groups = "drop") |>
  knitr::kable(caption = "All TFs: outcome in each method, with each method's own direction")
All TFs: outcome in each method, with each method’s own direction
MYC_gene_set method significant_on_all_genes still_active of_which_MYC_target_weighted lost too_few_targets_left emerged
CollecTRI MYC targets ULM 109 41 3 53 15 21
CollecTRI MYC targets MLM 68 30 0 27 11 15
CollecTRI MYC targets msVIPER 309 285 0 24 0 0
Hallmark MYC_TARGETS_V1 + V2 ULM 109 85 17 22 2 9
Hallmark MYC_TARGETS_V1 + V2 MLM 68 54 0 14 0 8
Hallmark MYC_TARGETS_V1 + V2 msVIPER 309 289 0 20 0 1 
summary_tf |>
  dplyr::count(n_active, name = "TFs") |>
  dplyr::left_join(summary_tf |> dplyr::filter(MYC_regulator) |> dplyr::count(n_active, name = "MYC_regulators"),
                   by = "n_active") |>
  dplyr::mutate(MYC_regulators = dplyr::coalesce(MYC_regulators, 0L)) |>
  dplyr::arrange(dplyr::desc(n_active)) |>
  knitr::kable(col.names = c(sprintf("Cells still active (of %d)", length(cell_cols)), "TFs", "of which MYC regulators"),
               caption = sprintf("%d TFs significant on all genes in at least one method", nrow(summary_tf)))
416 TFs significant on all genes in at least one method
Cells still active (of 6) TFs of which MYC regulators
6 1 1
5 4 1
4 18 6
3 14 3
2 283 42
1 59 15
0 37 10 
lab_dir    <- function(d) if (nrow(d)) paste0(d$TF, ifelse(d$direction == "up in tFL", " (up)", " (down)"),
                                              collapse = ", ") else "none"
ulm_active <- vapply(names(myc_sets), function(key)
  sum(ex_long$MYC_gene_set == set_labels[[key]] & ex_long$ulm_sig &
        grepl("^Survives|^Attenuated", ex_long$ulm_verdict)), integer(1))
untest_sig <- vapply(names(myc_sets), function(key)
  sum(ex_long$MYC_gene_set == set_labels[[key]] & ex_long$ulm_verdict %in% "Untestable after exclusion"),
  integer(1))
top_tf     <- summary_tf |> dplyr::filter(n_active == max(n_active))
reg_top    <- head(summary_tf, params$n_top) |> dplyr::filter(MYC_regulator)
n_conflict <- sum(summary_tf$direction_conflict)

share_removed <- ex_long |>
  dplyr::filter(!is.na(n_targets_removed)) |>
  dplyr::mutate(share = n_targets_removed / n_targets) |>
  dplyr::group_by(MYC_gene_set) |>
  dplyr::summarise(all_TFs = stats::median(share), MYC_regulators = stats::median(share[MYC_regulator]),
                   .groups = "drop")
med_removed     <- setNames(share_removed$all_TFs,        set_keys[share_removed$MYC_gene_set])
med_removed_reg <- setNames(share_removed$MYC_regulators, set_keys[share_removed$MYC_gene_set])

At a glance

  • TFs significant on all genes: 416 in at least one method (ULM 109, MLM 68, msVIPER 309).
  • Still active in ULM: 85 of 109 without the Hallmark genes, and 41 without the CollecTRI MYC targets. Without the CollecTRI targets, 126 TFs have too few targets left to score, 15 of them significant on all genes.
  • Most overlapping (6 of 6 cells): GATA3 (down).
  • MYC regulators among the top 40: GATA3 (down), RUNX1 (down), MYC (up), E2F1 (up), E2F4 (up), FOXM1 (up), MYCN (up), SPI1 (down), IKZF1 (down), EWSR1 (up), AR (up), TFDP1 (up), TFDP2 (up).
  • Direction conflicts: 12 TFs are significant in opposite directions in different methods. Those cells don’t count.
tile_cols <- c("Active" = col_survive, "Active, MYC-target-weighted" = col_atten, "Lost" = col_lost,
               "Untestable" = ink_muted, "Opposite direction" = col_grid, "Not significant on all genes" = col_grid,
               "Emerged" = col_grid, "Not in network" = col_surface)
tile_txt  <- c("Active" = "S", "Active, MYC-target-weighted" = "A", "Lost" = "L", "Untestable" = "U",
               "Opposite direction" = "O", "Not significant on all genes" = "", "Emerged" = "E",
               "Not in network" = "-")

plot_tiles <- function(tf_rows, title, subtitle) {
  mark <- !all(tf_rows$MYC_regulator)   # star MYC regulators unless every row is one
  lab  <- tf_rows |>
    dplyr::mutate(label = sprintf("%s (%s) %d/%d%s", TF, sub(" in tFL", "", direction), n_active, n_eligible,
                                  ifelse(mark & MYC_regulator, " *", "")))
  df <- cells |>
    dplyr::inner_join(lab |> dplyr::select(TF, label), by = "TF") |>
    dplyr::mutate(label        = factor(label, levels = rev(lab$label)),
                  method       = factor(method, levels = method_levels),
                  MYC_gene_set = factor(paste("Without", MYC_gene_set), levels = paste("Without", set_labels)),
                  txt          = unname(tile_txt[code]),
                  txt_col      = ifelse(code %in% c("Active", "Active, MYC-target-weighted", "Lost", "Untestable"),
                                        "white", ink_second))

  ggplot(df, aes(method, label)) +
    geom_tile(aes(fill = code), colour = col_surface, linewidth = 0.8) +
    geom_text(aes(label = txt, colour = txt_col), size = 2.8) +
    facet_grid(~ MYC_gene_set) +
    scale_fill_manual(values = tile_cols,
                      breaks = c("Active", "Active, MYC-target-weighted", "Lost", "Untestable",
                                 "Not significant on all genes")) +
    scale_colour_identity() +
    scale_x_discrete(position = "top") +
    guides(fill = guide_legend(nrow = 2, byrow = TRUE)) +
    labs(title = title, subtitle = subtitle, x = NULL, y = NULL,
         caption = paste0("S still active (significant, same sign); A still active but leaning on MYC target genes (ULM only); ",
                          "L lost; U too few targets left.\nO significant on all genes in the other direction (does not count); ",
                          "E significant only after the change; - not in that network.\n",
                          "Label: TF (overall direction) cells that count / cells where it was significant on all genes ",
                          "in that direction.", if (mark) " * regulates MYC in CollecTRI." else "")) +
    theme_sens() +
    theme(panel.grid.major = element_blank(), strip.placement = "outside")
}
plot_tiles(head(summary_tf, params$n_top),
           "Which TFs are still active without MYC target genes?",
           sprintf("Top %d of %d TFs significant on all genes in at least one method, most cells first",
                   params$n_top, nrow(summary_tf)))

reg_rows <- summary_tf |> dplyr::filter(MYC_regulator)
plot_tiles(head(reg_rows, params$n_top),
           "Which MYC regulators are still active without MYC target genes?",
           sprintf("Top %d of %d MYC regulators significant on all genes in at least one method, most cells first",
                   min(params$n_top, nrow(reg_rows)), nrow(reg_rows)))

summary_tf |>
  dplyr::select(rank, TF, MYC_regulator, direction, n_active, n_eligible, direction_conflict,
                dplyr::all_of(cell_cols), ulm_score, ulm_p, mlm_score, mlm_p, msviper_nes, msviper_p) |>
  show_table("TFs significant on all genes in at least one method, ranked by cells still active without MYC target genes")

7.4 MYC regulators in detail

The same results for MYC’s regulators, for each gene set, drawn like the ranked figure under MYC regulators ranked by activity: the same regulators (significant in ULM or MLM on all genes), in the same order. Each panel shows one method’s score without the MYC target genes, filled if the TF is still active, with a grey dot at its score on all genes. The ULM verdicts against the random-target null are in the table below the figures.

myc_rows <- ex_long |> dplyr::filter(TF == "MYC")

ulm_groups <- c("Still active", "Still active, MYC-target-weighted", "Lost", "Too few targets left")
ex_long |>
  dplyr::filter(MYC_regulator, ulm_sig, TF != "MYC") |>
  dplyr::mutate(ulm_group = dplyr::case_when(
                  grepl("^Survives", ulm_verdict)             ~ ulm_groups[1],
                  grepl("^Attenuated", ulm_verdict)           ~ ulm_groups[2],
                  ulm_verdict == "Untestable after exclusion" ~ ulm_groups[4],
                  TRUE                                        ~ ulm_groups[3]),
                ulm_group    = factor(ulm_group, levels = ulm_groups),
                MYC_gene_set = factor(MYC_gene_set, levels = set_labels)) |>
  dplyr::arrange(MYC_gene_set, ulm_group, ulm_p) |>
  dplyr::group_by(MYC_gene_set, ulm_group, .drop = FALSE) |>
  dplyr::summarise(n   = dplyr::n(),
                   TFs = paste0(TF, ifelse(ulm_score > 0, " (up)", " (down)"), collapse = ", "),
                   .groups = "drop") |>
  knitr::kable(col.names = c("MYC target genes removed", "ULM", "n", "Regulators"),
               caption = "MYC regulators significant in ULM (MYC itself not counted)")
MYC regulators significant in ULM (MYC itself not counted)
MYC target genes removed ULM n Regulators
Hallmark MYC_TARGETS_V1 + V2 Still active 12 E2F4 (up), IKZF1 (down), GATA3 (down), NFE2L2 (up), KLF5 (up), EWSR1 (up), ETV3 (down), RUNX1 (down), NFATC2 (down), SPI1 (down), SP3 (up), E2F5 (up)
Hallmark MYC_TARGETS_V1 + V2 Still active, MYC-target-weighted 11 E2F1 (up), TFDP1 (up), FOXM1 (up), MYCN (up), TBP (up), SP1 (up), STAT3 (up), JUN (up), ESR1 (up), AR (up), HIF1A (up)
Hallmark MYC_TARGETS_V1 + V2 Lost 12 RB1 (down), EZH2 (up), MYBL2 (up), KAT5 (up), CREBBP (up), AHR (up), FOSB (up), MYB (up), NFKB (up), PATZ1 (up), GLI2 (up), CTCFL (up)
Hallmark MYC_TARGETS_V1 + V2 Too few targets left 1 TFDP2 (up)
CollecTRI MYC targets Still active 12 E2F4 (up), FOXM1 (up), MYCN (up), IKZF1 (down), GATA3 (down), NFE2L2 (up), AR (up), RUNX1 (down), RB1 (down), NFATC2 (down), SPI1 (down), E2F5 (up)
CollecTRI MYC targets Still active, MYC-target-weighted 2 E2F1 (up), STAT3 (up)
CollecTRI MYC targets Lost 20 TFDP1 (up), TBP (up), SP1 (up), KLF5 (up), JUN (up), EWSR1 (up), ESR1 (up), ETV3 (down), HIF1A (up), EZH2 (up), MYBL2 (up), KAT5 (up), SP3 (up), CREBBP (up), AHR (up), FOSB (up), MYB (up), NFKB (up), GLI2 (up), CTCFL (up)
CollecTRI MYC targets Too few targets left 2 TFDP2 (up), PATZ1 (up)

MYC’s own CollecTRI score is 15.4 on all genes and 12.1 without the Hallmark genes. Without its own CollecTRI targets it has nothing left to score.

tf_cols <- c(edge_cols, "not a MYC regulator" = ink_second)

# one method, drawn like the ranked figure: a lollipop from zero to the score without the MYC target
# genes (filled = still active: significant on all genes and without them, same sign), plus a grey
# dot at the score on all genes (hollow if not significant there)
without_panel <- function(d, before, after, sig_before, active, title, xlab, first = FALSE) {
  d <- d |>
    dplyr::mutate(b = .data[[before]], a = .data[[after]], sb = .data[[sig_before]], act = .data[[active]])
  p <- ggplot(d, aes(y = label)) +
    geom_vline(xintercept = 0, colour = col_context, linewidth = 0.3) +
    geom_segment(data = ~ dplyr::filter(.x, !is.na(a)), aes(x = 0, xend = a, yend = label),
                 colour = col_context, linewidth = 0.4) +
    geom_point(data = ~ dplyr::filter(.x, !is.na(b), sb), aes(x = b),
               shape = 21, size = 1.6, fill = col_context, colour = col_surface, stroke = 0.3) +
    geom_point(data = ~ dplyr::filter(.x, !is.na(b), !sb), aes(x = b),
               shape = 21, size = 1.4, fill = col_surface, colour = col_context, stroke = 0.6) +
    geom_point(data = ~ dplyr::filter(.x, !is.na(a), act), aes(x = a, fill = colour_group),
               shape = 21, size = 2.6, colour = col_surface, stroke = 0.5) +
    geom_point(data = ~ dplyr::filter(.x, !is.na(a), !act), aes(x = a, colour = colour_group),
               shape = 21, size = 2.4, fill = col_surface, stroke = 1) +
    scale_y_discrete(limits = levels(d$label)) +   # the same rows in every panel
    scale_fill_manual(values = tf_cols, breaks = names(tf_cols), guide = if (first) "legend" else "none") +
    scale_colour_manual(values = tf_cols, guide = "none") +
    labs(title = title, x = xlab, y = NULL) +
    theme_sens()
  if (first) p else p + theme(axis.text.y = element_blank())
}

plot_without <- function(rows, key, title) {
  fig <- ex_long |>
    dplyr::filter(MYC_gene_set == set_labels[[key]], TF %in% rows) |>
    dplyr::mutate(colour_group = dplyr::coalesce(edge_to_MYC, "not a MYC regulator"),
                  ulm_active   = still_active %in% TRUE,
                  mlm_active   = mlm_call %in% "Survives",
                  vip_active   = msviper_call %in% "Survives",
                  label        = sprintf("%s (%d/%d)", TF, as.integer(n_targets_removed), as.integer(n_targets)))
  fig$label <- factor(fig$label, levels = fig$label[order(fig$ulm_score)])

  without_panel(fig, "ulm_score", "ulm_score_after", "ulm_sig", "ulm_active", "ULM", "ULM score", first = TRUE) +
    without_panel(fig, "mlm_score", "mlm_score_after", "mlm_sig", "mlm_active", "MLM", "MLM score") +
    without_panel(fig, "msviper_nes", "msviper_nes_after", "msviper_sig", "vip_active",
                  "msVIPER", "msVIPER NES (bcellViper)") +
    patchwork::plot_layout(widths = c(1.3, 1, 1)) +
    patchwork::plot_annotation(
      title    = paste0(title, ", without the ", set_labels[[key]]),
      subtitle = "Score without the MYC target genes (positive = more active in tFL); grey dot: score on all genes",
      caption  = paste0("Filled points: still active (significant on all genes and without the MYC target genes, same sign); ",
                        "hollow: not. Hollow grey dot: not significant on all genes.\n",
                        "No coloured point: too few targets left without the MYC target genes, or (msVIPER) no bcellViper regulon. ",
                        "Rows sorted by ULM score on all genes.\n",
                        "Label: TF (MYC target genes among its CollecTRI targets / targets). ",
                        "ULM verdicts against the random-target null are in the tables and the without_MYC_targets sheet."),
      theme    = theme_sens())
}

ranked_regs <- as.character(plot_df$TF)   # the regulators of the ranked figure
plot_without(ranked_regs, "hallmark",  "MYC regulators significant in ULM or MLM")

plot_without(ranked_regs, "collectri", "MYC regulators significant in ULM or MLM")

ex_long |>
  dplyr::filter(MYC_regulator, ulm_sig) |>
  dplyr::select(TF, MYC_gene_set, edge_to_MYC, consistent_with_MYC_up, still_active, n_targets, n_targets_removed,
                ulm_score, ulm_score_after, ulm_p_after, retention, null_lo, null_hi, ulm_verdict,
                mlm_call, msviper_call) |>
  show_table("MYC regulators significant in ULM on all genes, for each MYC target gene set")

7.5 All TFs in detail

The same figures for all TFs: the 40 TFs with the largest |ULM score| among those significant in ULM or MLM on all genes, sorted by ULM score. MYC regulators are coloured by their edge to MYC, other TFs dark grey.

top_all <- tf_base |>
  dplyr::filter(ulm_sig | mlm_sig) |>
  dplyr::slice_max(abs(ulm_score), n = params$n_top, with_ties = FALSE) |>
  dplyr::pull(TF)

plot_without(top_all, "hallmark",  sprintf("Top %d TFs significant in ULM or MLM", params$n_top))

plot_without(top_all, "collectri", sprintf("Top %d TFs significant in ULM or MLM", params$n_top))

8 Export

readme <- tibble::tribble(
  ~column, ~meaning,
  "rank", "Position when sorted by ULM score, highest first (positive = more active in tFL).",
  "edge_to_MYC / mor_to_MYC", "Whether CollecTRI lists the TF as an activator (+1) or repressor (-1) of MYC.",
  "consistent_with_MYC_up", "TRUE if the TF's change would push MYC up in tFL: an activator more active in tFL, or a repressor less active. Based on the ULM score sign.",
  "n_targets", "The TF's CollecTRI targets present in the signature (baseMean > 4.99).",
  "shared_with_MYC_targets", "Fraction of the TF's targets that are also MYC targets in CollecTRI. High values mean its activity partly reflects MYC's own downstream program.",
  "ulm_* / mlm_*", "decoupleR scores (t-values), nominal p and BH-FDR across all scored TFs, on the full network and signature.",
  "ulm_sig / mlm_sig", paste0(if (params$p_adjust == "BH") "BH-FDR" else "nominal p", " < ", params$p_cut, "."),
  "*_without_MYC", "The same scores with the MYC gene removed from the signature, so no TF uses MYC as a target.",
  "n_hallmark_MYC_targets", "How many of the TF's targets in the signature are Hallmark MYC target genes (MYC_TARGETS_V1 + V2).",
  "still_active_without_hallmark_MYC / still_active_without_collectri_MYC", "TFs significant in ULM only: TRUE if still significant in ULM with the same sign after removing that MYC target gene set (Attenuated included). Details in the without_MYC_targets sheet.",
  "n_cells_active_without_MYC_targets / n_cells_eligible", "Of the 6 summary cells (ULM, MLM and msVIPER, each without the Hallmark and without the CollecTRI MYC target genes): how many the TF is still active in, in its overall direction, and in how many it was significant on all genes in that direction. See the MYC_target_summary sheet.",
  "msviper_*", "msVIPER with shadow on bcellViper (NES, p, BH-FDR), if the TF has a bcellViper regulon. MYC_in_bcellviper_regulon: MYC is among its bcellViper targets.",
  "rho_MYC_samples / p_rho_samples", "Spearman correlation of per-sample TF activity level (MYC gene excluded) with MYC mRNA level, across all samples (samples with 0 MYC reads left out if exclude_myc_dropouts = TRUE).",
  "rho_MYC_tFL / p_rho_tFL", "The same level correlation across tFL samples only.",
  "rho_MYC_within_patients / p_rho_within", paste0("Spearman correlation of the tFL - FL change in activity with the change in MYC mRNA (",
                                                   n_pairs_used, " patients; those with a MYC dropout left out if exclude_myc_dropouts = TRUE)."),
  "rho_MYC_within_all_patients", "The same correlation using all 11 patients, including MYC dropouts.",
  "per_sample_support", "'tracks MYC as expected': within-patient rho has the sign of the edge (activator +, repressor -) and p < 0.05.",
  "change_tFL_vs_FL / p_change", "Paired FL -> tFL change in per-sample activity (activity ~ patient + FL/tFL).",
  "change_adj_for_MYC / p_change_adj / p_MYC_term", "The same change with MYC mRNA added to the model, and the p-value of the MYC term.",
  "deseq_*", "The TF's own DESeq2 result (tFL vs FL).",
  "MYC_target_summary sheet", "One row per TF significant on all genes in at least one method, ranked by n_active (ties: ULM p, then msVIPER p). direction: from ULM if significant there, otherwise MLM, otherwise msVIPER. direction_conflict: significant in the opposite direction in another method.",
  "hallmark_* / collectri_* (summary sheet)", "Outcome per cell: Active (still significant, same sign) / Active, MYC-target-weighted (ULM, below the random-target null) / Lost / Untestable (too few targets left) / Opposite direction (significant on all genes against the TF's overall direction; does not count) / Emerged (significant only after the change) / Not significant on all genes / Not in network.",
  "without_MYC_targets sheet", "One row per TF and MYC target gene set, for all TFs. *_after: scores without that gene set (FDR across all TFs in that run). The *_sig columns in this sheet and the summary sheet use all TFs.",
  "n_targets_removed", "How many of the TF's CollecTRI targets in the signature are in that MYC target gene set.",
  "retention / null_lo / null_med / null_hi / p_null", "ULM only: score kept after / before the exclusion, and the 2.5th / 50th / 97.5th percentile of the score kept when the same number of random targets is removed. p_null: share of random removals that keep as little or less.",
  "ulm_verdict", "Survives: robust / Survives (no MYC target genes among its targets) / Attenuated (still significant, but its MYC target genes carried more signal than random targets) / Lost: MYC-target-driven / Lost: no worse than random target loss / Lost: background shift only / Reversed / Untestable after exclusion / Emerged / Not significant.",
  "still_active", "TFs significant in ULM only: TRUE if still significant in ULM with the same sign after the exclusion (Attenuated included).",
  "mlm_call / msviper_call", "Survives (same sign, still significant) / Lost / Lost (direction reversed) / Untestable after exclusion / Emerged / Not significant; empty if the TF is not in that network.",
  "MYC_target_genes sheet", "Every gene in either MYC target set: Hallmark V1 / V2 membership, CollecTRI MYC target and its sign (collectri_MYC_mor), E2F_TARGETS / G2M_CHECKPOINT membership, whether it is in each signature, and its DESeq2 result."
)

run_info <- tibble::tibble(
  item  = c("date", "DESeq2 input", "p_cut", "p_adjust", "minsize", "CollecTRI TFs with MYC as a target", "scored",
            "signature genes", "samples with 0 MYC reads", "patients left out of the per-sample link",
            "Hallmark MYC target genes", "Hallmark MYC target genes: total / in DESeq2 file / in decoupleR signature / in msVIPER signature",
            "CollecTRI MYC targets: total / in decoupleR signature / in msVIPER signature",
            "n_perm / seed / n_top", "decoupleR", "viper", "R"),
  value = c(format(Sys.time()), deseq_file, params$p_cut, params$p_adjust, params$minsize,
            nrow(myc_regulators), sum(myc_regulators$scored), nrow(mat_stat),
            paste(sub("_normCounts$", "", myc_zero), collapse = ", "), paste(drop_patients, collapse = ", "),
            "MSigDB Hallmark MYC_TARGETS_V1 + MYC_TARGETS_V2 (human), gsea-msigdb.org, copied 2026-09-29",
            paste(nrow(myc_set_tbl), length(myc_set_genes), sets_tbl$in_decoupleR_signature[1],
                  sets_tbl$in_msVIPER_signature[1], sep = " / "),
            paste(nrow(collectri_myc), sets_tbl$in_decoupleR_signature[2], sets_tbl$in_msVIPER_signature[2], sep = " / "),
            paste(params$n_perm, params$seed, params$n_top, sep = " / "),
            as.character(packageVersion("decoupleR")), as.character(packageVersion("viper")), R.version.string))

# main sheet: the headline MYC-target results sit next to the MYC-gene check
mt_cols <- ex_long |>
  dplyr::filter(MYC_gene_set == set_labels[["hallmark"]]) |>
  dplyr::select(TF, n_hallmark_MYC_targets = n_targets_removed, still_active_without_hallmark_MYC = still_active) |>
  dplyr::left_join(ex_long |>
                     dplyr::filter(MYC_gene_set == set_labels[["collectri"]]) |>
                     dplyr::select(TF, still_active_without_collectri_MYC = still_active), by = "TF") |>
  dplyr::left_join(summary_all |> dplyr::select(TF, n_cells_active_without_MYC_targets = n_active,
                                                n_cells_eligible = n_eligible), by = "TF")
myc_sheet <- myc_table |>
  dplyr::left_join(mt_cols, by = "TF") |>
  dplyr::relocate(dplyr::all_of(setdiff(names(mt_cols), "TF")), .after = mlm_p_without_MYC)

openxlsx::write.xlsx(
  list(README              = readme,
       MYC_regulators      = myc_sheet,
       MYC_target_summary  = summary_tf,
       without_MYC_targets = ex_long,
       MYC_target_genes    = myc_gene_tbl,
       per_sample_activity = per_sample,
       not_scored          = myc_regulators |> dplyr::filter(!scored),
       run_info            = run_info),
  file = out_xlsx, overwrite = TRUE, firstRow = TRUE
)
cat("Wrote", out_xlsx, "\n")
## Wrote C:/Users/User/Desktop/MM Work/Main Work/fl_tFL analysis R/FL_TFL_MYC_regulators.xlsx

9 How to read the results

  • The activity score is not MYC’s. A TF’s score comes from all its targets. MYC is one target among many, so the scores barely change when MYC is removed (see At a glance).
  • A CollecTRI edge is literature evidence that the TF regulates MYC in some cell type. It is not proof that it does so in B cells.
  • Direction. An activator that is more active in tFL, or a repressor that is less active, is consistent with MYC going up. An inconsistent TF would, if anything, push MYC down.
  • Correlation with MYC mRNA does not show which way regulation runs. A TF can track MYC because it drives MYC, because MYC drives its targets (check shared_with_MYC_targets), or because both follow proliferation or tumour content.
  • MYC can rise without upstream TFs. Translocations and amplifications activate MYC directly, and MYC protein stability matters too. A weak link to upstream TFs is expected in patients with such events.
  • Small numbers. With 9 pairs, only |rho| above about 0.67 is significant. MYC mRNA rises in most pairs, so the MYC-adjusted model has limited power to separate the FL→tFL change from the MYC change.
  • MYC is lowly expressed (tens of reads per sample). The two samples with 0 MYC reads are the two smallest libraries, so they are treated as dropouts; set exclude_myc_dropouts: false to keep them.
  • p-values are nominal by default (FDR is in the table), and ULM and MLM can disagree for TFs with overlapping regulons.
  • Without MYC target genes. A TF that stays active has signal outside MYC’s own program. A TF that is lost is ambiguous: its call may have been MYC’s program, or it may co-regulate MYC’s target genes and have too little signal elsewhere. The random-target null tells the two apart only partly (Lost: MYC-target-driven vs Lost: no worse than random target loss).
  • Expect up-in-tFL TFs to fall below the null. The MYC target genes are up in tFL (see the gene-group table), so removing them costs more than removing random targets. What matters is whether the remaining targets still give a significant signal: Attenuated (yes) or Lost: MYC-target-driven (no).
  • The two gene sets are tests of different strength. The Hallmark set is MYC’s core program: the median TF has 0% of its targets in it (4% for MYC regulators), so most calls barely change. The CollecTRI set is every literature target of MYC: it removes a median 33% of a TF’s targets (37% for MYC regulators) and leaves 126 TFs with too few to score. The two sets share 56 genes in the signature.
  • The MYC target genes overlap E2F/G2M. Removing them also takes out part of the proliferation program: 50 of the 319 E2F/G2M genes with the Hallmark set (e.g. the MCM genes and PCNA), 80 with the CollecTRI set.
  • Overlap counts favour TFs in both networks. A TF only in bcellViper, or significant in a single method, can count in at most 2 of the 6 cells. Read n_active next to n_eligible.
  • Direction conflicts. msVIPER on bcellViper sometimes points the other way from ULM and MLM on CollecTRI (12 TFs). Those cells are marked O in the figures and don’t count.

10 Session info

sessionInfo()
## R version 4.4.0 (2024-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22000)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Europe/Budapest
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] openxlsx_4.2.9      patchwork_1.3.2     ggplot2_4.0.3      
##  [4] tidyr_1.3.2         dplyr_1.2.1         bcellViper_1.42.0  
##  [7] viper_1.40.0        Biobase_2.66.0      BiocGenerics_0.52.0
## [10] decoupleR_2.12.0   
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6        xfun_0.61           bslib_0.12.0       
##  [4] htmlwidgets_1.6.4   lattice_0.22-6      crosstalk_1.2.2    
##  [7] vctrs_0.7.3         tools_4.4.0         generics_0.1.4     
## [10] parallel_4.4.0      tibble_3.3.1        proxy_0.4-29       
## [13] pkgconfig_2.0.3     Matrix_1.7-6        KernSmooth_2.23-22 
## [16] data.table_1.17.4   RColorBrewer_1.1-3  S7_0.2.1           
## [19] lifecycle_1.0.5     compiler_4.4.0      farver_2.1.2       
## [22] stringr_1.6.0       mixtools_2.0.0.1    codetools_0.2-20   
## [25] htmltools_0.5.8.1   class_7.3-22        sass_0.4.10        
## [28] yaml_2.3.12         plotly_4.12.1       pillar_1.11.1      
## [31] jquerylib_0.1.4     MASS_7.3-60.2       DT_0.34.0          
## [34] BiocParallel_1.40.2 cachem_1.1.0        nlme_3.1-164       
## [37] parallelly_1.48.0   zip_2.3.3           tidyselect_1.2.1   
## [40] digest_0.6.35       stringi_1.8.9       purrr_1.2.2        
## [43] kernlab_0.9-33      labeling_0.4.3      splines_4.4.0      
## [46] fastmap_1.2.0       grid_4.4.0          cli_3.6.6          
## [49] magrittr_2.0.5      survival_3.5-8      e1071_1.7-17       
## [52] withr_3.0.3         scales_1.4.0        segmented_2.2-2    
## [55] rmarkdown_2.32      httr_1.4.9          otel_0.2.0         
## [58] evaluate_1.0.5      knitr_1.52          viridisLite_0.4.3  
## [61] mgcv_1.9-1          rlang_1.3.0         Rcpp_1.0.12        
## [64] glue_1.8.0          rstudioapi_0.19.0   jsonlite_2.0.0     
## [67] R6_2.6.1