Markers & expression

Fourth of five demonstration notebooks – this one covers marker-gene and pseudobulk expression plots: rank_genes_matrix_and_dot, annotation_marker_matrixplot, annotation_marker_stacked_violin, pseudobulk_boxplot, pseudobulk_multigene_boxplot, and ridgeline_plot.

Data: the same PDAC atlas subset used throughout this set of notebooks (Lucarelli et al., https://doi.org/10.64898/2026.03.19.712924). See Overview for dataset overview plots, Hierarchy & enrichment for hierarchy/enrichment plots, or API Reference for the complete parameter reference.

import anndata as ad
import scanpy as sc

import scplotkit as spk

sc.settings.verbosity = 1

ATLAS_PATH = "/data/Daniele/atlases/final_versions_topublish/Human_subsampled.zarr"

Load the atlas

adata = ad.read_zarr(ATLAS_PATH)
adata

AnnData object with n_obs × n_vars = 100000 × 39041
    obs: 'Sample_ID', 'Condition', 'Treatment', 'TreatmentType', 'TreatmentStatus', 'Tissue', 'Sex', 'Dataset', 'Technology', 'Level_1', 'Level_2', 'Level_3', 'Level_4', 'Age', 'Diabetes', 'Is_Core', 'EMT category', 'Dataset_ID', 'Cluster_Names'
    var: 'n_cells', 'ensembl_id', 'start', 'end', 'chromosome', 'gene_name_adata_sc', 'highly_variable_adata_sc', 'means_adata_sc', 'dispersions_adata_sc', 'dispersions_norm_adata_sc', 'highly_variable_nbatches_adata_sc', 'highly_variable_intersection_adata_sc', 'n_cells_by_counts_adata_sc', 'mean_counts_adata_sc', 'log1p_mean_counts_adata_sc', 'pct_dropout_by_counts_adata_sc', 'total_counts_adata_sc', 'log1p_total_counts_adata_sc', 'mito_adata_sc', 'n_cells_by_counts_adata_sn', 'mean_counts_adata_sn', 'log1p_mean_counts_adata_sn', 'pct_dropout_by_counts_adata_sn', 'total_counts_adata_sn', 'log1p_total_counts_adata_sn', 'Manual_Genes', 'mt', 'ribo', 'hb'
    uns: 'UMAP_0.85', 'Condition_colors', 'Is_Core_colors', 'Level_3_colors', 'UMAP_0.75', 'umap', 'Level_4_colors', 'neighbors', '_scvi_uuid', 'Level_2_colors', 'Level_1_colors', 'Level_4_All_colors', '_scvi_manager_uuid', 'Level_4_Final_colors'
    obsm: 'X_pca', 'pct_counts_mito', 'is_outlier_total_counts', 'scanvi_extended_atlas_emb', 'UMAP_0.85', 'scANVI_cross_species', 'log_counts', 'bin_edges', 'scanvi_L4_emb', '_scvi_batch', 'level0_leiden_subcluster', 'leiden_0.2_annotation', 'n_genes_by_counts', 'log1p_total_counts', 'outlier', 'total_counts', 'log1p_n_genes_by_counts', 'batch', 'total_counts_mito', 'n_genes', 'EMT_score_DL', 'leiden', 'UMAP_0.75', 'MALAT1_lognorm', 'UMAP_0.95', 'infercnv_score_malignant', 'Global_Leiden', 'X_umap', 'log1p_total_counts_mito', 'leiden_0.2', 'leiden_subcluster', 'EMT score', 'scANVI_emb_final', 'cnv_score_abs', 'n_counts', 'leiden_0.5', 'mt_frac', 'infercnv_score_malignant_refined', 'empty_droplet', '_scvi_labels'
    layers: 'counts', 'raw', 'log_norm'
    obsp: 'distances', 'connectivities'

Set up the plotter

config = spk.PlotConfig({
    "palettes": {
        "Level_2": {
            "Lymphoid": "#1f78b4",
            "Malignant Cell": "#e31a1c",
            "Stromal Cell": "#33a02c",
            "Myeloid": "#ff7f00",
            "Endothelial Cell": "#6a3d9a",
            "Exocrine Cell": "#b15928",
            "Endocrine Cell": "#a6cee3",
        },
        "TreatmentStatus": {"Untreated": "#4daf4a", "Treated": "#984ea3"},
    }
})
plotter = spk.ScPlotter(config=config, output_dir="figures")

rank_genes_matrix_and_dot

Data-driven markers via sc.tl.rank_genes_groups, shown as both a dotplot and a matrixplot. Two demonstrations: colored by log fold change (the default), and colored by mean expression via use_expression=True.

The config’s rank_genes defaults to a fixed vmin=-5, vmax=5 diverging scale, but the whole point of n_markers is to surface each group’s most discriminative genes – real log fold changes for those routinely run well past 5, so a fixed +-5 range clips them all to the same saturated end of the colormap. Passing vmin=None, vmax=None here overrides the config default and lets the colormap span the actual range of this atlas’s top marker logFCs instead.

plotter.rank_genes_matrix_and_dot(
    adata,
    groupby_column="Level_2",
    n_markers=4,
    layer="log_norm",
    vmin=None,
    vmax=None,
    save_name="rank_genes",
);

Computing marker genes for 'Level_2' using wilcoxon method...
Generating dotplot...
Figure saved to: figures/marker_analysis/rank_genes_dotplot.png
Generating matrixplot...
Figure saved to: figures/marker_analysis/rank_genes_matrixplot.png

use_expression=True colors dots/cells by mean expression (scaled per gene via standard_scale) instead of log fold change.

plotter.rank_genes_matrix_and_dot(
    adata,
    groupby_column="Level_2",
    n_markers=4,
    layer="log_norm",
    use_expression=True,
    save_name="rank_genes_expression",
);

Computing marker genes for 'Level_2' using wilcoxon method...
Generating dotplot...
Figure saved to: figures/marker_analysis/rank_genes_expression_dotplot.png
Generating matrixplot...
Figure saved to: figures/marker_analysis/rank_genes_expression_matrixplot.png

annotation_marker_matrixplot & annotation_marker_stacked_violin

Matrixplot/violin of a hand-picked, PDAC-relevant marker set rather than data-driven ones. Two demonstrations: the default coolwarm colormap as a matrixplot, and the same markers as a stacked violin with cmap="viridis".

known_markers = {
    "Lymphoid": ["CD3D", "CD3E", "MS4A1"],
    "Malignant Cell": ["S100A6", "TOP2A"],
    "Stromal Cell": ["COL1A1", "DCN", "PDGFRB"],
    "Myeloid": ["LYZ", "CD68"],
    "Endothelial Cell": ["PECAM1", "VWF"],
    "Exocrine Cell": ["PRSS1", "CPA1"],
    "Endocrine Cell": ["INS", "GCG"],
}
known_markers = {ct: [g for g in genes if g in adata.var_names] for ct, genes in known_markers.items()}

plotter.annotation_marker_matrixplot(
    adata,
    markers_dict=known_markers,
    groupby_column="Level_2",
    layer="log_norm",
    save_name="known_markers",
);

Figure saved to: figures/marker_analysis/known_markers_matrixplot_markers.png

plotter.annotation_marker_stacked_violin(
    adata,
    markers_dict=known_markers,
    groupby_column="Level_2",
    layer="log_norm",
    cmap="viridis",
    save_name="known_markers_violin",
);

Figure saved to: figures/marker_analysis/known_markers_violin_stacked_violin_markers.png

pseudobulk_boxplot

Per-sample pseudobulk (mean expression per sample) avoids pseudo-replication when comparing conditions across samples with many cells each. Two demonstrations: EPCAM in malignant cells across real treated/untreated patients, and a single-gene, no-group_by overall view restricted to a cell type with point shape mapped to Sex.

plotter.pseudobulk_boxplot(
    adata,
    gene="EPCAM",
    layer="log_norm",
    sample_column="Sample_ID",
    group_by="TreatmentStatus",
    celltype_column="Level_2",
    celltype_of_interest="Malignant Cell",
    save_name="epcam_pseudobulk",
);

Figure saved to: figures/pseudobulk/epcam_pseudobulk_box.png

Without group_by, a single overall box per gene – here PDGFRB in stromal cells, with point shape mapped to Sex.

plotter.pseudobulk_boxplot(
    adata,
    gene="PDGFRB",
    layer="log_norm",
    sample_column="Sample_ID",
    celltype_column="Level_2",
    celltype_of_interest="Stromal Cell",
    marker_column="Sex",
    save_name="pdgfrb_pseudobulk_overall",
);

Figure saved to: figures/pseudobulk/pdgfrb_pseudobulk_overall_box.png

pseudobulk_multigene_boxplot

The same idea extended to multiple genes at once (min-max scaled per gene), comparing CD4 and CD8 T cells (Level_3) across treatment status, with point shape mapped to Sex.

plotter.pseudobulk_multigene_boxplot(
    adata,
    genes=["CD8A", "IL7R", "FOXP3"],
    layer="log_norm",
    color_column="Diabetes",
    sample_column="Sample_ID",
    celltype_column="Level_3",
    marker_column="Sex",
    celltype_of_interest=["CD4+ T Cell", "CD8+ T Cell"],
    save_name="t_cell_pseudobulk",
);

Figure saved to: figures/pseudobulk/t_cell_pseudobulk.png

ridgeline_plot

Distribution of a continuous QC metric across categories, one ridge per group. Two demonstrations: total_counts across Level_2 compartments with one KDE per TreatmentStatus overlaid on each ridge, and the same metric restricted to the lymphoid compartment at Level_3, z-scored with a sharper bw_adjust.

adata.obs["total_counts"] = adata.obsm["total_counts"].ravel()

plotter.ridgeline_plot(
    adata,
    obs_key="total_counts",
    group_by="Level_2",
    condition_column="TreatmentStatus",
    save_name="total_counts_ridgeline",
);

Figure saved to: figures/ridgeline/total_counts_ridgeline.png

Restricted to the lymphoid compartment at Level_3, z-scored (scale=True) with a sharper bw_adjust for less-smoothed ridges.

plotter.ridgeline_plot(
    adata,
    obs_key="total_counts",
    group_by="Level_3",
    celltype_column="Level_2",
    celltype_of_interest="Lymphoid",
    scale=True,
    bw_adjust=0.5,
    save_name="total_counts_ridgeline_lymphoid",
);

Figure saved to: figures/ridgeline/total_counts_ridgeline_lymphoid.png

Next steps

See Hierarchy & enrichment for hierarchy & enrichment plots, or API Reference for the full parameter reference of every function used here.