Atera breast cancer Landscape + Clustergram + Enrich¶

This tutorial uses celldega>=0.26.0. It demonstrates the handoff from a precomputed Scanpy analysis plus published cell type annotations to Celldega: Scanpy provides normalized expression and UMAP coordinates, the cell types and colors come from 10x Genomics' own annotation of this dataset, and Celldega calculates set signatures, marker genes for each cell type, an all-gene hierarchy, and linked interactive widgets. The first run downloads the raw 10x file and runs the Scanpy filtering, normalization, PCA, neighbors, and UMAP steps once, saving the result. Later runs load the saved AnnData and go straight to Celldega.

Data source: broadinstitute/Celldega_Supporting_Data. Downloads use a direct Hugging Face /resolve/main/ URL through Python's standard library, so huggingface_hub is not a notebook dependency. Cell types are from 10x Genomics' WTA_Preview_FFPE_Breast_Cancer_cell_groups.csv.

In [1]:

from importlib.metadata import version
from pathlib import Path
import shutil
from urllib.request import Request, urlopen
import warnings

import celldega as dega
import pandas as pd
from pandas.errors import PerformanceWarning
import scanpy as sc

print(f"celldega {version('celldega')}")
print(f"scanpy {version('scanpy')}")

from importlib.metadata import version from pathlib import Path import shutil from urllib.request import Request, urlopen import warnings import celldega as dega import pandas as pd from pandas.errors import PerformanceWarning import scanpy as sc print(f"celldega {version('celldega')}") print(f"scanpy {version('scanpy')}")

celldega 0.26.0
scanpy 1.10.4

Data¶

If the processed AnnData is already on disk it is loaded directly. Otherwise the raw file is downloaded (if needed) and processed with Scanpy.

In [2]:

HF_RESOLVE_BASE = (
    "https://huggingface.co/datasets/"
    "broadinstitute/Celldega_Supporting_Data/resolve/main"
)
RAW_REMOTE_FILENAME = "WTA_Preview_FFPE_Breast_Cancer_outs_AnnData.h5"
CELL_GROUPS_URL = (
    "https://cf.10xgenomics.com/samples/atera/dev/WTA_Preview_FFPE_Breast_Cancer/"
    "WTA_Preview_FFPE_Breast_Cancer_cell_groups.csv"
)

# Store data in the repository's notebooks/data folder when run from a checkout,
# otherwise beside the notebook.
working_dir = Path.cwd().resolve()
project_root = next(
    (
        d
        for d in (working_dir, *working_dir.parents)
        if (d / "pyproject.toml").is_file() and (d / "notebooks").is_dir()
    ),
    None,
)
data_dir = (project_root / "notebooks" if project_root else working_dir) / "data"
data_dir = data_dir / "atera_preview_data_AnnDatas"

raw_path = data_dir / RAW_REMOTE_FILENAME
processed_path = data_dir / "WTA_Preview_FFPE_Breast_Cancer_outs_AnnData_processed.h5ad"
set_collection_path = data_dir / "Atera_cell_type_SetCollection.h5mu"

HF_RESOLVE_BASE = ( "https://huggingface.co/datasets/" "broadinstitute/Celldega_Supporting_Data/resolve/main" ) RAW_REMOTE_FILENAME = "WTA_Preview_FFPE_Breast_Cancer_outs_AnnData.h5" CELL_GROUPS_URL = ( "https://cf.10xgenomics.com/samples/atera/dev/WTA_Preview_FFPE_Breast_Cancer/" "WTA_Preview_FFPE_Breast_Cancer_cell_groups.csv" ) # Store data in the repository's notebooks/data folder when run from a checkout, # otherwise beside the notebook. working_dir = Path.cwd().resolve() project_root = next( ( d for d in (working_dir, *working_dir.parents) if (d / "pyproject.toml").is_file() and (d / "notebooks").is_dir() ), None, ) data_dir = (project_root / "notebooks" if project_root else working_dir) / "data" data_dir = data_dir / "atera_preview_data_AnnDatas" raw_path = data_dir / RAW_REMOTE_FILENAME processed_path = data_dir / "WTA_Preview_FFPE_Breast_Cancer_outs_AnnData_processed.h5ad" set_collection_path = data_dir / "Atera_cell_type_SetCollection.h5mu"

In [3]:

def download_if_missing(remote_filename, destination):
    """Download one public Hub file, with no HF client dependency."""
    if destination.is_file():
        return destination

    destination.parent.mkdir(parents=True, exist_ok=True)
    partial = destination.with_name(f"{destination.name}.part")
    url = f"{HF_RESOLVE_BASE}/{remote_filename}?download=true"
    request = Request(url, headers={"User-Agent": "celldega-example-notebook"})
    print(f"Downloading {url}")
    try:
        with urlopen(request) as response, partial.open("wb") as output:
            shutil.copyfileobj(response, output, length=8 * 1024 * 1024)
        partial.replace(destination)
    except Exception:
        partial.unlink(missing_ok=True)
        raise
    return destination

def download_if_missing(remote_filename, destination): """Download one public Hub file, with no HF client dependency.""" if destination.is_file(): return destination destination.parent.mkdir(parents=True, exist_ok=True) partial = destination.with_name(f"{destination.name}.part") url = f"{HF_RESOLVE_BASE}/{remote_filename}?download=true" request = Request(url, headers={"User-Agent": "celldega-example-notebook"}) print(f"Downloading {url}") try: with urlopen(request) as response, partial.open("wb") as output: shutil.copyfileobj(response, output, length=8 * 1024 * 1024) partial.replace(destination) except Exception: partial.unlink(missing_ok=True) raise return destination

Load or compute the Scanpy preprocessing¶

Preprocessing keeps all genes and saves raw counts in layers["counts"] before log normalization, so SetCollection aggregation and marker ranking both start from untransformed counts. Cell types are added from the 10x annotation in the next section.

In [4]:

if processed_path.is_file():
    print(f"Using existing {processed_path}")
    adata = sc.read_h5ad(processed_path)
else:
    adata = sc.read_10x_h5(download_if_missing(RAW_REMOTE_FILENAME, raw_path))
    sc.pp.filter_cells(adata, min_counts=10)
    sc.pp.filter_genes(adata, min_cells=5)
    adata.layers["counts"] = adata.X.copy()

    sc.pp.normalize_total(adata, target_sum=10_000)
    sc.pp.log1p(adata)
    sc.pp.pca(adata)
    sc.pp.neighbors(adata)
    sc.tl.umap(adata, random_state=0)
    adata.write_h5ad(processed_path)
    print(f"Saved processed AnnData to {processed_path}")

adata

if processed_path.is_file(): print(f"Using existing {processed_path}") adata = sc.read_h5ad(processed_path) else: adata = sc.read_10x_h5(download_if_missing(RAW_REMOTE_FILENAME, raw_path)) sc.pp.filter_cells(adata, min_counts=10) sc.pp.filter_genes(adata, min_cells=5) adata.layers["counts"] = adata.X.copy() sc.pp.normalize_total(adata, target_sum=10_000) sc.pp.log1p(adata) sc.pp.pca(adata) sc.pp.neighbors(adata) sc.tl.umap(adata, random_state=0) adata.write_h5ad(processed_path) print(f"Saved processed AnnData to {processed_path}") adata

Using existing /Users/feni/Documents/celldega/notebooks/data/atera_preview_data_AnnDatas/WTA_Preview_FFPE_Breast_Cancer_outs_AnnData_processed.h5ad

Out[4]:

AnnData object with n_obs × n_vars = 170022 × 18028
    obs: 'n_counts'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells'
    uns: 'log1p', 'neighbors', 'pca', 'umap'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'counts'
    obsp: 'connectivities', 'distances'

10x cell types¶

The 10x cell groups CSV is small (one row per cell: cell_id, group, color), so it is read straight from the URL and propagated onto adata.obs by cell id. The 10x colors are stored as adata.uns["cell_type_colors"] (Scanpy's palette convention), so Scanpy, the Landscape, and the Clustergram all share them.

In [5]:

# The 10x CDN rejects Python's default User-Agent, so send the same one used above.
cell_groups = pd.read_csv(
    CELL_GROUPS_URL,
    index_col="cell_id",
    storage_options={"User-Agent": "celldega-example-notebook"},
)

# A handful of 10x cells were removed by the min_counts filter above; keep only
# AnnData cells that have a 10x label.
labeled = adata.obs_names.isin(cell_groups.index)
print(f"{labeled.sum():,} of {adata.n_obs:,} cells have a 10x cell group")
adata = adata[labeled].copy()

adata.obs["cell_type"] = pd.Categorical(cell_groups.loc[adata.obs_names, "group"])

# One color per group, in category order. A group listed with more than one
# color in the CSV uses its most common one.
palette = cell_groups.groupby("group")["color"].agg(lambda c: c.mode()[0])
adata.uns["cell_type_colors"] = palette[adata.obs["cell_type"].cat.categories].tolist()

adata.obs["cell_type"].value_counts()

# The 10x CDN rejects Python's default User-Agent, so send the same one used above. cell_groups = pd.read_csv( CELL_GROUPS_URL, index_col="cell_id", storage_options={"User-Agent": "celldega-example-notebook"}, ) # A handful of 10x cells were removed by the min_counts filter above; keep only # AnnData cells that have a 10x label. labeled = adata.obs_names.isin(cell_groups.index) print(f"{labeled.sum():,} of {adata.n_obs:,} cells have a 10x cell group") adata = adata[labeled].copy() adata.obs["cell_type"] = pd.Categorical(cell_groups.loc[adata.obs_names, "group"]) # One color per group, in category order. A group listed with more than one # color in the CSV uses its most common one. palette = cell_groups.groupby("group")["color"].agg(lambda c: c.mode()[0]) adata.uns["cell_type_colors"] = palette[adata.obs["cell_type"].cat.categories].tolist() adata.obs["cell_type"].value_counts()

170,022 of 170,022 cells have a 10x cell group

Out[5]:

cell_type
11q13 High Grade DCIS Cells                 64036
CAFs, Low Grade DCIS Associated             24442
Luminal-like Amorphous DCIS Cells           13028
Basal-like Structured DCIS Cells             8760
Endothelial Cells                            8624
Myoepithelial Cells                          8438
Pericytes                                    8087
T Lymphocytes                                8018
CXCL14+ Fibroblasts                          4913
Dendritic Cells                              4008
CAFs, High Grade DCIS Associated             4001
Macrophages                                  3861
Plasma & Mast Cell Mixture                   2512
11q13 High Grade DCIS Cells (Mitotic)        2462
11q13 High Grade DCIS Tumor Cells (G1/S)     2374
Plasma Cells                                  873
Myeloid Cells                                 813
Apocrine Cells                                403
Mast Cells                                    369
Name: count, dtype: int64

In [6]:

sc.pl.umap(adata, color="cell_type")

sc.pl.umap(adata, color="cell_type")

Landscape¶

In [7]:

base_url = (
    "https://raw.githubusercontent.com/cornhundred/"
    "DegaFiles_WTA_Preview_FFPE_Breast_Cancer_outs/main"
)

landscape = dega.viz.Landscape(
    technology="Xenium",
    base_url=base_url,
    height=600,
    adata=adata,
    cell_attr=["cell_type"],
    cluster_attr="cell_type",
)

base_url = ( "https://raw.githubusercontent.com/cornhundred/" "DegaFiles_WTA_Preview_FFPE_Breast_Cancer_outs/main" ) landscape = dega.viz.Landscape( technology="Xenium", base_url=base_url, height=600, adata=adata, cell_attr=["cell_type"], cluster_attr="cell_type", )

SetCollection expression, fraction expressing, and marker ranks¶

Each 10x cell type becomes a set, and Celldega finds the marker genes for each one.

expression.X is the log1p-CPM pseudobulk signature. expression.layers["fraction_expressing"] is the dot-size channel. Marker ranks follow Scanpy's recommendation and are calculated from log-normalized counts (log1p(normalize_total(counts)), computed automatically), and persist in the SetCollection's expression modality.

In [8]:

setc = dega.SetCollection(adata, set_col="cell_type", name="cell_type")

with warnings.catch_warnings():
    warnings.filterwarnings(
        "ignore",
        category=PerformanceWarning,
        module=r"scanpy\.tools\._rank_genes_groups",
    )
    setc.calc_signature(
        adata,
        modality_name="expression",
        layer="counts",
        aggregate="sum",
        normalization="log1p_cpm",
        fraction_expressing_layer="fraction_expressing",
        rank_genes_groups=True,
        rank_genes_groups_kwargs={"method": "t-test"},
    )

setc.write(set_collection_path)
setc.mod["expression"]

setc = dega.SetCollection(adata, set_col="cell_type", name="cell_type") with warnings.catch_warnings(): warnings.filterwarnings( "ignore", category=PerformanceWarning, module=r"scanpy\.tools\._rank_genes_groups", ) setc.calc_signature( adata, modality_name="expression", layer="counts", aggregate="sum", normalization="log1p_cpm", fraction_expressing_layer="fraction_expressing", rank_genes_groups=True, rank_genes_groups_kwargs={"method": "t-test"}, ) setc.write(set_collection_path) setc.mod["expression"]

calc_signature('expression'): sum of adata.layers['counts'] (normalization='log1p_cpm')
  layers['fraction_expressing']: fraction of cells with adata.layers['counts'] > 0.0

  uns['rank_genes_groups']: t-test on log1p(normalize_total(adata.layers['counts'])), grouped by 'cell_type'

Out[8]:

AnnData object with n_obs × n_vars = 19 × 18028
    obs: 'cell_type', 'n_cells', 'set_source', 'color'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'gene', 'cell_type_marker', 'entity_type'
    uns: 'feature_type', 'aggregate', 'normalization', 'layer', 'expr_threshold', 'fraction_expressing_layer', 'rank_genes_groups', 'axis_entities'
    layers: 'fraction_expressing'

Clustergram over all genes¶

The full expression matrix is hierarchically clustered. Marker-ranked views are additional slider levels; they do not replace or prefilter the all-gene matrix.

Columns are colored by 10x cell type with the same palette as the Landscape. Rows are colored by a gene-level grouping derived from those cell types: each gene is assigned to the cell type where it scores highest in rank_genes_groups (expression.var["cell_type_marker"]). The SetCollection uses the cell types to define associated gene groups, kept alongside the signature. Double-click a category name to reorder by it.

In [9]:

marker_levels = [1, 2, 3, 5, 10, 25, 50, 75]

mat = dega.clust.Matrix(
    collection=setc,
    color_by="expression",
    size_by_layer="fraction_expressing",
    # Color columns by cell type and rows by each gene's best-scoring cell type,
    # using the same 10x palette as the Landscape.
    col_attr=["cell_type"],
    row_attr=["cell_type_marker"],
)

# Transform only the Matrix copy; the persisted SetCollection remains unchanged.
mat.norm("row", by="zscore")
mat.cluster(view="rank_genes_groups", levels=marker_levels)

cgm = dega.viz.Clustergram(
    matrix=mat,
    width=450,
    height=500,
    manual_col_cat="cell_type",
    manual_row_cat="info",
    rank_dim=3,
)

marker_levels = [1, 2, 3, 5, 10, 25, 50, 75] mat = dega.clust.Matrix( collection=setc, color_by="expression", size_by_layer="fraction_expressing", # Color columns by cell type and rows by each gene's best-scoring cell type, # using the same 10x palette as the Landscape. col_attr=["cell_type"], row_attr=["cell_type_marker"], ) # Transform only the Matrix copy; the persisted SetCollection remains unchanged. mat.norm("row", by="zscore") mat.cluster(view="rank_genes_groups", levels=marker_levels) cgm = dega.viz.Clustergram( matrix=mat, width=450, height=500, manual_col_cat="cell_type", manual_row_cat="info", rank_dim=3, )

Linked Landscape, Clustergram, and Enrich¶

In [10]:

linked_widgets = dega.viz.spatial_clustergram(
    landscape,
    cgm,
    enrich=True,
)

linked_widgets = dega.viz.spatial_clustergram( landscape, cgm, enrich=True, )

In [11]:

linked_widgets

linked_widgets

Out[11]: