FL vs tFL TF activity: proliferation sensitivity analysis

Do the TF activity calls survive removing Hallmark E2F_TARGETS ∪ G2M_CHECKPOINT genes?

1 About this analysis

This file reruns the TF activity tests from Decoupler analysis.Rmd twice on the same DESeq2 signature: once with all genes (before), and once with the Hallmark E2F_TARGETS ∪ G2M_CHECKPOINT genes removed (after). It then compares the two runs TF by TF, to see which TF calls hold up once the proliferation genes are gone.

• Input: DESeq2_mraz.tsv (11 FL vs 11 tFL; set with deseq_file in the YAML params). The signature is the DESeq2 Wald stat, filtered as in Decoupler analysis.Rmd. Positive = higher in tFL (see the sign check under Input signature).
• Tests (unchanged):
  • CollecTRI with decoupleR ULM and MLM (minsize = 5).
  • bcellViper with classic msVIPER with shadow (pleiotropy correction).
• Exclusion: the proliferation genes are removed from the signature itself, so they count as neither TF targets nor background genes. The networks are not edited; each method matches its regulons against the smaller signature.
• Still active: significant before, and still significant after exclusion with the same sign. For ULM, a size-matched random-target null (see ULM: size-matched random-target null) splits these into robust and Attenuated: an Attenuated TF’s proliferation targets carried more signal than random targets would. Both count as still active, as in Decoupler MYC regulators.Rmd. MLM and msVIPER have no null, because it would take hours to compute.
• Summary across methods: each TF gets 3 cells, ULM, MLM and msVIPER. A cell counts if the TF is still active there and points in the TF’s overall direction. TFs are ranked by how many cells count (see Which TFs are still active).
• Significance: nominal p < 0.05. You can change this in the YAML params.

2 Libraries and paths

suppressPackageStartupMessages({
  library(decoupleR)
  library(viper)
  library(bcellViper)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(ggrepel)
  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: three categorical slots (validated all-pairs for scatter), greys for context
col_survive <- "#2a78d6"
col_atten   <- "#eb6834"
col_lost    <- "#1baf7a"
col_context <- "#c3c2b7"
col_grid    <- "#e1e0d9"
col_surface <- "#fcfcfb"
ink_primary <- "#0b0b0b"
ink_second  <- "#52514e"
ink_muted   <- "#898781"

group_cols <- c("Survives"                       = col_survive,
                "Attenuated (still significant)" = col_atten,
                "Lost / reversed"                = col_lost,
                "Not significant before"         = col_context)

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",
          legend.position       = "top",
          legend.justification  = "left",
          legend.title          = element_blank(),
          legend.text           = element_text(colour = ink_second))
}

3 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)

Sign check. Proliferation markers should be higher in tFL, and their stat should be positive.

norm_fl  <- grep("^FL_.*_normCounts$",  colnames(deseq_unfiltered), value = TRUE)
norm_tfl <- grep("^tFL_.*_normCounts$", colnames(deseq_unfiltered), value = TRUE)

sign_check <- deseq_unfiltered |>
  dplyr::filter(gene_name %in% c("MKI67", "TOP2A", "CDK1", "MYC")) |>
  dplyr::mutate(mean_norm_FL  = rowMeans(dplyr::pick(dplyr::all_of(norm_fl))),
                mean_norm_tFL = rowMeans(dplyr::pick(dplyr::all_of(norm_tfl)))) |>
  dplyr::transmute(gene_name, stat, mean_norm_FL, mean_norm_tFL,
                   higher_in = ifelse(mean_norm_tFL > mean_norm_FL, "tFL", "FL"))
knitr::kable(sign_check, digits = 2,
             caption = sprintf("%d FL and %d tFL samples", length(norm_fl), length(norm_tfl)))

11 FL and 11 tFL samples

gene_name	stat	mean_norm_FL	mean_norm_tFL	higher_in
CDK1	7.07	45.32	195.16	tFL
TOP2A	6.77	198.41	766.19	tFL
MKI67	5.56	143.44	484.68	tFL
MYC	3.67	20.65	87.95	tFL

if (!all((sign_check$stat > 0) == (sign_check$higher_in == "tFL")))
  warning("DESeq2 stat does not follow tFL > FL for these genes; check the contrast direction.")

4 Proliferation gene set

The lists below are the MSigDB Hallmark gene sets E2F_TARGETS and G2M_CHECKPOINT (200 genes each), human release v2026.1.Hs, copied from the gene set pages on gsea-msigdb.org on 2026-09-24. Gene names are current HGNC symbols. Your DESeq2 annotation still uses some older names (e.g. H2AFZ, PAPD7), so each renamed gene is also matched on the original symbol that MSigDB gives for it.

This list replaces the 327-symbol vector in Dano_Data/prolif_residual_screen.R, which lacked KIF23 and SS18 and had two clone IDs (AC027237.1, AC091021.1) in their place.

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

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_e2f)) == 200, length(unique(hallmark_g2m)) == 200)

# current HGNC symbol -> original MSigDB source symbol, for renamed genes
msigdb_alias <- c(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)
}

prolif_tbl <- tibble::tibble(msigdb_symbol = union(hallmark_e2f, hallmark_g2m)) |>
  dplyr::mutate(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),
                in_decoupler_signature = !is.na(data_symbol) & data_symbol %in% rownames(mat_stat),
                in_msviper_signature   = !is.na(data_symbol) & data_symbol %in% names(sign_C)) |>
  dplyr::left_join(deseq_lookup |> dplyr::select(gene_name, deseq_baseMean, deseq_stat),
                   by = c("data_symbol" = "gene_name")) |>
  dplyr::arrange(msigdb_symbol)

prolif_genes <- sort(unique(stats::na.omit(prolif_tbl$data_symbol)))

tibble::tibble(
  gene_set = c("E2F_TARGETS", "G2M_CHECKPOINT", "Union", "  found in DESeq2 file",
               "  in decoupleR signature (baseMean > 4.99)", "  in msVIPER signature (baseMean > 10)"),
  genes = c(length(hallmark_e2f), length(hallmark_g2m), nrow(prolif_tbl), length(prolif_genes),
            sum(prolif_tbl$in_decoupler_signature), sum(prolif_tbl$in_msviper_signature))) |>
  knitr::kable()

gene_set	genes
E2F_TARGETS	200
G2M_CHECKPOINT	200
Union	327
found in DESeq2 file	323
in decoupleR signature (baseMean > 4.99)	319
in msVIPER signature (baseMean > 10)	316

renamed <- prolif_tbl |> dplyr::filter(!is.na(data_symbol), data_symbol != msigdb_symbol)
cat("Matched on the older symbol:",
    paste(renamed$msigdb_symbol, renamed$data_symbol, sep = " -> ", collapse = ", "), "\n")

## Matched on the older symbol: H2AX -> H2AFX, H2AZ1 -> H2AFZ, H2AZ2 -> H2AFV, H2BC12 -> HIST1H2BK, PRP4K -> PRPF4B, TENT4A -> PAPD7

cat("Not in the DESeq2 file:",
    paste(prolif_tbl$msigdb_symbol[is.na(prolif_tbl$data_symbol)], collapse = ", "), "\n")

## Not in the DESeq2 file: HOXC10, NME1, RBM14, SPAG5

stat_df <- tibble::tibble(gene = rownames(mat_stat), stat = mat_stat[, 1]) |>
  dplyr::mutate(group = ifelse(gene %in% prolif_genes, "E2F/G2M genes", "All other genes"))
stat_means <- stat_df |>
  dplyr::group_by(group) |>
  dplyr::summarise(mean = mean(stat), n = dplyr::n(), .groups = "drop")

ggplot(stat_df, aes(stat, colour = group)) +
  geom_vline(xintercept = 0, colour = col_context, linewidth = 0.3) +
  geom_density(linewidth = 0.8, key_glyph = "path") +
  geom_vline(data = stat_means, aes(xintercept = mean, colour = group), linewidth = 0.5,
             show.legend = FALSE) +
  scale_colour_manual(values = c("All other genes" = ink_muted, "E2F/G2M genes" = col_survive)) +
  labs(title = "Proliferation genes are shifted up in tFL",
       subtitle = paste(sprintf("%s: n = %d, mean stat %+.2f", stat_means$group, stat_means$n,
                                stat_means$mean), collapse = "   |   "),
       x = "DESeq2 stat (positive = higher in tFL)", y = "Density",
       caption = "Vertical lines: group means. decoupleR signature (baseMean > 4.99).") +
  theme_sens()

5 Signatures after exclusion

mat_stat_ex <- mat_stat[!rownames(mat_stat) %in% prolif_genes, , drop = FALSE]
sign_C_ex   <- sign_C[!names(sign_C) %in% prolif_genes]

tibble::tibble(signature      = c("decoupleR (CollecTRI)", "msVIPER (bcellViper)"),
               genes_before   = c(nrow(mat_stat), length(sign_C)),
               prolif_removed = c(nrow(mat_stat) - nrow(mat_stat_ex), length(sign_C) - length(sign_C_ex)),
               genes_after    = c(nrow(mat_stat_ex), length(sign_C_ex))) |>
  knitr::kable()

signature	genes_before	prolif_removed	genes_after
decoupleR (CollecTRI)	19627	319	19308
msVIPER (bcellViper)	17093	316	16777

6 CollecTRI: decoupleR ULM + MLM, before and after

The decouple() call is identical to section 2a of Decoupler analysis.Rmd. It runs once on each signature.

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

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
  )
}

res_collectri_before <- run_collectri(mat_stat)
res_collectri_after  <- run_collectri(mat_stat_ex)

dplyr::bind_rows(before = res_collectri_before, after = res_collectri_after, .id = "run") |>
  dplyr::count(statistic, run) |>
  tidyr::pivot_wider(names_from = run, values_from = n) |>
  knitr::kable(caption = "TFs scored (regulons with >= minsize targets in the signature)")

TFs scored (regulons with >= minsize targets in the signature)

statistic	after	before
mlm	662	700
ulm	662	700

Regulon composition: how many of each TF’s CollecTRI targets (within the signature) are proliferation genes.

collectri_comp <- net_collectri |>
  dplyr::filter(target %in% rownames(mat_stat)) |>
  dplyr::distinct(source, target) |>
  dplyr::group_by(source) |>
  dplyr::summarise(n_targets_before = dplyr::n(),
                   n_prolif_targets = sum(target %in% prolif_genes),
                   prolif_targets   = paste(sort(target[target %in% prolif_genes]), collapse = ", "),
                   .groups = "drop") |>
  dplyr::mutate(n_targets_after = n_targets_before - n_prolif_targets,
                prolif_frac     = n_prolif_targets / n_targets_before)

7 ULM: size-matched random-target null

When a regulon loses targets, its ULM t-value falls even if those targets were ordinary, because t grows with regulon size. So “significant before, weaker after” is not evidence on its own.

For each TF with m proliferation targets out of k, the null removes m random targets of that TF instead, 1,000 times. Everything else matches the real exclusion. In particular, every proliferation gene that is not a target of this TF is removed from the background in each permutation too, so the only difference is which targets are removed. If the retained score falls below the null’s 2.5th percentile, the TF’s proliferation targets carried more signal than its average target.

The null uses a closed-form ULM. It computes the same t-value as decoupleR::run_ulm, which is checked against the decoupleR results before the null is run.

y     <- mat_stat[, 1]
genes <- rownames(mat_stat)
is_p  <- genes %in% prolif_genes

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)))
}

t_cf_before <- setNames(ulm_closed(y, W), tfs)
t_cf_after  <- setNames(ulm_closed(y[!is_p], W[!is_p, , drop = FALSE]), tfs)
t_cf_after[Matrix::colSums(W[!is_p, , drop = FALSE] != 0) < params$minsize] <- NA

ulm_dc_before <- res_collectri_before |> dplyr::filter(statistic == "ulm")
ulm_dc_after  <- res_collectri_after  |> dplyr::filter(statistic == "ulm")

validation <- tibble::tibble(
  run         = c("before", "after"),
  tfs_decoupleR   = c(nrow(ulm_dc_before), nrow(ulm_dc_after)),
  tfs_closed_form = c(sum(!is.na(t_cf_before)), sum(!is.na(t_cf_after))),
  same_tfs    = c(setequal(ulm_dc_before$source, names(t_cf_before)[!is.na(t_cf_before)]),
                  setequal(ulm_dc_after$source,  names(t_cf_after)[!is.na(t_cf_after)])),
  max_abs_diff = c(max(abs(ulm_dc_before$score - t_cf_before[ulm_dc_before$source])),
                   max(abs(ulm_dc_after$score  - t_cf_after[ulm_dc_after$source]))))
knitr::kable(validation, digits = 16, caption = "Closed-form ULM vs decoupleR::run_ulm")

Closed-form ULM vs decoupleR::run_ulm

run	tfs_decoupleR	tfs_closed_form	same_tfs	max_abs_diff
before	700	700	TRUE	1.6e-14
after	662	662	TRUE	6.2e-15

if (!all(validation$same_tfs) || any(validation$max_abs_diff > 1e-6))
  stop("Closed-form ULM does not reproduce decoupleR; the random-target null would not be valid.")

set.seed(params$seed)
n_all <- length(y); sy_all <- sum(y);       syy_all <- sum(y^2)
n_p   <- sum(is_p); sy_p   <- sum(y[is_p]); syy_p   <- sum(y[is_p]^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]
  pj  <- is_p[idx]
  k   <- length(idx); m <- sum(pj)
  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 proliferation gene that is NOT a target of this TF removed
  nA   <- n_all   - (n_p - m)
  syA  <- sy_all  - (sy_p  - sum(yj[pj]))
  syyA <- syy_all - (syy_p - sum(yj[pj]^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[pj]), sum(yj[pj]^2), sum(w[pj]), sum(w[pj]^2), sum(w[pj] * yj[pj]))
  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
}

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

# the real exclusion inside null_one must equal the decoupleR 'after' run
chk <- ulm_null |>
  dplyr::inner_join(ulm_dc_after |> dplyr::select(source, score), by = "source")
stopifnot(max(abs(chk$t_after - chk$score)) < 1e-6)

8 bcellViper: msVIPER with shadow, before and after

Same call as section 3b of Decoupler analysis.Rmd. Both runs use the same regulon set (TFs expressed at baseMean > 10), so a TF whose own gene is a proliferation gene is still scored after exclusion. msVIPER uses its default minimum of 25 targets.

data(bcellViper, package = "bcellViper", envir = environment())
regulon_filtered <- regulon[intersect(names(regulon), deg$gene_name)]

mrs_before <- viper::msviper(sign_C,    regulon_filtered, pleiotropy = TRUE, verbose = FALSE)
mrs_after  <- viper::msviper(sign_C_ex, regulon_filtered, pleiotropy = TRUE, verbose = FALSE)

viper_es <- function(mrs) {
  tibble::tibble(source  = names(mrs$es$nes),
                 score   = as.numeric(mrs$es$nes),
                 p_value = as.numeric(mrs$es$p.value),
                 size    = as.numeric(mrs$es$size))
}
viper_before <- viper_es(mrs_before)
viper_after  <- viper_es(mrs_after)

viper_comp <- tibble::tibble(
  source           = names(regulon_filtered),
  n_targets_before = vapply(regulon_filtered, function(r) sum(names(r$tfmode) %in% names(sign_C)), integer(1)),
  n_prolif_targets = vapply(regulon_filtered, function(r) {
    g <- names(r$tfmode); sum(g %in% names(sign_C) & g %in% prolif_genes)
  }, integer(1)),
  prolif_targets   = vapply(regulon_filtered, function(r) {
    g <- names(r$tfmode); paste(sort(g[g %in% names(sign_C) & g %in% prolif_genes]), collapse = ", ")
  }, character(1))) |>
  dplyr::mutate(n_targets_after = n_targets_before - n_prolif_targets,
                prolif_frac     = n_prolif_targets / n_targets_before)

tibble::tibble(run = c("before", "after"),
               regulons_scored = c(nrow(viper_before), nrow(viper_after))) |>
  knitr::kable(caption = sprintf("bcellViper regulons expressed: %d; median share of proliferation targets: %.1f%%",
                                 length(regulon_filtered), 100 * median(viper_comp$prolif_frac)))

bcellViper regulons expressed: 529; median share of proliferation targets: 4.3%

run	regulons_scored
before	529
after	529

9 Which TFs are still active

Each TF gets 3 cells: ULM, MLM and msVIPER. A cell counts if the TF is still active there, meaning significant after exclusion with the same sign as on all genes (Attenuated included), 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. This is the same summary as in Decoupler MYC regulators.Rmd.

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

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")
}

# one row per TF: <method>_score, <method>_p, <method>_fdr (FDR across all TFs in the run)
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))
}
viper_wide <- function(v, suffix = "") {
  v |>
    dplyr::transmute(source, msviper_nes = score, msviper_p = p_value,
                     msviper_fdr = stats::p.adjust(p_value, "BH")) |>
    dplyr::rename_with(~ paste0(.x, suffix), -source)
}

# TFs whose canonical program IS proliferation. Removing E2F/G2M targets removes the genes
# that define their activity, so a loss does not mean the call was an artifact. Flagged only.
prolif_program_tfs <- c(paste0("E2F", 1:8), "TFDP1", "TFDP2", "TFDP3", "FOXM1", "MYBL1", "MYBL2",
                        "LIN9", "LIN54", "RB1", "RBL1", "RBL2", "HCFC1", "MYC", "MYCN", "MYCL")

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)
}

# one row per TF, on all genes and without the proliferation genes, every method side by side
tf_tab <- to_wide(res_collectri_before) |>
  dplyr::full_join(viper_wide(viper_before), by = "source") |>
  dplyr::left_join(to_wide(res_collectri_after, "_after"), by = "source") |>
  dplyr::left_join(viper_wide(viper_after, "_after"), by = "source") |>
  dplyr::left_join(collectri_comp |> dplyr::select(source, n_targets = n_targets_before,
                                                   n_targets_removed = n_prolif_targets, prolif_targets),
                   by = "source") |>
  dplyr::left_join(ulm_null |> dplyr::select(source, retention, null_lo, null_med, null_hi, p_null), by = "source") |>
  dplyr::left_join(viper_comp |> dplyr::select(source, bcellviper_n_targets = n_targets_before,
                                               bcellviper_n_removed = n_prolif_targets), by = "source") |>
  dplyr::mutate(
    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"),
    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 proliferation targets)",
      ulm_call == "Survives" & retention >= null_lo   ~ "Survives: robust",
      ulm_call == "Survives"                          ~ "Attenuated: still significant, proliferation-weighted",
      ulm_call == "Lost (direction reversed)"         ~ "Reversed: significant in the opposite direction after exclusion",
      ulm_call == "Lost" & n_targets_removed == 0     ~ "Lost: background shift only (no proliferation targets)",
      ulm_call == "Lost" & retention < null_lo        ~ "Lost: proliferation-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)),
    prolif_program_tf = source %in% prolif_program_tfs,
    tf_in_prolif_set  = source %in% prolif_genes) |>
  dplyr::left_join(deseq_lookup, by = c("source" = "gene_name")) |>
  dplyr::rename(TF = source) |>
  dplyr::select(TF, prolif_program_tf, tf_in_prolif_set, 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,
                n_targets_removed, prolif_targets, 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,
                bcellviper_n_targets, bcellviper_n_removed, dplyr::starts_with("deseq_"))

# one row per TF and method; a cell counts if still active in the TF's overall direction
method_levels <- c("ULM", "MLM", "msVIPER")
cells <- dplyr::bind_rows(
  tf_tab |> dplyr::transmute(TF, method = "ULM",     outcome = ulm_verdict,
                             sig_before = ulm_sig,     sign_before = sign(ulm_score)),
  tf_tab |> dplyr::transmute(TF, method = "MLM",     outcome = mlm_call,
                             sig_before = mlm_sig,     sign_before = sign(mlm_score)),
  tf_tab |> dplyr::transmute(TF, method = "msVIPER", outcome = msviper_call,
                             sig_before = msviper_sig, sign_before = sign(msviper_nes))) |>
  dplyr::left_join(tf_tab |> 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, proliferation-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, proliferation-weighted"),
    eligible = sig_before & agrees,
    cell     = paste0("prolif_", method))

cell_cols <- paste0("prolif_", method_levels)

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_tab, 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, prolif_program_tf, direction, n_active, n_eligible, direction_conflict,
                dplyr::all_of(cell_cols), n_targets, n_targets_removed,
                ulm_score, ulm_p, mlm_score, mlm_p, msviper_nes, msviper_p)

9.1 Outcome counts

cells |>
  dplyr::filter(!is.na(outcome)) |>
  dplyr::mutate(method = factor(method, levels = method_levels)) |>
  dplyr::group_by(method) |>
  dplyr::summarise(significant_on_all_genes        = sum(sig_before),
                   still_active                    = sum(sig_before & grepl("^Survives|^Attenuated", outcome)),
                   of_which_proliferation_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

method	significant_on_all_genes	still_active	of_which_proliferation_weighted	lost	too_few_targets_left	emerged
ULM	109	50	17	50	9	13
MLM	68	48	0	17	3	14
msVIPER	309	286	0	23	0	0

tf_list <- function(tf, direction, max_n = 30) {
  lab <- paste0(tf, ifelse(direction == "up in tFL", " (up)", " (down)"))
  paste0(paste(head(lab, max_n), collapse = ", "),
         if (length(lab) > max_n) sprintf(", ... (+%d more)", length(lab) - max_n) else "")
}
ulm_levels <- c("Survives: robust", "Survives (no proliferation targets)",
                "Attenuated: still significant, proliferation-weighted",
                "Lost: proliferation-driven", "Lost: no worse than random target loss",
                "Lost: background shift only (no proliferation targets)",
                "Reversed: significant in the opposite direction after exclusion",
                "Untestable after exclusion")
call_levels <- c("Survives", "Lost", "Lost (direction reversed)", "Untestable after exclusion")

dplyr::bind_rows(
  "ULM (CollecTRI)"      = tf_tab |> dplyr::filter(ulm_sig) |>
                             dplyr::transmute(TF, outcome = factor(ulm_verdict, levels = ulm_levels), p = ulm_p,
                                              dir = ifelse(ulm_score > 0, "up in tFL", "down in tFL")),
  "MLM (CollecTRI)"      = tf_tab |> dplyr::filter(mlm_sig) |>
                             dplyr::transmute(TF, outcome = factor(mlm_call, levels = call_levels), p = mlm_p,
                                              dir = ifelse(mlm_score > 0, "up in tFL", "down in tFL")),
  "msVIPER (bcellViper)" = tf_tab |> dplyr::filter(msviper_sig) |>
                             dplyr::transmute(TF, outcome = factor(msviper_call, levels = call_levels), p = msviper_p,
                                              dir = ifelse(msviper_nes > 0, "up in tFL", "down in tFL")),
  .id = "method") |>
  dplyr::arrange(method, outcome, p) |>
  dplyr::group_by(method, outcome) |>
  dplyr::summarise(n = dplyr::n(), TFs = tf_list(TF, dir), .groups = "drop") |>
  knitr::kable(caption = "Outcome of every TF significant on all genes, per method (most significant first)")

Outcome of every TF significant on all genes, per method (most significant first)

method	outcome	n	TFs
MLM (CollecTRI)	Untestable after exclusion	3	TOX3 (down), MTF2 (up), ARID1A (down)
MLM (CollecTRI)	Survives	48	MYC (up), E2F4 (up), E2F1 (up), FOXO3 (down), MYCN (up), ASCL1 (up), GRHL3 (down), ATF6 (up), NANOG (down), GATA3 (down), AP1 (down), RUNX1 (down), POU4F1 (down), NFIL3 (down), TP63 (down), FOXC1 (down), IRF1 (down), ZKSCAN7 (up), NPM1 (down), MSX1 (down), SRSF2 (up), SOX6 (up), IKZF1 (down), NONO (up), NKX2-2 (down), STAT5B (up), FOXE1 (down), BHLHE41 (up), HOXA9 (down), IRF2 (up), … (+18 more)
MLM (CollecTRI)	Lost	17	OLIG2 (up), KLF10 (down), TFDP1 (up), AR (up), SMARCC1 (up), TCF3 (down), HES1 (up), E2F2 (up), ZNF143 (up), PGR (down), TLX1 (down), SSRP1 (down), BATF (down), FOXM1 (up), ZFPM1 (up), RBPJ (down), ELK4 (up)
ULM (CollecTRI)	Survives: robust	25	FOXO3 (down), IKZF1 (down), TFAM (up), GATA3 (down), ATF6 (up), HSF1 (up), NFE2L2 (up), DDIT3 (up), SREBF2 (up), FOXP3 (down), RUNX1 (down), HIF1A (up), TBX21 (down), PPARGC1A (up), NFYA (up), ID4 (down), SREBF1 (up), CEBPD (up), THRB (up), HSF4 (up), TBX1 (down), ELK1 (up), KLF7 (down), NFYC (up), FOXN1 (up)
ULM (CollecTRI)	Survives (no proliferation targets)	8	BATF (down), NFIL3 (down), MSX1 (down), ZFPM2 (down), HOPX (up), CREBZF (up), NRL (down), MTA3 (down)
ULM (CollecTRI)	Attenuated: still significant, proliferation-weighted	17	MYC (up), E2F1 (up), E2F4 (up), E2F3 (up), TFDP1 (up), FOXM1 (up), MYCN (up), ZNF143 (up), SRSF2 (up), OLIG2 (up), IRF1 (down), TBP (up), SP1 (up), JUN (up), NKX2-1 (up), AR (up), NRF1 (up)
ULM (CollecTRI)	Lost: proliferation-driven	28	E2F2 (up), ZHX2 (down), KLF10 (down), STAT3 (up), MSC (down), KCNIP3 (up), KLF5 (up), EWSR1 (up), NCOA3 (up), ESR1 (up), SMARCA1 (up), HDAC1 (down), EZH2 (up), MYBL2 (up), KAT5 (up), ATF1 (up), SP3 (up), NKX6-1 (up), CREBBP (up), AHR (up), HMGA2 (up), KMT2B (up), MYB (up), HR (up), POU2F1 (up), NFKB (up), DNMT1 (up), GLI2 (up)
ULM (CollecTRI)	Lost: no worse than random target loss	22	ETV3 (down), RB1 (down), ZNF331 (down), NFATC2 (down), HMGB2 (up), GTF3A (up), NR5A2 (up), SPI1 (down), ASCL1 (up), NONO (up), STAT6 (down), TBX15 (up), HES6 (up), FOSB (up), BHLHE41 (up), KLF3 (down), PATZ1 (up), NFIB (down), NKX2-2 (down), TP63 (down), E2F5 (up), CTCFL (up)
ULM (CollecTRI)	Untestable after exclusion	9	HCFC1 (up), ARID3A (up), STOX1 (up), TFDP2 (up), TFDP3 (down), E2F7 (down), TRIM28 (down), TOX3 (down), CTBP1 (down)
msVIPER (bcellViper)	Survives	286	HNRNPAB (up), MRPL28 (up), PHF1 (down), TOP2A (up), FOXM1 (up), PRKDC (up), ILF2 (up), NR1D2 (down), NOLC1 (up), FOXJ2 (down), TFEB (down), RBL2 (down), ZFP36L2 (down), STAT5A (down), ZNF862 (down), SP100 (down), ZFP36L1 (down), HMGA1 (up), MYBL2 (up), HDGF (up), ZNF264 (down), ZNF510 (down), TFAP4 (up), ZNF266 (down), ZBTB20 (down), PLAGL1 (down), MEF2A (down), CREBBP (down), STAT3 (down), TAF5 (up), … (+256 more)
msVIPER (bcellViper)	Lost	23	KLF10 (up), AEBP1 (down), BATF (down), ARNT2 (down), NFIB (down), DRAP1 (down), NR2F2 (down), MEOX1 (down), CNOT8 (down), TRAFD1 (down), HIF1A (down), ZNF467 (down), ATF2 (down), ZFX (down), ZNF24 (down), SPEN (down), CEBPD (down), FOXC1 (up), SNAI2 (down), MLXIP (down), THOC1 (up), ZSCAN12 (down), ZNF354A (down)

summary_tf |>
  dplyr::count(n_active, name = "TFs") |>
  dplyr::left_join(summary_tf |> dplyr::filter(prolif_program_tf) |> dplyr::count(n_active, name = "prolif_program_TFs"),
                   by = "n_active") |>
  dplyr::mutate(prolif_program_TFs = dplyr::coalesce(prolif_program_TFs, 0L)) |>
  dplyr::arrange(dplyr::desc(n_active)) |>
  knitr::kable(col.names = c(sprintf("Cells still active (of %d)", length(cell_cols)), "TFs",
                             "of which proliferation-program TFs"),
               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 3)	TFs	of which proliferation-program TFs
3	5	1
2	20	5
1	321	8
0	70	4

lab_dir    <- function(d) if (nrow(d)) paste0(d$TF, ifelse(d$direction == "up in tFL", " (up)", " (down)"),
                                              collapse = ", ") else "none"
n_still    <- function(sig, outcome) sum(tf_tab[[sig]] & grepl("^Survives|^Attenuated", tf_tab[[outcome]]))
top_tf     <- summary_tf |> dplyr::filter(n_active == max(n_active))
prog_act   <- summary_tf |> dplyr::filter(prolif_program_tf, n_active > 0)
n_conflict <- sum(summary_tf$direction_conflict)

At a glance

• TFs significant on all genes: 416 in at least one method (ULM 109, MLM 68, msVIPER 309).
• Still active without the E2F/G2M genes: ULM 50 (17 of them proliferation-weighted), MLM 48, msVIPER 286. In ULM, 9 TFs have too few targets left to score.
• Most overlapping (3 of 3 cells): MYC (up), FOXO3 (down), GATA3 (down), IRF1 (down), RUNX1 (down).
• Proliferation-program TFs still active in at least one cell: MYC (up), E2F1 (up), E2F4 (up), E2F3 (up), FOXM1 (up), MYCN (up), TFDP1 (up), E2F2 (up), HCFC1 (up), TFDP2 (up), MYBL2 (up), E2F6 (up), RBL2 (down), RBL1 (down).
• Direction conflicts: 12 TFs are significant in opposite directions in different methods. Those cells don’t count.

tile_cols <- c("Active" = col_survive, "Active, proliferation-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, proliferation-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$prolif_program_tf)   # star proliferation-program TFs 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 & prolif_program_tf, " *", "")))
  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),
                  set     = "Without E2F_TARGETS + G2M_CHECKPOINT genes",
                  txt     = unname(tile_txt[code]),
                  txt_col = ifelse(code %in% c("Active", "Active, proliferation-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(~ set) +
    scale_fill_manual(values = tile_cols,
                      breaks = c("Active", "Active, proliferation-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 proliferation genes (ULM only);\n",
                          "L lost; U too few targets left; O significant on all genes in the other direction (does not count);\n",
                          "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) "\n* proliferation-program TF." 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 the E2F/G2M 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)))

prog_rows <- summary_tf |> dplyr::filter(prolif_program_tf)
plot_tiles(prog_rows,
           "Proliferation-program TFs without the E2F/G2M genes",
           sprintf("All %d proliferation-program TFs significant on all genes in at least one method", nrow(prog_rows)))

summary_tf |>
  dplyr::select(rank, TF, prolif_program_tf, 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 the E2F/G2M genes")

9.2 Top TFs in detail

The 40 TFs with the largest |ULM score| without the E2F/G2M genes, sorted by that score, drawn like the figures in Decoupler MYC regulators.Rmd. Each panel shows one method’s score without the E2F/G2M genes, filled if the TF is still active, with a grey dot at its score on all genes. A TF that becomes significant in ULM only without the E2F/G2M genes can make the list; its ULM grey dot is hollow. Proliferation-program TFs are blue, other TFs dark grey. The ULM verdicts against the random-target null are in the table below.

tf_cols <- c("proliferation-program TF" = col_survive, "other TF" = ink_second)

# TFs ranked by |ULM score without the E2F/G2M genes|, among those significant in ULM on all genes
# or without the E2F/G2M genes; the per-sample plots use the same order
ranked_after <- tf_tab |>
  dplyr::filter(ulm_sig | is_sig(ulm_p_after)) |>
  dplyr::arrange(dplyr::desc(abs(ulm_score_after))) |>   # untestable after exclusion (NA) last
  dplyr::pull(TF)

fig_df <- tf_tab |>
  dplyr::filter(TF %in% head(ranked_after, params$n_top)) |>
  dplyr::mutate(colour_group = ifelse(prolif_program_tf, "proliferation-program TF", "other TF"),
                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_df$label <- factor(fig_df$label, levels = fig_df$label[order(fig_df$ulm_score_after, na.last = FALSE)])

# one method, drawn like the ranked figures of the MYC report: a lollipop from zero to the score
# without the E2F/G2M 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())
}

without_panel(fig_df, "ulm_score", "ulm_score_after", "ulm_sig", "ulm_active", "ULM", "ULM score", first = TRUE) +
  without_panel(fig_df, "mlm_score", "mlm_score_after", "mlm_sig", "mlm_active", "MLM", "MLM score") +
  without_panel(fig_df, "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    = sprintf("Top %d TFs by ULM score without the E2F/G2M genes", params$n_top),
    subtitle = "Score without the E2F_TARGETS + G2M_CHECKPOINT 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 E2F/G2M genes, same sign); ",
                      "hollow: not. Hollow grey dot: not significant on all genes.\n",
                      "In the ULM panel, a hollow grey dot marks a TF that is significant only without the E2F/G2M genes. ",
                      "Rows sorted by ULM score without the E2F/G2M genes.\n",
                      "No coloured point: too few targets left without the E2F/G2M genes, or (msVIPER) no bcellViper regulon. ",
                      "Label: TF (E2F/G2M genes among its CollecTRI targets / targets).\n",
                      "ULM verdicts against the random-target null are in the table below and the without_prolif_genes sheet."),
    theme    = theme_sens())

tf_tab |>
  dplyr::filter(ulm_sig) |>
  dplyr::select(TF, prolif_program_tf, direction, still_active, ulm_verdict, ulm_score, ulm_score_after,
                ulm_p, ulm_p_after, retention, null_lo, null_hi, n_targets, n_targets_removed,
                mlm_call, msviper_call) |>
  show_table("ULM: TFs significant on all genes")

9.3 Before vs after, per method

plot_group <- function(outcome, sig_before) {
  dplyr::case_when(grepl("^Survives", outcome)   ~ "Survives",
                   grepl("^Attenuated", outcome) ~ "Attenuated (still significant)",
                   sig_before                    ~ "Lost / reversed",
                   TRUE                          ~ "Not significant before")
}

# one method's columns from tf_tab, in the shape the scatter plot needs
method_tab <- function(score, score_after, sig, outcome) {
  tf_tab |>
    dplyr::filter(!is.na(.data[[score]])) |>
    dplyr::transmute(source = TF, score_before = .data[[score]], score_after = .data[[score_after]],
                     sig_before = .data[[sig]], call = .data[[outcome]])
}

plot_before_after <- function(tab, method_label, score_label) {
  df <- tab |>
    dplyr::filter(!is.na(score_after)) |>
    dplyr::mutate(group = factor(plot_group(call, sig_before), levels = names(group_cols)))
  lab <- dplyr::bind_rows(
    df |> dplyr::filter(sig_before) |> dplyr::slice_max(abs(score_before), n = params$n_label),
    df |> dplyr::filter(call == "Emerged after exclusion") |> dplyr::slice_max(abs(score_after), n = 5)) |>
    dplyr::distinct(source, .keep_all = TRUE)
  n_untestable <- sum(tab$call == "Untestable after exclusion", na.rm = TRUE)

  ggplot(df, aes(score_before, score_after)) +
    geom_hline(yintercept = 0, colour = col_context, linewidth = 0.3) +
    geom_vline(xintercept = 0, colour = col_context, linewidth = 0.3) +
    geom_abline(slope = 1, intercept = 0, colour = ink_muted, linewidth = 0.3) +
    geom_point(data = ~ dplyr::filter(.x, group == "Not significant before"),
               aes(fill = group), shape = 21, size = 1.6, colour = col_surface, stroke = 0.3) +
    geom_point(data = ~ dplyr::filter(.x, group != "Not significant before"),
               aes(fill = group), shape = 21, size = 2.6, colour = col_surface, stroke = 0.5) +
    ggrepel::geom_text_repel(data = lab, aes(label = source), size = 3, colour = ink_second,
                             segment.colour = col_context, segment.size = 0.3, max.overlaps = Inf,
                             min.segment.length = 0, box.padding = 0.4, point.padding = 0.25,
                             seed = params$seed) +
    scale_fill_manual(values = group_cols, breaks = names(group_cols)) +
    coord_equal() +
    labs(title = paste0(method_label, ": before vs after removing E2F/G2M genes"),
         x = paste(score_label, "before exclusion (all genes)"),
         y = paste(score_label, "after exclusion"),
         caption = paste0("Diagonal: unchanged score. Positive = more active in tFL. Labels: top ",
                          params$n_label, " by |score| before, plus up to 5 TFs that emerged after exclusion.\n",
                          n_untestable, " TF(s) significant before cannot be scored after exclusion and are not shown.")) +
    theme_sens()
}

plot_before_after(method_tab("ulm_score", "ulm_score_after", "ulm_sig", "ulm_verdict"), "ULM (CollecTRI)", "ULM score")

plot_before_after(method_tab("mlm_score", "mlm_score_after", "mlm_sig", "mlm_call"), "MLM (CollecTRI)", "MLM score")

plot_before_after(method_tab("msviper_nes", "msviper_nes_after", "msviper_sig", "msviper_call"), "msVIPER (bcellViper)", "NES")

TFs that emerged after exclusion (not significant on all genes, significant without the E2F/G2M genes):

dplyr::bind_rows(
  "ULM"     = tf_tab |> dplyr::filter(ulm_verdict %in% "Emerged after exclusion") |>
                dplyr::transmute(TF, score = ulm_score, score_after = ulm_score_after, p = ulm_p, p_after = ulm_p_after),
  "MLM"     = tf_tab |> dplyr::filter(mlm_call %in% "Emerged after exclusion") |>
                dplyr::transmute(TF, score = mlm_score, score_after = mlm_score_after, p = mlm_p, p_after = mlm_p_after),
  "msVIPER" = tf_tab |> dplyr::filter(msviper_call %in% "Emerged after exclusion") |>
                dplyr::transmute(TF, score = msviper_nes, score_after = msviper_nes_after, p = msviper_p,
                                 p_after = msviper_p_after),
  .id = "method") |>
  dplyr::left_join(tf_tab |> dplyr::select(TF, n_targets, n_targets_removed), by = "TF") |>
  show_table()

10 Per-sample activity and expression

One figure per TF, FL to tFL, with a line per patient:

• Left: the TF’s activity in each sample: ULM on that sample’s log2(normalized count + 1), each gene centred on its mean across the 22 samples. All genes are kept.
• Right: the TF’s own mRNA, log2(normalized count + 1). A line’s rise is that patient’s log2 fold change. CollecTRI complexes such as NFKB and AP1 have no single gene, so their mRNA panel is empty.
• Subtitle: the TF’s ULM score on all genes and without the E2F/G2M genes, its outcome without them in each method, and its DESeq2 log2 fold change (from the input file, shrunken) and padj.

Choose the TFs with per_sample_tfs: in the YAML header: a number such as 40 (the top N by |ULM score| without the E2F/G2M genes, as in Top TFs in detail), "all" with the quotes (every TF significant in ULM on all genes or without the E2F/G2M genes, in the same order), or a list of names such as ["GATA3", "RUNX1"].

ps_nc   <- grep("normCounts$", colnames(deseq_unfiltered), value = TRUE)
ps_rows <- deseq[!duplicated(deseq$gene_name), ]
ps_rows <- ps_rows[!is.na(ps_rows$stat), ]
ps_X    <- log2(as.matrix(ps_rows[, ps_nc]) + 1)
rownames(ps_X) <- ps_rows$gene_name
ps_X    <- ps_X[apply(ps_X, 1, sd) > 0 & rowSums(ps_X > 0) >= 3, ]

ps_samples <- tibble::tibble(sample  = colnames(ps_X),
                             group   = ifelse(grepl("^tFL_", colnames(ps_X)), "tFL", "FL"),
                             patient = sub("^t?FL_(\\d+)_normCounts$", "\\1", colnames(ps_X)))

# activity per sample: ULM on each sample's centred log2 expression
ps_act <- decoupleR::run_ulm(ps_X - rowMeans(ps_X), network = net_collectri,
                             .source = "source", .target = "target", .mor = "mor",
                             minsize = params$minsize) |>
  dplyr::select(TF = source, sample = condition, activity = score) |>
  dplyr::left_join(ps_samples, by = "sample")

# which TFs to plot: a number = top N by |ULM score without the E2F/G2M genes| (ranked_after, as in
# Top TFs in detail), "all" = every one of them, or TF names (any TF scored per sample)
sel    <- params$per_sample_tfs
ps_tfs <- if (is.numeric(sel)) head(ranked_after, sel) else if (identical(sel, "all")) ranked_after else
            intersect(sel, unique(ps_act$TF))

ps_cols <- c(FL = "#898781", tFL = "#1baf7a")   # the FL / tFL colours of the MYC report

# a TF's outcome without the E2F/G2M genes, in a few words
short_call <- function(x) dplyr::case_when(
  is.na(x)                          ~ "n/a",
  grepl("^Survives", x)             ~ "still active",
  grepl("^Attenuated", x)           ~ "still active (proliferation-weighted)",
  x == "Untestable after exclusion" ~ "too few targets left",
  grepl("^Emerged", x)              ~ "newly significant",
  grepl("^Lost|^Reversed", x)       ~ "lost",
  TRUE                              ~ "not significant")   # neither on all genes nor without them

# FL to tFL, one line per patient
paired_panel <- function(df, y, title) {
  w <- tidyr::pivot_wider(df, id_cols = patient, names_from = group, values_from = dplyr::all_of(y))
  ggplot(df, aes(group, .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 = ps_cols) +
    scale_x_discrete(expand = expansion(add = 0.4)) +
    labs(title = title,
         subtitle = sprintf("Higher in tFL in %d of %d patients", sum(w$tFL > w$FL, na.rm = TRUE),
                            sum(!is.na(w$tFL - w$FL))),
         x = NULL, y = NULL) +
    theme_sens()
}

plot_per_sample <- function(tf) {
  d   <- ps_act |> dplyr::filter(TF == tf)
  row <- match(tf, deseq_unfiltered$gene_name)   # first row per gene, as in deseq_lookup
  de  <- deseq_lookup[match(tf, deseq_lookup$gene_name), ]
  tt  <- tf_tab[match(tf, tf_tab$TF), ]

  p_act  <- paired_panel(d, "activity", "Activity per sample (ULM t-value)")
  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 = ps_nc, tf_mrna = log2(as.numeric(unlist(deseq_unfiltered[row, ps_nc])) + 1)) |>
      dplyr::left_join(ps_samples, by = "sample") |>
      paired_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))

  scores   <- sprintf("ULM score %+.1f on all genes, %s without the E2F/G2M genes.", tt$ulm_score,
                      if (is.na(tt$ulm_score_after)) "n/a" else sprintf("%+.1f", tt$ulm_score_after))
  outcomes <- sprintf("Without the E2F/G2M genes: ULM %s, MLM %s, msVIPER %s.", short_call(tt$ulm_verdict),
                      short_call(tt$mlm_call), short_call(tt$msviper_call))

  p_act + p_mrna +
    patchwork::plot_annotation(
      title    = if (isTRUE(tt$prolif_program_tf)) paste(tf, "(proliferation-program TF)") else tf,
      subtitle = paste(c(scores, strwrap(outcomes, 110), paste("mRNA:", mrna_txt)), collapse = "\n"),   # wrapped to fit
      theme    = theme_sens())
}

for (tf in ps_tfs) print(plot_per_sample(tf))

11 Export

readme <- tibble::tribble(
  ~column, ~meaning,
  "prolif_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", "The TF's overall direction: from ULM if significant there, otherwise MLM, otherwise msVIPER.",
  "prolif_ULM / prolif_MLM / prolif_msVIPER", "Outcome per cell: Active (still significant, same sign) / Active, proliferation-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 exclusion) / Not significant on all genes / Not in network.",
  "n_active / n_eligible", "Cells the TF is still active in, in its overall direction, and cells where it was significant on all genes in that direction.",
  "direction_conflict", "Significant on all genes in the opposite direction in another method.",
  "without_prolif_genes sheet", "One row per TF (CollecTRI and bcellViper), every method side by side.",
  "ulm_* / mlm_* / msviper_*", "Scores on all genes (ULM, MLM: t-value; msVIPER: NES), nominal p and BH-FDR across all TFs scored in that run.",
  "*_after", "The same scores with the E2F_TARGETS / G2M_CHECKPOINT genes removed from the signature.",
  "*_sig", paste0(if (params$p_adjust == "BH") "BH-FDR" else "nominal p", " < ", params$p_cut, ", across all TFs."),
  "n_targets / n_targets_removed / prolif_targets", "The TF's CollecTRI targets in the signature, how many of them are E2F/G2M genes, and which.",
  "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 (E2F/G2M genes that are not targets of the TF are removed in every permutation). p_null: share of random removals that keep as little or less.",
  "ulm_verdict", "Survives: robust / Survives (no proliferation targets) / Attenuated (still significant, but its proliferation targets carried more signal than random targets) / Lost: proliferation-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.",
  "prolif_program_tf", "TF whose canonical program is proliferation (E2F, TFDP, FOXM1, MYBL, MuvB, RB family, HCFC1, MYC family). Losing these after exclusion is expected and does not mean the call was an artifact.",
  "tf_in_prolif_set", "The TF's own gene is an E2F/G2M gene. This does not affect its activity score, which comes from its targets.",
  "bcellviper_n_targets / bcellviper_n_removed", "The TF's bcellViper targets in the msVIPER signature, and how many of them are E2F/G2M genes.",
  "deseq_*", "The TF's own DESeq2 result (tFL vs FL).",
  "prolif_genes sheet", "The E2F_TARGETS / G2M_CHECKPOINT genes, the symbol matched in the data, and whether each is in the signatures.",
  "per_sample_activity sheet", "Every TF's activity in every sample (ULM t-value on log2(normalized count + 1), each gene centred across the 22 samples, all genes kept), with patient and FL/tFL."
)

run_info <- tibble::tibble(
  item  = c("date", "DESeq2 input", "p_cut", "p_adjust", "minsize (decoupleR)", "minsize (msVIPER)", "n_perm", "seed", "n_top",
            "decoupleR signature genes before / after", "msVIPER signature genes before / after",
            "proliferation gene source", "decoupleR", "viper", "bcellViper", "R"),
  value = c(format(Sys.time()), deseq_file, params$p_cut, params$p_adjust, params$minsize, 25, params$n_perm, params$seed,
            params$n_top, paste(nrow(mat_stat), "/", nrow(mat_stat_ex)), paste(length(sign_C), "/", length(sign_C_ex)),
            "MSigDB Hallmark E2F_TARGETS + G2M_CHECKPOINT, v2026.1.Hs",
            as.character(packageVersion("decoupleR")), as.character(packageVersion("viper")),
            as.character(packageVersion("bcellViper")), R.version.string))

openxlsx::write.xlsx(
  list(README               = readme,
       prolif_summary       = summary_tf,
       without_prolif_genes = tf_tab,
       prolif_genes         = prolif_tbl,
       per_sample_activity  = ps_act |>
                                dplyr::mutate(sample = sub("_normCounts$", "", sample)) |>
                                dplyr::select(TF, sample, patient, group, activity),
       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_TF_activity_prolif_sensitivity.xlsx

12 How to read the results

• Why ULM has a null. A smaller regulon means less power, so comparing p-values before and after overstates how much signal was lost. The random-target null takes that shrinkage into account. Below null means the TF’s proliferation targets carried more signal than its average target.
• Expect up-in-tFL TFs with proliferation targets to fall below the null. Proliferation genes are strongly up in tFL (see the density plot), so removing them costs more than removing random targets. What matters is whether the remaining targets still give a significant signal: Attenuated (yes, still active) or Lost: proliferation-driven (no).
• Proliferation-program TFs (prolif_program_tf): their canonical targets are the removed genes. If one is lost, the test cannot tell “confounded by proliferation” apart from “regulates proliferation”. If one stays active, its signal reaches beyond the cell-cycle program.
• TFs with no proliferation targets can still shift slightly, because removing the proliferation genes also changes the background.
• Emerged TFs were masked by the proliferation genes, often through the background. Treat them as leads rather than findings.
• Overlap counts favour TFs in both networks. A TF only in bcellViper, or significant in a single method, can count in at most 1 of the 3 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.
• p-values: ULM and MLM treat genes as independent, which regulon targets are not, so nominal p-values are anti-conservative. Rank by score, and use p_adjust: "BH" for a stricter cut.
• msVIPER: bcellViper regulons contain few E2F/G2M genes (median 4.3%), so smaller changes are expected there than with CollecTRI.

13 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     ggrepel_0.9.8      
##  [4] ggplot2_4.0.3       tidyr_1.3.2         dplyr_1.2.1        
##  [7] bcellViper_1.42.0   viper_1.40.0        Biobase_2.66.0     
## [10] BiocGenerics_0.52.0 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] rlang_1.3.0         Rcpp_1.0.12         glue_1.8.0         
## [64] rstudioapi_0.19.0   jsonlite_2.0.0      R6_2.6.1