Figure S5: EBMF vs matched-RQVI factor comparison across level2 clusters

workflowr

• Summary
• Checks
• Past versions

Last updated: 2026-07-24

Checks: 7 0

Knit directory: immgenT-GP-analysis/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.7.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(1)

The command set.seed(1) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: b138063

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version b138063. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .DS_Store
    Ignored:    analysis/.DS_Store
    Ignored:    analysis/.Rhistory
    Ignored:    analysis/assets/.DS_Store
    Ignored:    code/.DS_Store
    Ignored:    data
    Ignored:    experiments/.DS_Store
    Ignored:    figures/.DS_Store
    Ignored:    figures/final-selected/.DS_Store
    Ignored:    figures/final-selected/bits/.DS_Store
    Ignored:    figures/final-selected/bits/Figure 1/.DS_Store
    Ignored:    figures/final-selected/bits/Figure 2/.DS_Store
    Ignored:    figures/final-selected/bits/Figure 3/.DS_Store
    Ignored:    figures/final-selected/bits/Figure 4/.DS_Store
    Ignored:    figures/final-selected/bits/Figure 6/.DS_Store
    Ignored:    figures/final-selected/bits/Figure 7/.DS_Store
    Ignored:    figures/final-selected/bits/Figure S1/.DS_Store
    Ignored:    figures/final-selected/bits/Figure S2/.DS_Store
    Ignored:    figures/final-selected/bits/Figure S3/.DS_Store
    Ignored:    figures/final-selected/bits/Figure S6/.DS_Store
    Ignored:    figures/final-selected/bits/Figure S7/.DS_Store
    Ignored:    figures/generated/.DS_Store
    Ignored:    figures/generated/Figure 1/.DS_Store
    Ignored:    figures/generated/Figure 2/.DS_Store
    Ignored:    tmp/

Untracked files:
    Untracked:  .claude/
    Untracked:  code/other/topic_flashier_20250212.R
    Untracked:  code/other/topic_wrapper_20250215_alldata_backfit.sh
    Untracked:  experiments/active_metric_comparison/
    Untracked:  experiments/figure6b_sparsity_ordered/
    Untracked:  experiments/gene_umap_gp_space/
    Untracked:  experiments/protein_threshold_gp_loading_diff/
    Untracked:  experiments/surface_protein_activation_scatter/
    Untracked:  figures/generated/Figure 1/1B.pdf
    Untracked:  figures/generated/Figure 3/3A.pdf
    Untracked:  figures/generated/Figure 3/3B.pdf
    Untracked:  figures/generated/Figure 3/3H.pdf
    Untracked:  figures/generated/Figure 3/3I.pdf
    Untracked:  figures/generated/Figure 3/3J.pdf
    Untracked:  figures/generated/Figure 3/3K.pdf
    Untracked:  figures/generated/Figure 3/3L.pdf
    Untracked:  figures/generated/Figure 3/3M.pdf
    Untracked:  figures/generated/Figure 4/4f.pdf
    Untracked:  figures/generated/Figure 4/4g.pdf
    Untracked:  figures/generated/Figure 5/
    Untracked:  figures/generated/Figure 6/6a.pdf
    Untracked:  figures/generated/Figure 7/
    Untracked:  log/2026-07-13-surface-protein-activation-scatter.md
    Untracked:  log/2026-07-15-protein-threshold-gp-loading-diff-tables.md
    Untracked:  plan/2026-07-13-active-cell-gene-metrics-exploration.md
    Untracked:  plan/2026-07-13-surface-protein-activation-scatter.md
    Untracked:  plan/2026-07-13_figure6b_triangular_order_plan.md
    Untracked:  plan/2026-07-13_figure6b_wide_gp_columns_formal_plan.md
    Untracked:  plan/2026-07-14_figure6b_heatmap_typography_plan.md
    Untracked:  plan/2026-07-15-protein-threshold-gp-loading-diff-tables.md
    Untracked:  tables/

Unstaged changes:
    Modified:   analysis/Methods_FlashierFit.Rmd
    Modified:   figures/generated/Figure 1/1A.pdf
    Modified:   figures/generated/Figure 1/1C.pdf
    Modified:   figures/generated/Figure 1/1D.pdf
    Deleted:    figures/generated/Figure 1/1E.pdf
    Deleted:    figures/generated/Figure 1/1F.pdf
    Deleted:    figures/generated/Figure 1/1G.pdf
    Deleted:    figures/generated/Figure 1/1H.pdf
    Deleted:    figures/generated/Figure 1/1I.pdf
    Deleted:    figures/generated/Figure 1/hist_active_cells_per_GP.pdf
    Deleted:    figures/generated/Figure 1/hist_active_genes_prop_per_GP.pdf
    Deleted:    figures/generated/Figure 1/scatter_e_active_cells_vs_genes.pdf
    Modified:   figures/generated/Figure 2/2A.pdf
    Modified:   figures/generated/Figure 2/2B.pdf
    Modified:   figures/generated/Figure 2/2C.pdf
    Modified:   figures/generated/Figure 2/2D.pdf
    Modified:   figures/generated/Figure 2/2E.pdf
    Modified:   figures/generated/Figure 2/2F.pdf
    Deleted:    figures/generated/Figure 2/2G.pdf
    Deleted:    figures/generated/Figure 2/2H.pdf
    Deleted:    figures/generated/Figure 2/2I.pdf
    Deleted:    figures/generated/Figure 2/2J.pdf
    Deleted:    figures/generated/Figure 2/2K.pdf
    Deleted:    figures/generated/Figure 2/2L.pdf
    Deleted:    figures/generated/Figure 2/2M.pdf
    Modified:   figures/generated/Figure 3/3c.pdf
    Modified:   figures/generated/Figure 3/3d.pdf
    Modified:   figures/generated/Figure 3/3e.pdf
    Modified:   figures/generated/Figure 3/3f.pdf
    Modified:   figures/generated/Figure 3/3g.pdf
    Deleted:    figures/generated/Figure 4/4a.pdf
    Deleted:    figures/generated/Figure 4/4b.pdf
    Modified:   figures/generated/Figure 4/4c.pdf
    Modified:   figures/generated/Figure 4/4d.pdf
    Modified:   figures/generated/Figure 4/4e.pdf
    Modified:   figures/generated/Figure S1/S1A.pdf
    Modified:   figures/generated/Figure S1/S1B.pdf
    Modified:   figures/generated/Figure S1/S1C.pdf
    Modified:   figures/generated/Figure S1/S1D.pdf
    Modified:   figures/generated/Figure S1/S1E.pdf
    Modified:   figures/generated/Figure S3/s3c.pdf
    Modified:   figures/generated/Figure S3/s3d.pdf
    Modified:   figures/generated/Figure S3/s3e.pdf
    Modified:   figures/generated/Figure S3/s3f.pdf
    Modified:   figures/generated/Figure S3/s3g.pdf
    Modified:   figures/generated/Figure S3/s3h.pdf

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/FigureS5.Rmd) and HTML (docs/FigureS5.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	b138063	Ziang Zhang	2026-07-24	Reformat Figure S5 page: lead with the figure, concise methods, link
html	9398c72	Ziang Zhang	2026-07-23	Publish Figure S5 workflowr page
Rmd	b9f4f58	Ziang Zhang	2026-07-23	Add Figure S5: EBMF vs matched-RQVI level2-cluster comparison

---

Figure S5 compares the cluster-level activity of the 200 EBMF gene programs (our flashier fit) with 200 matched RQVI programs from our collaborator (Tianze, TianzeCompbio/RQVI_GP_figures). It is produced by script/FigureS5.R, script/FigureS5_rematch.py, and script/FigureS5_plot.py, in the collaborator’s plotting style: two white-to-blue heatmaps sharing row and column order, per-factor relative loading on [0, 1], broad lineage annotations above the columns, and a single shared colorbar.

Cells are our L_pm_filtered cells (those that passed the iterative total-loading filter), restricted to non-thymocytes (annotation_level1 != "thymocyte", all conditions), intersected with the cells carrying RQVI loadings. EBMF and RQVI loadings are averaged within annotation_level2 on those common cells, and each RQVI program is placed on the row of its paired EBMF factor. The collaborator’s EBMF factor F_k was verified equal to our flashier GP_k at the cell level (Pearson r = 1.0 for all 200).

Figure S5

Past versions of S5_ebmf_rqvi_level2_comparison.png

Version	Author	Date
9398c72	Ziang Zhang	2026-07-23

Fig. S5. Cluster-level activity profiles of corresponding EBMF and RQVI factors. EBMF cell loadings (our flashier fit) and matched RQVI program loadings were averaged within 107 fine-grained annotation_level2 clusters, using all non-thymocyte cells present in both the filtered EBMF loading matrix and the RQVI analysis (n = 629,551 cells). The left heatmap shows 200 EBMF factors ordered by hierarchical clustering of their cluster-level profiles. The right heatmap shows, on the same row and column order, the RQVI program assigned to each EBMF factor by a global one-to-one maximum-correlation assignment over all 2,560 candidates (10 seeds x 256 programs), re-derived on this clustering using all common cells. Loadings were independently rescaled to [0, 1] within each factor for display; broad lineage annotations are shown above the columns. The median per-factor Pearson correlation across the displayed clusters was 0.766, and 92.5% of factors had r >= 0.5. Standalone left/right panels and the shared colorbar are saved under figures/generated/Figure S5/S5_subfigures/.

How the figure is made

Three steps, run from the repository root:

Rscript script/FigureS5.R          # cluster-mean matrices + column/palette metadata
python  script/FigureS5_rematch.py # one-to-one EBMF<->RQVI matching on our basis
python  script/FigureS5_plot.py    # Tianze-style heatmaps

The Python steps use an environment with matplotlib, pandas, scipy, h5py, and numpy. The collaborator’s cell-level loading package (data/rqvi_loading/RQVI_EBMF_heatmap_data_v1/) provides the RQVI and EBMF cell-level loadings; data/ is a git-ignored symlink and is not tracked here.

Cluster means

script/FigureS5.R builds the raw annotation_level2 cluster-mean matrices for the EBMF factors and the matched RQVI programs on the common non-thymocyte cells, plus the column order (Figure-1 level1 lineage order, alphabetical within lineage) and the lineage color palette.

# Figure S5 (data step). EBMF vs matched-RQVI level2-cluster comparison.
#
# Design (version B, Tianze-style plotting done in script/FigureS5_plot.py):
#   * Cells: OUR L_pm_filtered cells (passed the iterative total-loading filter),
#     restricted to non-thymocytes (annotation_level1 != "thymocyte"), ALL
#     conditions (no healthy-only restriction), intersected with the cells that
#     have RQVI loadings -> "common cells".
#   * EBMF matrix: our flashier loadings L_pm_filtered (GP1..GP200), averaged
#     within annotation_level2 on the common cells.
#   * RQVI matrix: Tianze's 200 matched RQVI programs (raw cell loadings),
#     averaged within annotation_level2 on the SAME common cells. Each matched
#     program is placed under the column of its paired EBMF factor. The pairing
#     is Tianze's one-to-one match table; F_k == our GP_k (verified at cell level,
#     cell-level Pearson r = 1.0 for all 200).
#   * Columns (level2 clusters) are ordered by the Figure-1 level1 lineage order
#     and alphabetically within each lineage.
#
# This script writes raw cluster-mean matrices + column/palette metadata.
# Row ordering (hierarchical clustering of EBMF), per-factor 0-1 scaling, and
# the heatmaps are produced by script/FigureS5_plot.py in Tianze's visual style.

suppressPackageStartupMessages({
  library(arrow)
  library(data.table)
  library(ZemmourLib)
})

if (!file.exists("code/R/setup_data.R")) {
  stop("Run this script from the immgenT-GP-analysis repository root.")
}
source("code/R/setup_data.R")

outdir <- "figures/generated/Figure S5"
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

PKG <- "data/rqvi_loading/RQVI_EBMF_heatmap_data_v1/data"
parquet_path <- file.path(PKG, "rqvi_matched_200_cell_loadings.parquet")
matches_path <- file.path(PKG, "ebmf_rqvi_multiseed_level2_one_to_one_matches.csv")
level1_order <- c("CD8", "CD4", "Treg", "gdT", "CD8aa", "Tz", "DN", "DP")

## ---- our EBMF loadings + metadata ----
gp <- load_gp_data()
L  <- gp$L_pm_filtered
colnames(L) <- paste0("GP", seq_len(ncol(L)))
meta <- gp$seurat_meta_filtered
stopifnot(identical(rownames(L), rownames(meta)), ncol(L) == 200L)

nonthy <- meta$annotation_level1 != "thymocyte"
if (anyNA(nonthy)) stop("annotation_level1 has missing values.")

## ---- matched RQVI cell loadings (200 programs) ----
pq   <- as.data.frame(arrow::read_parquet(parquet_path))
prog <- setdiff(colnames(pq), "cell_id")
stopifnot(length(prog) == 200L)
matches <- fread(matches_path)
f_for <- matches$ebmf_factor[match(prog, matches$rqvi_candidate)]   # "F<k>" per program
if (anyNA(f_for) || length(unique(f_for)) != 200L) {
  stop("Could not map every RQVI program column to a unique EBMF factor.")
}

## ---- common cells (non-thymocyte, in L_pm, and with RQVI loading) ----
our_nonthy_ids <- rownames(L)[nonthy]
common_ids <- our_nonthy_ids[our_nonthy_ids %in% pq$cell_id]
grp <- factor(as.character(meta$annotation_level2[match(common_ids, rownames(meta))]))
grp <- droplevels(grp)
if (anyNA(grp) || any(as.character(grp) == "")) stop("Missing level2 labels on common cells.")

message(sprintf("common non-thymocyte cells: %d (healthy %d / non-healthy %d); level2 clusters: %d",
                length(common_ids),
                sum(meta$condition_broad[match(common_ids, rownames(meta))] == "healthy"),
                sum(meta$condition_broad[match(common_ids, rownames(meta))] != "healthy"),
                nlevels(grp)))

## ---- EBMF cluster means (columns F1..F200) ----
Lc <- L[common_ids, , drop = FALSE]
colnames(Lc) <- paste0("F", seq_len(ncol(Lc)))
esum <- rowsum(Lc, grp)
ecnt <- as.integer(table(grp)[rownames(esum)])
ebmf_means <- sweep(esum, 1L, ecnt, "/")                 # K x 200

## ---- matched-RQVI cluster means (columns F1..F200, same pairing) ----
pqc <- as.matrix(pq[match(common_ids, pq$cell_id), prog, drop = FALSE])
colnames(pqc) <- f_for
pqc <- pqc[, paste0("F", seq_len(200)), drop = FALSE]     # reorder to F1..F200
rsum <- rowsum(pqc, grp)
rcnt <- as.integer(table(grp)[rownames(rsum)])
rqvi_means <- sweep(rsum, 1L, rcnt, "/")

stopifnot(identical(rownames(ebmf_means), rownames(rqvi_means)),
          identical(ecnt, rcnt),                           # same cells -> same counts
          nrow(ebmf_means) == 107L, ncol(ebmf_means) == 200L)

## ---- level2 column order (Figure-1 level1 order, alphabetical within) ----
l2  <- rownames(ebmf_means)
lmap <- unique(data.frame(level2 = as.character(meta$annotation_level2),
                          level1 = as.character(meta$annotation_level1),
                          stringsAsFactors = FALSE))
if (anyDuplicated(lmap$level2)) stop("A level2 label maps to multiple level1 labels.")
lin <- lmap$level1[match(l2, lmap$level2)]
if (anyNA(lin) || !all(lin %in% level1_order)) stop("Unexpected level1 lineage among clusters.")
ord <- order(match(lin, level1_order), l2)
cluster_order <- data.frame(
  level2_cluster = l2[ord],
  level1         = lin[ord],
  display_column = seq_along(ord) - 1L,
  n_cells        = ecnt[ord],
  stringsAsFactors = FALSE
)

## ---- level1 palette (our canonical colors) as hex ----
l1pal <- ZemmourLib::immgent_colors$level1
hex <- vapply(l1pal, function(cc) {
  v <- grDevices::col2rgb(cc)
  grDevices::rgb(v[1], v[2], v[3], maxColorValue = 255)
}, character(1))
pal_df <- data.frame(level1 = names(hex), color = unname(hex), stringsAsFactors = FALSE)

## ---- write outputs ----
# cell -> level1/level2 for all L_pm_filtered cells; consumed by FigureS5_rematch.py
# to define common cells and clusters (avoids any machine-specific path).
fwrite(data.frame(
  cellID            = rownames(L),
  annotation_level1 = as.character(meta$annotation_level1),
  annotation_level2 = as.character(meta$annotation_level2),
  stringsAsFactors  = FALSE
), file.path(outdir, "S5_cell_metadata.csv.gz"))

fwrite(data.frame(level2_cluster = rownames(ebmf_means), ebmf_means, check.names = FALSE),
       file.path(outdir, "S5_ebmf_raw_means_level2.csv"))
fwrite(data.frame(level2_cluster = rownames(rqvi_means), rqvi_means, check.names = FALSE),
       file.path(outdir, "S5_rqvi_matched_raw_means_level2.csv"))
fwrite(cluster_order, file.path(outdir, "S5_cluster_order.csv"))
fwrite(pal_df,        file.path(outdir, "S5_level1_palette.csv"))
fwrite(data.frame(
  metric = c("common_cells", "healthy_cells", "nonhealthy_cells", "level2_clusters",
             "ebmf_factors", "rqvi_matched_programs"),
  value  = c(length(common_ids),
             sum(meta$condition_broad[match(common_ids, rownames(meta))] == "healthy"),
             sum(meta$condition_broad[match(common_ids, rownames(meta))] != "healthy"),
             nlevels(grp), 200L, 200L)
), file.path(outdir, "S5_build_summary.csv"))

message("Wrote Fig S5 data inputs to ", normalizePath(outdir))

Matching

The matching is our collaborator’s method applied unchanged to the re-aligned data (only the inputs differ: our updated annotation_level2 clustering and our common cell set). Each factor’s mean-loading profile is z-scored across clusters; the signed Pearson correlation between every EBMF factor and all 2,560 seed-specific RQVI candidates (10 seeds x 256 programs) is EBMF_z^T @ RQVI_z / K; candidates with a constant profile are excluded; and a maximum-weight one-to-one assignment (scipy.optimize.linear_sum_assignment) pairs each EBMF factor with a distinct RQVI candidate so that the total signed correlation is maximized. The matching uses all common cells (L_pm_filtered intersect RQVI, 108 level2 clusters including the thymocyte cluster); the figure displays the 107 non-thymocyte clusters. Median per-factor r = 0.766, 92.5% >= 0.5. Full code: script/FigureS5_rematch.py.

Plotting

script/FigureS5_plot.py orders the EBMF rows by hierarchical clustering (average linkage, correlation distance, optimal leaf ordering), rescales every factor to [0, 1] across clusters, and draws the two heatmaps with a shared colorbar; the matched RQVI panel reuses the EBMF row order. Full code: script/FigureS5_plot.py.

Session information

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6.1

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_CA/en_CA/en_CA/C/en_CA/en_CA

time zone: Asia/Shanghai
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] vctrs_0.7.3     cli_3.6.6       knitr_1.50      rlang_1.2.0    
 [5] xfun_0.55       stringi_1.8.7   otel_0.2.0      promises_1.5.0 
 [9] jsonlite_2.0.0  workflowr_1.7.2 glue_1.8.1      rprojroot_2.1.1
[13] git2r_0.36.2    htmltools_0.5.9 httpuv_1.6.16   sass_0.4.10    
[17] rmarkdown_2.30  evaluate_1.0.5  jquerylib_0.1.4 tibble_3.3.0   
[21] fastmap_1.2.0   yaml_2.3.12     lifecycle_1.0.5 whisker_0.4.1  
[25] stringr_1.6.0   compiler_4.5.1  fs_1.6.6        Rcpp_1.1.1-1.1 
[29] pkgconfig_2.0.3 later_1.4.4     digest_0.6.39   R6_2.6.1       
[33] pillar_1.11.1   magrittr_2.0.5  bslib_0.9.0     tools_4.5.1    
[37] cachem_1.1.0