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. Scanpy calculates expression signatures and marker rankings in plain AnnData objects; Celldega builds the 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 anndata as ad
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')}")

warnings.filterwarnings("ignore", category=PerformanceWarning, module=r"scanpy\.tools\._rank_genes_groups")

from importlib.metadata import version from pathlib import Path import shutil from urllib.request import Request, urlopen import warnings import anndata as ad 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')}") warnings.filterwarnings("ignore", category=PerformanceWarning, module=r"scanpy\.tools\._rank_genes_groups")

celldega 0.26.1
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"
data_dir.mkdir(parents=True, exist_ok=True)

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" data_dir.mkdir(parents=True, exist_ok=True)

In [3]:

if not processed_path.is_file() and not raw_path.is_file():
    request = Request(f"{HF_RESOLVE_BASE}/{RAW_REMOTE_FILENAME}", headers={"User-Agent": "celldega-example-notebook"})
    partial = raw_path.with_suffix(".h5.part")
    with urlopen(request) as response, partial.open("wb") as output:
        shutil.copyfileobj(response, output)
    partial.replace(raw_path)

if not processed_path.is_file() and not raw_path.is_file(): request = Request(f"{HF_RESOLVE_BASE}/{RAW_REMOTE_FILENAME}", headers={"User-Agent": "celldega-example-notebook"}) partial = raw_path.with_suffix(".h5.part") with urlopen(request) as response, partial.open("wb") as output: shutil.copyfileobj(response, output) partial.replace(raw_path)

Load or compute the Scanpy preprocessing¶

Keep all genes and preserve raw counts in layers["counts"] for aggregation and detection fractions. Normalize and log-transform X for UMAP and marker testing. The saved preprocessing is also used by the neighborhood exploration tutorial; signatures and marker rankings are calculated below on each run. 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(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(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

Saved processed AnnData to /Users/nicolasfernandez/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', 'pca', 'neighbors', 'umap'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'
    layers: 'counts'
    obsp: 'distances', 'connectivities'

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.value_counts(["group", "color"]).reset_index(name="n")
           .sort_values("n", ascending=False).drop_duplicates("group").set_index("group")["color"])
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.value_counts(["group", "color"]).reset_index(name="n") .sort_values("n", ascending=False).drop_duplicates("group").set_index("group")["color"]) 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¶

Start with a full-width, 700px-high tissue view. landscape_overview is a separate widget instance from landscape, which is reserved for the linked Landscape–Clustergram–Enrich view below. Linking and resizing that view therefore leaves this overview at its original size.

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

In [8]:

landscape_overview = dega.viz.Landscape(
    name="atera_overview",
    technology="Xenium",
    base_url=base_url,
    height=700,
    width=0,  # Full available width.
    adata=adata,
    cell_attr=["cell_type"],
    cluster_attr="cell_type",
)
landscape_overview.layout.width = "100%"

landscape_overview = dega.viz.Landscape( name="atera_overview", technology="Xenium", base_url=base_url, height=700, width=0, # Full available width. adata=adata, cell_attr=["cell_type"], cluster_attr="cell_type", ) landscape_overview.layout.width = "100%"

In [9]:

landscape_overview

landscape_overview

Out[9]:

AnnData expression, fraction expressing, and marker ranks¶

sc.get.aggregate summarizes raw counts per 10x cell type. Normalize the sums to log1p-CPM for dot colors and divide nonzero-cell counts by group size for dot sizes. Run sc.tl.rank_genes_groups directly on the log-normalized cell-level adata.X. Its standard Scanpy results are stored in signatures.uns["rank_genes_groups"]; the signatures AnnData contains all genes, their best-scoring cell type, and the same category colors as the Landscape.

In [10]:

signatures = sc.get.aggregate(adata, by="cell_type", func=["sum", "count_nonzero"], layer="counts")
signatures.obs["n_cells"] = signatures.obs["cell_type"].map(adata.obs["cell_type"].value_counts()).astype(int)
signatures.X = signatures.layers["sum"].astype("float32")
signatures.layers["fraction_expressing"] = (
    signatures.layers["count_nonzero"] / signatures.obs["n_cells"].to_numpy()[:, None]
)
sc.pp.normalize_total(signatures, target_sum=1_000_000)
sc.pp.log1p(signatures)
palette = pd.Series(adata.uns["cell_type_colors"], index=adata.obs["cell_type"].cat.categories)
signatures.obs["color"] = signatures.obs["cell_type"].map(palette)
signatures.uns["cell_type_colors"] = palette.loc[signatures.obs["cell_type"]].tolist()

# Test cell-level log-normalized expression, retaining ranks for every gene.
sc.tl.rank_genes_groups(adata, groupby="cell_type", method="t-test",
                       use_raw=False, n_genes=adata.n_vars)
signatures.uns["rank_genes_groups"] = adata.uns["rank_genes_groups"]
markers = sc.get.rank_genes_groups_df(adata, group=None)
# Color each gene by the cluster where its marker score is highest.
best_group = markers.sort_values("scores", ascending=False).drop_duplicates("names").set_index("names")["group"]
signatures.var["cell_type_marker"] = pd.Categorical(
    best_group.reindex(signatures.var_names), categories=adata.obs["cell_type"].cat.categories,
)
signatures

signatures = sc.get.aggregate(adata, by="cell_type", func=["sum", "count_nonzero"], layer="counts") signatures.obs["n_cells"] = signatures.obs["cell_type"].map(adata.obs["cell_type"].value_counts()).astype(int) signatures.X = signatures.layers["sum"].astype("float32") signatures.layers["fraction_expressing"] = ( signatures.layers["count_nonzero"] / signatures.obs["n_cells"].to_numpy()[:, None] ) sc.pp.normalize_total(signatures, target_sum=1_000_000) sc.pp.log1p(signatures) palette = pd.Series(adata.uns["cell_type_colors"], index=adata.obs["cell_type"].cat.categories) signatures.obs["color"] = signatures.obs["cell_type"].map(palette) signatures.uns["cell_type_colors"] = palette.loc[signatures.obs["cell_type"]].tolist() # Test cell-level log-normalized expression, retaining ranks for every gene. sc.tl.rank_genes_groups(adata, groupby="cell_type", method="t-test", use_raw=False, n_genes=adata.n_vars) signatures.uns["rank_genes_groups"] = adata.uns["rank_genes_groups"] markers = sc.get.rank_genes_groups_df(adata, group=None) # Color each gene by the cluster where its marker score is highest. best_group = markers.sort_values("scores", ascending=False).drop_duplicates("names").set_index("names")["group"] signatures.var["cell_type_marker"] = pd.Categorical( best_group.reindex(signatures.var_names), categories=adata.obs["cell_type"].cat.categories, ) signatures

Out[10]:

AnnData object with n_obs × n_vars = 19 × 18028
    obs: 'cell_type', 'n_cells', 'color'
    var: 'gene_ids', 'feature_types', 'genome', 'n_cells', 'cell_type_marker'
    uns: 'log1p', 'cell_type_colors', 'rank_genes_groups'
    layers: 'sum', 'count_nonzero', '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 (signatures.var["cell_type_marker"]). Double-click a category name to reorder by it.

In [11]:

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

mat = dega.clust.Matrix(
    signatures, col_attr=["cell_type"], row_attr=["cell_type_marker"],
    col_entity={"entity": "cell", "attr": "cell_type"},
)
mat.set_size_matrix(ad.AnnData(
    signatures.layers["fraction_expressing"], obs=signatures.obs.copy(), var=signatures.var.copy(),
))

# Z-score only the Matrix copy.
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( signatures, col_attr=["cell_type"], row_attr=["cell_type_marker"], col_entity={"entity": "cell", "attr": "cell_type"}, ) mat.set_size_matrix(ad.AnnData( signatures.layers["fraction_expressing"], obs=signatures.obs.copy(), var=signatures.var.copy(), )) # Z-score only the Matrix copy. 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 [12]:

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

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

In [13]:

linked_widgets

linked_widgets

Out[13]:

Yearbook: visually query individual cells¶

Yearbook displays a grid of tissue portraits centered on individual cells, so you can compare their morphology and local surroundings side by side. Use its query controls to select cells by cluster or cell type, gene expression, or both. A cell-type query samples cells from that group; adding a gene ranks the matching cells by expression, highest first.

This example starts with Myeloid Cells ranked by LST1 expression. Eight portraits are shown per page, with up to 48 cells available through the paging controls. Change the cell type or gene in the Yearbook to explore another population. These queries run in the browser and also work in the saved documentation page.

In [14]:

yearbook = dega.viz.Yearbook(
    name="atera_yearbook",
    base_url=base_url,
    adata=adata,
    cell_attr=["cell_type"],
    cluster_attr="cell_type",
    front_end_query={
        "cluster": {"attr": "cell_type", "value": "Myeloid Cells"},
        "gene": "LST1",
        "max_cells": 48,
    },
    rows=2,
    cols=4,
    portrait_size_um=50,
    height=700,
    width=0,
)
yearbook.layout.width = "100%"

yearbook = dega.viz.Yearbook( name="atera_yearbook", base_url=base_url, adata=adata, cell_attr=["cell_type"], cluster_attr="cell_type", front_end_query={ "cluster": {"attr": "cell_type", "value": "Myeloid Cells"}, "gene": "LST1", "max_cells": 48, }, rows=2, cols=4, portrait_size_um=50, height=700, width=0, ) yearbook.layout.width = "100%"

In [15]:

yearbook

yearbook

Out[15]: