Broad-Allen Developing Mouse Brain E15.5: Viz¶

This tutorial uses celldega>=0.26.1. It visualizes the 3D alignment from Broad-Allen_Dev-Mouse-Brain_E15.5_Align:

1. Visualize in 2D: a Landscape of one slice at full resolution.
2. Explore the clusters in 3D: a wide Composition plot of the Leiden clusters in every slice introduces the clusters along the front-to-back axis. A CellCloud of all 3.8 million aligned cells, linked to a Clustergram of the 480-gene signatures of the clusters, then lets you explore them.

Using the linked widgets. Drag to orbit around the brain and scroll to zoom. Click a column of the Clustergram (a cluster) to highlight that cluster in 3D, or a row (a gene) to color the cells by its expression.

Watch on YouTube

The experiment¶

This dataset is being generated in collaboration with the Allen Institute. It is whole-brain spatial transcriptomics of a mouse embryo at E15.5. The brain was cut into 44 coronal sections, ordered from the front of the brain (slice 1, anterior) to the back (slice 44, posterior), with a section taken every 100 µm, and profiled with the 10x Xenium platform and a 480-gene panel. Cells were segmented from the images, which gives each cell a position (in µm) and a gene-count vector. After quality control the dataset holds about 3.8 million cells.

All cells from all slices were clustered together with Leiden using rapids-singlecell on Manifold.ai, so a cluster means the same thing in every slice. The clustered AnnData (raw counts, cell positions, slice numbers and Leiden labels) is hosted on Hugging Face: Allen-Broad_Dev-Mouse-Brain/E15.5/xenium_E15.5_clustered.h5ad (about 1.9 GB).

The DegaFiles (the 3D CellCloud and the 2D Landscape of slice 35) are hosted on GitHub at broadinstitute/celldega_xenium_E15.5_cell_cloud.

In [1]:

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

import anndata as ad
import numpy as np
import pandas as pd
from pandas.errors import PerformanceWarning

import celldega as dega


# scanpy's rank_genes_groups builds its result table column by column
warnings.filterwarnings("ignore", category=PerformanceWarning)

dega.__version__

from pathlib import Path import shutil from urllib.request import Request, urlopen import warnings import anndata as ad import numpy as np import pandas as pd from pandas.errors import PerformanceWarning import celldega as dega # scanpy's rank_genes_groups builds its result table column by column warnings.filterwarnings("ignore", category=PerformanceWarning) dega.__version__

Out[1]:

'0.26.1'

Data¶

The widgets read their DegaFiles straight from GitHub. The clustered AnnData, used for the composition and the gene signatures, is downloaded from Hugging Face unless it is already on disk.

In [2]:

GITHUB_BASE = (
    "https://raw.githubusercontent.com/broadinstitute/celldega_xenium_E15.5_cell_cloud/main"
)
CELL_CLOUD_URL = f"{GITHUB_BASE}/xenium_E15.5_cell_cloud"  # all aligned cells, 3D
LANDSCAPE_URL = f"{GITHUB_BASE}/output-XETG00323__0092186__35__20260814__172429"  # slice 35, 2D
CLUSTER_KEY = "leiden_1.0"

# Initial camera of the 3D views. rotation_x tilts from side-on (0) to top-down (90);
# rotation_orbit swings the camera around the vertical axis. Anterior is on the left.
ROTATION_X = 20
ROTATION_ORBIT = -25

GITHUB_BASE = ( "https://raw.githubusercontent.com/broadinstitute/celldega_xenium_E15.5_cell_cloud/main" ) CELL_CLOUD_URL = f"{GITHUB_BASE}/xenium_E15.5_cell_cloud" # all aligned cells, 3D LANDSCAPE_URL = f"{GITHUB_BASE}/output-XETG00323__0092186__35__20260814__172429" # slice 35, 2D CLUSTER_KEY = "leiden_1.0" # Initial camera of the 3D views. rotation_x tilts from side-on (0) to top-down (90); # rotation_orbit swings the camera around the vertical axis. Anterior is on the left. ROTATION_X = 20 ROTATION_ORBIT = -25

To view DegaFiles you wrote yourself (Align notebook, section 8), serve them from a local server and point CELL_CLOUD_URL there instead:

CELL_CLOUD_URL = f"http://localhost:{dega.viz.get_local_server()}/DegaFiles/xenium_E15.5_cell_cloud"

In [3]:

HF_RESOLVE_BASE = (
    "https://huggingface.co/datasets/broadinstitute/Celldega_Supporting_Data/resolve/main"
)
H5AD_REMOTE = "Allen-Broad_Dev-Mouse-Brain/E15.5/xenium_E15.5_clustered.h5ad"

DATA_DIR = Path("data/Allen_Dev_Mouse_Brain/E15.5")
H5AD_PATH = DATA_DIR / "xenium_E15.5_clustered.h5ad"


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


download_if_missing(H5AD_REMOTE, H5AD_PATH)

HF_RESOLVE_BASE = ( "https://huggingface.co/datasets/broadinstitute/Celldega_Supporting_Data/resolve/main" ) H5AD_REMOTE = "Allen-Broad_Dev-Mouse-Brain/E15.5/xenium_E15.5_clustered.h5ad" DATA_DIR = Path("data/Allen_Dev_Mouse_Brain/E15.5") H5AD_PATH = DATA_DIR / "xenium_E15.5_clustered.h5ad" 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 download_if_missing(H5AD_REMOTE, H5AD_PATH)

Out[3]:

PosixPath('data/Allen_Dev_Mouse_Brain/E15.5/xenium_E15.5_clustered.h5ad')

In [4]:

adata = ad.read_h5ad(H5AD_PATH)  # raw counts in X (sparse)

# The same cells as the 3D stack: slice 44 and the extra piece of tissue on slice 13 are left out
in_stack = ((adata.obs["slice"] != 44) & ~adata.obs["extra_tissue_piece"]).to_numpy()
cells = adata.obs[in_stack]
cluster_colors = dict(
    zip(adata.obs[CLUSTER_KEY].cat.categories, adata.uns[f"{CLUSTER_KEY}_colors"], strict=True)
)
print(f"{len(cells):,} cells, {cells['slice'].nunique()} slices, {adata.n_vars} genes")

adata = ad.read_h5ad(H5AD_PATH) # raw counts in X (sparse) # The same cells as the 3D stack: slice 44 and the extra piece of tissue on slice 13 are left out in_stack = ((adata.obs["slice"] != 44) & ~adata.obs["extra_tissue_piece"]).to_numpy() cells = adata.obs[in_stack] cluster_colors = dict( zip(adata.obs[CLUSTER_KEY].cat.categories, adata.uns[f"{CLUSTER_KEY}_colors"], strict=True) ) print(f"{len(cells):,} cells, {cells['slice'].nunique()} slices, {adata.n_vars} genes")

3,787,366 cells, 43 slices, 480 genes

1. Visualize in 2D: a Landscape of one slice¶

A Landscape shows one slice at full resolution: the DAPI and cell-boundary images, the cell outlines, and the transcripts of every gene. Cells are colored by their Leiden cluster, and the UMAP of the slice's cells is available from the control panel. Here we show slice 35, in the posterior part of the brain.

The cells in the AnnData are named <slice>_<cell id>. cell_name_prefix=True drops the slice prefix so the names match the cell ids of the slice's Xenium output.

In [5]:

landscape = dega.viz.Landscape(
    technology="Xenium",
    base_url=LANDSCAPE_URL,
    adata=adata[adata.obs["slice"] == 35],
    cell_attr=[CLUSTER_KEY],
    cluster_attr=CLUSTER_KEY,
    cell_name_prefix=True,
    rotate=180,
    height=650,
)
landscape

landscape = dega.viz.Landscape( technology="Xenium", base_url=LANDSCAPE_URL, adata=adata[adata.obs["slice"] == 35], cell_attr=[CLUSTER_KEY], cluster_attr=CLUSTER_KEY, cell_name_prefix=True, rotate=180, height=650, ) landscape

Out[5]:

2. Explore the clusters in 3D¶

Cluster composition across slices¶

The composition plot shows, for every slice (one bar each, anterior to posterior), the share of cells in each Leiden cluster. It gives an overview of all clusters at once and shows which ones belong to which part of the brain along the front-to-back axis.

In [6]:

counts = pd.crosstab(cells["slice"], cells[CLUSTER_KEY])
counts.index = counts.index.astype(str)
counts.columns = counts.columns.astype(str)
props = counts.div(counts.sum(axis=1), axis=0)

composition = dega.viz.Composition(
    props,
    colors=cluster_colors,
    col_weights=counts.sum(axis=1).to_dict(),  # cells per slice, for the COUNTS view
    cluster=False,  # keep the bars in slice order
    width=1200,
    height=260,
)
composition

counts = pd.crosstab(cells["slice"], cells[CLUSTER_KEY]) counts.index = counts.index.astype(str) counts.columns = counts.columns.astype(str) props = counts.div(counts.sum(axis=1), axis=0) composition = dega.viz.Composition( props, colors=cluster_colors, col_weights=counts.sum(axis=1).to_dict(), # cells per slice, for the COUNTS view cluster=False, # keep the bars in slice order width=1200, height=260, ) composition

Out[6]:

Gene-expression signatures of the clusters¶

A SetCollection treats each Leiden cluster as a set of cells and computes its gene-expression signature: log1p counts per million summed over the cluster, the fraction of the cluster's cells expressing each gene, and the top marker genes of each cluster (rank_genes_groups).

Marker ranking works on log-normalized expression of every cell, which for 3.8 million cells needs far more memory than a laptop has. The signatures are therefore computed from a random sample of up to 20,000 cells per cluster (all cells of smaller clusters). For this dataset the sampled signatures correlate with those of all cells at r > 0.999 for every cluster.

In [7]:

%%time
MAX_CELLS_PER_CLUSTER = 20_000

rng = np.random.default_rng(0)
rows = pd.Series(np.flatnonzero(in_stack), index=cells.index)
picked = np.sort(
    np.concatenate(
        [
            rng.choice(r.to_numpy(), min(len(r), MAX_CELLS_PER_CLUSTER), replace=False)
            for _, r in rows.groupby(cells[CLUSTER_KEY], observed=True)
        ]
    )
)
sample = adata[picked].copy()
sample.layers["counts"] = sample.X  # raw counts
print(f"{sample.n_obs:,} cells sampled")

setc = dega.SetCollection(sample, set_col=CLUSTER_KEY, name="xenium_cluster")
setc.calc_signature(
    sample,
    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.mod["expression"]

%%time MAX_CELLS_PER_CLUSTER = 20_000 rng = np.random.default_rng(0) rows = pd.Series(np.flatnonzero(in_stack), index=cells.index) picked = np.sort( np.concatenate( [ rng.choice(r.to_numpy(), min(len(r), MAX_CELLS_PER_CLUSTER), replace=False) for _, r in rows.groupby(cells[CLUSTER_KEY], observed=True) ] ) ) sample = adata[picked].copy() sample.layers["counts"] = sample.X # raw counts print(f"{sample.n_obs:,} cells sampled") setc = dega.SetCollection(sample, set_col=CLUSTER_KEY, name="xenium_cluster") setc.calc_signature( sample, 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.mod["expression"]

756,381 cells sampled

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 'leiden_1.0'

CPU times: user 23.2 s, sys: 9.07 s, total: 32.2 s
Wall time: 18.5 s

Out[7]:

AnnData object with n_obs × n_vars = 38 × 480
    obs: 'leiden_1.0', 'n_cells', 'set_source', 'color'
    var: 'gene', 'leiden_1.0_marker', 'entity_type'
    uns: 'feature_type', 'aggregate', 'normalization', 'layer', 'expr_threshold', 'fraction_expressing_layer', 'rank_genes_groups', 'axis_entities'
    layers: 'fraction_expressing'

Each column of the Clustergram is a Leiden cluster and each row is a gene. Color is the row-z-scored signature and the size of a square is the fraction of the cluster's cells that express the gene. Rows and columns are hierarchically clustered. The DIM slider keeps only the top 1, 2, 3, 5, 10 or 25 marker genes per cluster, so you can go from a handful of markers to the whole panel.

In [8]:

mat = dega.clust.Matrix(
    collection=setc,
    color_by="expression",
    size_by_layer="fraction_expressing",
    col_attr=[CLUSTER_KEY],
    global_colors=cluster_colors,
)
mat.norm("row", by="zscore")
mat.cluster(view="rank_genes_groups", levels=[1, 2, 3, 5, 10, 25])
clustergram = dega.viz.Clustergram(matrix=mat, width=480, height=500, rank_dim=3)

mat = dega.clust.Matrix( collection=setc, color_by="expression", size_by_layer="fraction_expressing", col_attr=[CLUSTER_KEY], global_colors=cluster_colors, ) mat.norm("row", by="zscore") mat.cluster(view="rank_genes_groups", levels=[1, 2, 3, 5, 10, 25]) clustergram = dega.viz.Clustergram(matrix=mat, width=480, height=500, rank_dim=3)

Linked 3D view¶

A CellCloud of every aligned cell, linked to the Clustergram. Click a column of the Clustergram to highlight that cluster in 3D, or a row to color the cells by that gene's expression.

The CellCloud holds both alignments of the same cells. Pass alignment="centroids" to see the result of step 1 alone.

In [9]:

cell_cloud = dega.viz.CellCloud(
    base_url=CELL_CLOUD_URL,
    alignment="neighborhood",
    cluster_attr=CLUSTER_KEY,
    rotation_x=ROTATION_X,
    rotation_orbit=ROTATION_ORBIT,
    height=700,
    cell_size=10,
)
dega.viz.spatial_clustergram(cell_cloud, clustergram, width="780px", height="700px")

cell_cloud = dega.viz.CellCloud( base_url=CELL_CLOUD_URL, alignment="neighborhood", cluster_attr=CLUSTER_KEY, rotation_x=ROTATION_X, rotation_orbit=ROTATION_ORBIT, height=700, cell_size=10, ) dega.viz.spatial_clustergram(cell_cloud, clustergram, width="780px", height="700px")

Out[9]: