Generating pseudobulk from single-cell and single-nucleus RNA-seq • powerup

Single-cell RNA-seq or single-nucleus RNA-seq can be used with POWERUP after cells are aggregated into pseudobulk profiles.

This vignette describes a workflow for generating pseudobulk expression profiles:

raw UMI counts -> QC -> clustering -> sum counts by group -> CPM -> log2(CPM + 1)

The result is a sample-by-gene expression matrix that can be supplied to prepare_powerup_data() as the user_matrix input.

This example begins from 10x HDF5 count matrices (.h5 files). Please note that the parameters used here for processing the data (e.g. QC, clustering) are for demonstration only.

Preserve the full raw-count matrix

First, we want to keep the original UMI counts separate from the processed Scanpy object used for clustering.

So in this example, we use two objects:

adata_counts contains the complete original UMI count matrix and full gene set.
adata is a working copy used for QC, normalization, dimensionality reduction, and clustering.

The final pseudobulk profiles are generated from adata_counts.X, restricted to cells that survive QC in adata and grouped using labels stored in adata.obs.

You may also use adata.raw as the source for pseudobulk counts if you have explicitly verified that it still contains the untouched full raw-count matrix.

Load raw count matrices

The example below combines three samples, each with its own count matrix, generated from two batches. Cell barcodes are prefixed with the sample identifier so they remain unique after concatenation.

import anndata as ad
import scanpy as sc
from pathlib import Path

work_dir = Path("single_cell_powerup")
work_dir.mkdir(parents=True, exist_ok=True)

samples = {
    "sample_A": {"file": "data/sample_A_filtered.h5", "batch": "batch_1"},
    "sample_B": {"file": "data/sample_B_filtered.h5", "batch": "batch_1"},
    "sample_C": {"file": "data/sample_C_filtered.h5", "batch": "batch_2"}
}

adatas = []

for sample_id, info in samples.items():
    sample_adata = sc.read_10x_h5(info["file"])
    sample_adata.var_names_make_unique()

    sample_adata.obs["barcode"] = sample_adata.obs_names
    sample_adata.obs_names = [f"{sample_id}_{x}" for x in sample_adata.obs_names]
    sample_adata.obs["sample"] = sample_id
    sample_adata.obs["batch"] = info["batch"]

    adatas.append(sample_adata)

adata_counts = ad.concat(adatas, join="outer", fill_value=0)

At this point, adata_counts.X should contain untouched UMI counts. Keep this object unchanged.

Perform expression QC

Create a separate working copy and perform expression-based QC in that object.

adata = adata_counts.copy()

sc.pp.filter_genes(adata, min_cells=100)
sc.pp.filter_cells(adata, min_genes=300)

adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

adata = adata[
    (adata.obs["n_genes_by_counts"] < 6000)
    & (adata.obs["pct_counts_mt"] < 20)
].copy()

These thresholds are examples, not POWERUP requirements. The important result is that adata.obs_names now defines the cells that will contribute to the final pseudobulk profiles. Gene filtering here applies only to the processed analysis object; the original genes remain available in adata_counts for pseudobulk aggregation.

Normalize and cluster the cells

Highly variable genes can be used to identify transcriptional groups without restricting the gene space used for the final POWERUP matrix.

sc.pp.highly_variable_genes(adata, n_top_genes=5000, flavor="seurat_v3")

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

adata = adata[:, adata.var["highly_variable"]].copy()

sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata)

If the dataset contains technical batches, batch correction can be applied to the PCA representation. For example, Harmony can be used as follows:

import scanpy.external as sce

sce.pp.harmony_integrate(adata, "batch", max_iter_harmony=20)
adata.obsm["X_pca"] = adata.obsm["X_pca_harmony"]

Construct the neighbor graph and identify Leiden clusters. Here, the clusters define the groups that will become pseudobulk POWERUP samples.

n_neighbors = round(adata.n_obs ** 0.5)
sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=30)

sc.tl.leiden(adata, flavor="igraph", n_iterations=2, resolution=2, key_added="leiden")
sc.tl.umap(adata, min_dist=0.25, spread=5, random_state=123)

sc.pl.umap(adata, color=["sample", "leiden"])

Leiden clusters are only one possible grouping. A different biologically meaningful label can be used for pseudobulk aggregation as long as it is stored for the QC-retained cells in adata.obs.

In a real application, the resolution of clustering will likely be smaller (i.e. fewer clusters). For POWERUP predictions, my preference is to over-cluster slightly so that variability in predictions could be evaluated. So you may have leiden for the real clustering resolution, and leiden_2 for a second, higher resolution.

Save separate analysis and raw-count checkpoints

You can always save the processed and raw-count objects separately to make the pseudobulk step reproducible without relying on objects remaining in memory.

analysis_file = work_dir / "single_cell_annotated.h5ad"
counts_file = work_dir / "single_cell_raw_counts.h5ad"

adata.write_h5ad(analysis_file, compression="gzip")
adata_counts.write_h5ad(counts_file, compression="gzip")

The processed object supplies the retained cells and their grouping labels. The raw-count object supplies the complete original count matrix.

Generate the POWERUP pseudobulk matrix

When ready to generate the pseudobulk matrix, all you need to do is select the QC-retained cells from the raw-count object, and sum their original UMI counts within each group. In this example, we are using the leiden clusters as the grouping variable.

import numpy as np
import pandas as pd
import scanpy as sc
from pathlib import Path
from scipy import sparse

work_dir = Path("single_cell_powerup")

analysis_file = work_dir / "single_cell_annotated.h5ad"
counts_file = work_dir / "single_cell_raw_counts.h5ad"

adata = sc.read_h5ad(analysis_file)
adata_counts = sc.read_h5ad(counts_file)

cells = adata.obs_names
counts = adata_counts[cells].X
clusters = adata.obs["leiden"].astype(str)

cluster_levels = sorted(clusters.unique(), key=int)
cluster_map = {cluster: i for i, cluster in enumerate(cluster_levels)}
cluster_codes = clusters.map(cluster_map).to_numpy()

group_matrix = sparse.csr_matrix(
    (np.ones(len(cells)), (cluster_codes, np.arange(len(cells)))),
    shape=(len(cluster_levels), len(cells))
)

pseudobulk_counts = group_matrix @ counts

pseudobulk = pd.DataFrame(
    pseudobulk_counts.toarray(),
    index=[f"cluster_{int(cluster):02d}" for cluster in cluster_levels],
    columns=adata_counts.var_names
)

library_size = pseudobulk.sum(axis=1)
bulk_matrix = np.log2(pseudobulk.div(library_size, axis=0) * 1e6 + 1)
bulk_matrix.index.name = "cell_line"

output_file = work_dir / "powerup_pseudobulk.tsv"
bulk_matrix.to_csv(output_file, sep="\t")

print("Pseudobulk samples:", bulk_matrix.shape[0])
print("Genes:", bulk_matrix.shape[1])
print("Output:", output_file)
bulk_matrix.head()

Each row of bulk_matrix is now one pseudobulk sample (e.g. a cluster) and each column is a gene (using the complete raw-count feature space not just the highly variable genes).

Note on expression units

Unlike bulk RNAseq, UMI-based sn/scRNA-seq data contain molecule counts.

For each pseudobulk group, raw UMI counts are first summed by gene and then converted to counts per million:

$\mathrm{CPM}_g = \frac{\mathrm{summed\ UMI\ counts}_g}{\mathrm{total\ UMI\ counts}} \times 10^6$

The final POWERUP expression value is:

$\log_2(\mathrm{CPM}_g + 1)$

The resulting matrix should therefore be described as log2(CPM + 1). We accept these CPMs as approximations to TPMs.

Use the pseudobulk matrix in POWERUP

We can now read the generated TSV in R and supply it to the normal POWERUP data-preparation step. Here, expression and response are the reference datasets and desired_targets contains the responses (perturbations) we wish to model.

user_matrix <- read.delim("single_cell_powerup/powerup_pseudobulk.tsv", check.names = FALSE)

prepared <- prepare_powerup_data(
  gene_expression = expression,
  response = response,
  targets = desired_targets,
  user_matrix = user_matrix
)

POWERUP uses the features shared between the reference expression data and the pseudobulk user matrix. Feature names are matched case-insensitively, and terminal Entrez annotations such as TSPAN8 (7105) are ignored by default unless duplicate base feature names require Entrez-based disambiguation. Prefixes remain part of the feature name and are not stripped.

Note on feature selection for training models

For sc/snRNAseq work, you may be more interested in differential response between clusters (e.g. tumor vs stroma). In that case, you can set feature_selection_source to "user_matrix" in prepare_powerup_data() to select the variable features from the pseudobulk matrix instead of the reference. This will allow POWERUP to focus on the genes that vary between your clusters of interest. The performance of some predictive models may drop because their important features are either not expressed or not variable among the sc/snRNAseq clusters, but the expectation is that such models would have predicted similar response between clusters in that scenario.

Summary

For single-cell or single-nucleus RNA-seq, the recommended preparation is:

raw UMI counts -> preserve full raw-count object -> QC and cluster a separate working object -> select QC-retained cells from raw counts -> sum counts by group -> CPM -> log2(CPM + 1) -> POWERUP user matrix

The key separation between the raw-count and processed objects allows Scanpy analysis to use normalization and highly variable gene selection without changing or restricting the counts used to construct the final POWERUP pseudobulk profiles.

Next step

Once the pseudobulk matrix has been generated, continue with Preparing data for POWERUP to align the user matrix with the reference data and prepare it for model training and prediction.