Human and Mammalian Brain Atlas (HMBA) - Marmoset Subcortex Xenium Spatial Transcriptomics: Gene Expression

The Marmoset Subcortical Cell Atlas (MSCA), is a multiomic atlas of the common marmoset (Callithrix jacchus) subcortex that integrates single nucleus transcriptomics, chromatin accessibility, and spatial transcriptomics. Part of the Human and Mammalian Brain Atlas (HMBA) effort within the NIH’s BRAIN Initiative® Cell Atlas Network (BICAN), this dataset comprehensively samples key subcortical structures in the non-human primate, including the hippocampal formation, amygdala, olfactory areas, basal forebrain, thalamus, hypothalamus, midbrain, and major hindbrain structures (excluding cerebellar cortex).

This cellular-resolution Xenium spatial transcriptomic dataset includes 5.6 million cells profiled from 47 coronal hemisections from a single adult male marmoset donor, sampling the rostral-caudal extent of subcortical structures. The spatial cells were mapped to the 1,737 snRNA-seq Clusters defined in the MSCA reference cell type taxonomy, and the resulting spatial localization of the Clusters provided anatomical anchoring and nomenclature to the Group and Class hierarchy levels of the MSCA taxonomy.

For more information on this dataset and the associated 10X Multiome data, please refer to the publication: Dan, Turner, DeBerardine, et al. (2026). Spatial and regulatory organization of cell identity across the primate subcortex. Cell. https://doi.org/10.1016/j.cell.2026.10.007.

---

In this notebook, we build on the our previous investigation into the cell level metadata, combining the metadata gene expressions for select genes that we extract from the released h5ad files.

You need to be connected to the internet to run this notebook and have run through the getting started notebook.

import json
import geojson
from typing import List, Tuple, Optional

import pandas as pd
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

from abc_atlas_access.abc_atlas_cache.abc_project_cache import AbcProjectCache
from abc_atlas_access.abc_atlas_cache.anndata_utils import get_gene_data

We will interact with the data using the AbcProjectCache. This cache object tracks which data has been downloaded and serves the path to the requested data on disk. For metadata, the cache can also directly serve a up a Pandas DataFrame. See the getting_started notebook for more details on using the cache including installing the cache if it has not already been.

Change the download_base variable to where you would like to download the data in your system.

download_base = Path('../../data/abc_atlas')
abc_cache = AbcProjectCache.from_cache_dir(
    cache_dir=download_base,
)

abc_cache.current_manifest

/Users/chris.morrison/src/abc_atlas_access/src/abc_atlas_access/abc_atlas_cache/cloud_cache.py:665: UserWarning: The manifest version recorded as last used for this cache -- releases/20261130/manifest.json-- is not a valid manifest for this dataset. Loading latest version -- releases/20261101/manifest.json -- instead.
  warnings.warn(msg, UserWarning)
/Users/chris.morrison/src/abc_atlas_access/src/abc_atlas_access/abc_atlas_cache/cloud_cache.py:697: OutdatedManifestWarning: You are loading
releases/20261101/manifest.json
which is newer than the most recent manifest file you have previously been working with
releases/20261130/manifest.json
It is possible that some data files have changed between these two data releases, which will force you to re-download those data files (currently downloaded files will not be overwritten). To continue using releases/20261130/manifest.json, run
type.load_manifest('releases/20261130/manifest.json')
  warnings.warn(msg, OutdatedManifestWarning)

'releases/20261130/manifest.json'

Cell metadata

The Marmoset Subcortext spatial dataset has it’s own h5ad files and associated cell metadata tables. These metadata contain the basic identifiers for each cell and forms the basis by which we will join all other cell level and associated data. In the previous tutorial, we merged the cell metadata tables with their associated specimen metadatadata and then with the rotated coordinates we use to plot the cell locations.

Below we load the base metadata needed for extracting gene expressions.

# Load the cell metadata and specimen metadata
# and join them together with the rotated coordinates.

cell_metadata = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-Xenium-Marmoset',
    file_name='cell_metadata',
    dtype={"cell_label": str},
).set_index('cell_label')

specimen = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-Xenium-Marmoset',
    file_name='specimen_metadata',
    dtype={"cell_label": str}
).set_index('brain_section_label')

cell_metadata = cell_metadata.join(specimen, on='brain_section_label', rsuffix='_specimen')

rotated_coordinates = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-Xenium-Marmoset',
    file_name='rotated_coordinates',
    dtype={"cell_label": str}
).set_index('cell_label')
cell_metadata = cell_metadata.join(
    rotated_coordinates[['x_rotated_mm', 'y_rotated_mm']],
    on='cell_label'
)

del specimen
del rotated_coordinates

cell_metadata.head()

	brain_section_label	brain_section_barcode	x_experiment	y_experiment	segmentation_job_id	segmentation_method	cell_area	nucleus_area	n_nuclei_in_cell	total_counts	...	hemisphere	brain_division_label	brain_division_ordinal	donor_label	slab_plane_ordinal	slab_plane_label	slab_plane_order	slab_z_order	x_rotated_mm	y_rotated_mm
cell_label
1382042951_aaondhbe-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4459.1006	4524.2420	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	89.81579	60.735157	1	339	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.651601	-5.898891
1382042951_aaoneabm-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4668.6323	4653.2130	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	105.98172	74.327190	1	762	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.508436	-6.098991
1382042951_aaonekin-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4593.0610	4594.6304	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	129.50813	74.598130	1	406	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.555394	-6.015698
1382042951_aaoofhco-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4529.3900	4348.2266	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	136.37189	72.475784	1	948	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.518584	-5.763877
1382042951_aaopgffc-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4524.6533	4382.6323	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	129.41782	78.391250	1	941	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.536286	-5.793757

5 rows × 40 columns

Gene expression data

Now we can extract expression for specific genes from the h5ad files and pair it with our gene metadata. In general, we can use the function get_gene_data to extract the expression of specific genes for all the cells in a given dataset or a subset of cells.. For more details on using this convenience function, see the Accessing 10x RNA-seq gene expression data tutorial notebook.

For this tutorial we will use precomputed tables of the expression for specific genes to compare expression both spatially and by the taxonomy.

To use the gene expression data, we first load the set of genes for the subcortex data.

genes = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-Xenium-Marmoset',
    file_name='gene',
).set_index('gene_identifier')
genes.head()

gene.csv: 100%|██████████| 22.4k/22.4k [00:00<00:00, 172kMB/s] 

	gene_symbol	description	molecular_type
gene_identifier
NCBIGene:100393956	ABI3BP	ABI family member 3 binding protein	protein_coding
NCBIGene:100388523	ABL2	ABL proto-oncogene 2, non-receptor tyrosine ki...	protein_coding
NCBIGene:100389219	ABLIM1	actin binding LIM protein 1	protein_coding
NCBIGene:100400228	ACAP2	ArfGAP with coiled-coil, ankyrin repeat and PH...	protein_coding
NCBIGene:100392257	ACHE	acetylcholinesterase (Cartwright blood group)	protein_coding

We define a set of genes to load via their symbol. We reuse the smaller overlapping subset of genes we used in the previous 10X gene expression notebook.

gene_names = ['AQP4', 'DRD1', 'DRD2', 'SLC17A6']

In the cell below, we show how to use the get_gene_data function to extract the expression of specific genes from the h5ad metadata. In the commented code, we combine our abc_cache object, cell metadata, gene metadata, and gene symbols we would like select in our get_gene_data function that downloads and extracts gene expression data from the h5ad file. Note that this function is an example and should not be used on small systems to extract a large number of genes for all cells at once. For a small subset of cells, it can be used to pull expression data from the h5ad file for all genes.

The methods below will pull data from the natural-log, batch-normalized data. These are data where the raw counts are normalized by scaling all cells within each section (per brain_section_barcode) to the mean total counts of a reference cell population, in this case the counts in the cell type Group, Astrocyte. This is done because their population count correlates with the total neuron count across sections and they are distributed across diverse anatomical regions using the following qsudo-code:

section_factor = global_reference_cells_mean / section_reference_cells_mean
batch_normalized_counts = section_factor * raw_counts

For this tutorial we will skip using these functions and instead use a pre-computed gene expression to save time. For specifics on how to use the below functions see this tutorial.

"""
gene_data = get_gene_data(
    abc_atlas_cache=abc_cache,
    all_cells=cell_metadata,
    all_genes=genes,
    selected_genes=gene_names,
    data_type='ln'
)"""

"\ngene_data = get_gene_data(\n    abc_atlas_cache=abc_cache,\n    all_cells=cell_metadata,\n    all_genes=genes,\n    selected_genes=gene_names,\n    data_type='ln'\n)"

We load the precomputed gene expressions for each species below instead of running the processing.

gene_data = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-Xenium-Marmoset',
    file_name='example_gene_expression',
).set_index('cell_label')
gene_data.head()

example_gene_expression.csv: 100%|██████████| 291M/291M [00:43<00:00, 6.64MMB/s]   

	AQP4	DRD1	DRD2	SLC17A6
cell_label
1382042951_aaondhbe-1-20250227	0.000000	0.926574	0.000000	0.0
1382042951_aaoneabm-1-20250227	0.000000	0.000000	0.566971	0.0
1382042951_aaonekin-1-20250227	0.566971	0.566971	0.000000	0.0
1382042951_aaoofhco-1-20250227	0.000000	0.566971	0.000000	0.0
1382042951_aaopgffc-1-20250227	0.000000	0.000000	0.000000	0.0

Now that we’ve loaded our expression matrix tables, indexed by cell_label, we can merge them into our cell metadata tables.

cell_metadata_with_genes = cell_metadata.join(gene_data, on='cell_label')

# clean up memory
del gene_data

Geojson polygons

The dataset includes polygons outlining the hemisphere boundary of each brain section, which can be plotted over rotated coordinates to allow us to orient the subcortical cells in their full hemisection contex. These polygons are provided in geojson format, and here we use the associated geojson python package to load the jsons for use in our visualization.

Additional, manually drawn polygons outlining select regions of the Basal Ganglia (STR, GP, GPe, GPi, VeP, STH, and SN), as well as polygons identifying the cells that are also part of the HMBA-BG dataset (basal_ganglia_and_adjacent), are included as part of the HMBA-BG-Xenium-Marmoset dataset and are described and can be accessed in that dataset’s equivalent tutorial notebook.

First, we load the term set for these data. All cells in this HMBA-Subcortex-Xenium-Marmoset dataset receive the term set label “Subcortical”.

value_sets = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-Xenium-Marmoset',
    file_name='value_sets'
).set_index('label')
pred = value_sets['field'] == 'parcellation_term_symbol'
parcellation_terms = value_sets[pred].copy()
parcellation_terms.head()

	table	field	description	order	external_identifier	parent_label	color_hex_triplet
label
Hemisphere	cell_anatomical_annotations	parcellation_term_symbol	Brain Hemisphere	1	hemisphere	NaN	#000000
Subcortical	cell_anatomical_annotations	parcellation_term_symbol	Subcortical regions	2	subcortical	Hemisphere	#008080
Basal ganglia and adjacent	cell_anatomical_annotations	parcellation_term_symbol	Basal ganglia and adjacent	3	basal_ganglia_and_adjacent	Subcortical	#BEBEBE
STR	cell_anatomical_annotations	parcellation_term_symbol	striatum	4	HOMBA:10333	Basal ganglia and adjacent	#cd86cc
GP	cell_anatomical_annotations	parcellation_term_symbol	globus pallidus	5	HOMBA:10342	Basal ganglia and adjacent	#7AE361

Below we define two functions that will be used in our plotting of the geojsons. The first takes in the data and a matplotlib axis object and plots the geojson data on the provided axis, extracting color and name information from the value/term set. The second loads the data based on the appropriate slab and species provided.

def plot_geo_data(geo_data: dict, ax: plt.Axes, terms=parcellation_terms):
    """Plot a geojson polygon over the input axis.

    Parameters
    ----------
    geo_data: dict
        Geojson formatted dictionary containing Polygons to plot.
    ax: matplotlib.pyplot.axis
        Axis object to plot the geojson data over.
    terms: pandas.DataFrame
        Term list containing the colors and orders for the parcellation, brain region terms. 
    """
    if geo_data is None:
        return None
    # If the geo_data is a FeatureCollection, we need to iterate over its
    # features and plot each one.
    if geo_data['type'] == 'FeatureCollection':
        for sub_data in geo_data['features']:
            plot_geo_data(sub_data, ax)
    elif geo_data['type'] == 'Feature':
        geometry_data = geo_data['geometry']
        name = geo_data['properties']['label']
        # color = terms.loc[name, 'color_hex_triplet'] # Get the color for this object from the term list.
        color = 'black' # Default color for the polygons.
        # If the geometry is a GeometryCollection, we need to iterate over its
        # geometries that make up the complete feature.
        if geometry_data['type'] == 'GeometryCollection':
            for idx, geometry in enumerate(geometry_data['geometries']):
                coords = np.array(list(geojson.coords(geometry)))
                # Only label the first geometry in the collection. Others
                # will be plotted without a label to avoid duplicate legend entries.
                if idx == 0:
                    ax.fill(coords[:, 0],
                            coords[:, 1],
                            edgecolor=color,
                            fill=False,
                            linewidth=3,
                            label=name)
                else:
                    ax.fill(coords[:, 0],
                            coords[:, 1],
                            edgecolor=color,
                            fill=False,
                            linewidth=3)
        else:
            coords = np.array(list(geojson.coords(geometry_data)))
            try:
                ax.fill(coords[:, 0],
                        coords[:, 1],
                        edgecolor=color,
                        fill=False,
                        linewidth=3,
                        label=name)
            except ValueError as e:
                import pdb; pdb.set_trace()
    else:
        print("unknown data", geo_data)

def get_polygons(slab_name: str) -> dict:
    """Load the geojson polygons for the specified slab and species.

    Parameters
    ----------
    slab_name: str
        Name of the slab to load.

    Returns
    -------
    dict
        Geojson formatted dictionary containing the polygons for the specified slab and species.
    """

    base_dir = 'HMBA-Subcortex-Xenium-Marmoset'

    file_path = abc_cache.get_file_path(
        directory=base_dir,
        file_name=slab_name,
    )
    if not file_path.exists():
        return None

    with open(file_path) as jfile:
        data = json.load(jfile)
    
    return data

Plotting the Slabs

We define a helper function plot_slabs to visualize the cells for a specified set of brain sections either by colorized metadata or gene expression. This function can be used to plot individual brain sections as well.

def plot_slabs(
    df: pd.DataFrame,
    feature: str,
    slab_list: List[str],
    cmap: Optional[plt.Colormap] = None,
    fig_width: float = 40,
    fig_height: float = 8,
    plot_polygons: bool = False,
    point_size: float = 1.0,
    fix_bounds: bool = False
) -> Tuple[plt.Figure, np.ndarray]:
    """Plot spatial transcriptomics data.

    Parameters
    ----------
    df: pandas.DataFrame
        cell metadata data frame with x,y locations to plot.
    feature: str
        Expression value or color column to plot.
    slab_list: List[str]
        List of slab_plane ids to plot.
    cmap: matplotlib.colors.Colormap
        Colormap for expression features.
    fig_width: int
        Width of the plotted figure.
    fig_height: int
        Height of the plotted figure.
    plot_polygons: bool
        Overplot geojson polygons. Note that these are note available for human
        and this should be set to False.
    point_size: float
        Size of the points to plot.
    fix_bounds: bool
        If True, fix the x and y bounds across all slabs to be the same.
    """
    if isinstance(slab_list, str):
        slab_list = [slab_list]

    fig, ax = plt.subplots(1, len(slab_list))
    fig.set_size_inches(fig_width, fig_height)

    if len(slab_list) == 1:
        ax = [ax]

    if 'x' in df.columns:
        x_col = 'x'
    elif 'x_rotated_mm' in df.columns:
        x_col = 'x_rotated_mm'
    if 'y' in df.columns:
        y_col = 'y'
    elif 'y_rotated_mm' in df.columns:
        y_col = 'y_rotated_mm'

    if fix_bounds:
        x_min = df[x_col].min() * 0.8
        x_max = df[x_col].max() * 1.2
        y_min = df[y_col].min() * 0.8
        y_max = df[y_col].max() * 1.2

    for idx, slab in enumerate(slab_list):

        if 'slab_plane_label' in df.columns:
            filtered = df[df['slab_plane_label'] == slab]
        elif 'brain_section_label' in df.columns:
            filtered = df[df['brain_section_label'] == slab]
        
        xx = filtered[x_col]
        yy = filtered[y_col]
        vv = filtered[feature]
        
        if cmap is not None:
            ax[idx].scatter(xx, yy, s=point_size, c=vv, marker='.', cmap=cmap)
        else:
            
            ax[idx].scatter(xx, yy, s=point_size, c=vv, marker=".")

        if plot_polygons:
            plot_geo_data(get_polygons(slab), ax[idx])
            
        ax[idx].axis('equal')
        if fix_bounds:
            ax[idx].set_ylim(y_min, y_max)
            ax[idx].set_xlim(x_min, x_max)
        ylim = ax[idx].get_ylim()
        ax[idx].set_xticks([])
        ax[idx].set_yticks([])
        
        ax[idx].set_title(slab)
        ax[idx].legend(loc=0)
        
    plt.subplots_adjust(wspace=0.01, hspace=0.01)
    return fig, ax

Below we plot the spatial expression of our genes across the different species. We select a minium expression value of which effectively masks out zero expression cells from the plots, showing a bit more contrast. Feel free to experiment with different vmin and vmax values in investigate the gene expressions more in depth.

gene_names = cell_metadata_with_genes.columns[
    cell_metadata_with_genes.columns.isin(gene_names)
].to_numpy()
slab_list = [
    'CJ23.56.004.CX.42.Plane.05', 'CJ23.56.004.CX.42.Plane.13', 'CJ23.56.004.CX.43.Plane.05',
    'CJ23.56.004.CX.43.Plane.15', 'CJ23.56.004.CX.44.Plane.03'
]
for gene_symbol in gene_names:
    fig, ax = plot_slabs(
        cell_metadata_with_genes,
        gene_symbol,
        slab_list,
        cmap=plt.cm.magma_r,
        plot_polygons=True,
    )
    fig.suptitle(f'Gene Expression: {gene_symbol}')
    plt.show()

Taxonomy mapping

Finally, let’s load the taxonomy information so that we can plot gene expression versus various taxons in a heatmap for each species. Below, we load and merge all the information we need to plot cell colors and taxonomy as done in previous HMBA-Subcortex tutorial.

Note that the taxonomy mapping also contains neurotransmitter and dominant brain region locations.

Below we load the subcortex taxonomy. For more details on this taxonomy see the notebook here.

cell_to_cluster_membership = abc_cache.get_metadata_dataframe(
    'HMBA-Subcortex-Xenium-Marmoset', 'cell_to_cluster_membership'
).set_index('cell_label')

cluster = abc_cache.get_metadata_dataframe(
    'HMBA-Subcortex-taxonomy-Marmoset', 'cluster'
).set_index('label')
cluster_annotation_term = abc_cache.get_metadata_dataframe(
    'HMBA-Subcortex-taxonomy-Marmoset', 'cluster_annotation_term'
).set_index('label')

cluster_annotation_term_set = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-taxonomy-Marmoset',
    file_name='cluster_annotation_term_set'
).set_index('label')

cluster_to_cluster_annotation_membership = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-taxonomy-Marmoset',
    file_name='cluster_to_cluster_annotation_membership'
).set_index('cluster_annotation_term_label')

membership_with_cluster_info = cluster_to_cluster_annotation_membership.join(
    cluster.reset_index().set_index('cluster_alias')[['number_of_cells']],
    on='cluster_alias'
)

membership_with_cluster_info = membership_with_cluster_info.join(cluster_annotation_term, rsuffix='_anno_term').reset_index()

membership_groupby = membership_with_cluster_info.groupby(
    ['cluster_alias', 'cluster_annotation_term_set_name']
)

cluster_details = membership_groupby['cluster_annotation_term_name'].first().unstack()
cluster_order = membership_groupby['term_order'].first().unstack()

cluster_colors = membership_groupby['color_hex_triplet'].first().unstack()
cluster_colors = cluster_colors[cluster_annotation_term_set['name']]

cluster_to_dominant_region = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-taxonomy-Marmoset',
    file_name='cluster_to_dominant_region',
).set_index('cluster_alias')
cluster_to_neurotransmitter = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-taxonomy-Marmoset',
    file_name='cluster_to_neurotransmitter',
).set_index('cluster_alias')

cell_metadata_with_genes = cell_metadata_with_genes.join(cell_to_cluster_membership)
cell_metadata_with_genes = cell_metadata_with_genes.join(cluster_details, on='cluster_alias')
cell_metadata_with_genes = cell_metadata_with_genes.join(cluster_colors, on='cluster_alias', rsuffix='_color')
cell_metadata_with_genes = cell_metadata_with_genes.join(cluster_order, on='cluster_alias', rsuffix='_order')
cell_metadata_with_genes = cell_metadata_with_genes.join(cluster_to_dominant_region, on='cluster_alias', rsuffix='_region')
cell_metadata_with_genes = cell_metadata_with_genes.join(cluster_to_neurotransmitter, on='cluster_alias', rsuffix='_neurotransmitter')

del cell_to_cluster_membership

cell_metadata.head()

	brain_section_label	brain_section_barcode	x_experiment	y_experiment	segmentation_job_id	segmentation_method	cell_area	nucleus_area	n_nuclei_in_cell	total_counts	...	hemisphere	brain_division_label	brain_division_ordinal	donor_label	slab_plane_ordinal	slab_plane_label	slab_plane_order	slab_z_order	x_rotated_mm	y_rotated_mm
cell_label
1382042951_aaondhbe-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4459.1006	4524.2420	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	89.81579	60.735157	1	339	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.651601	-5.898891
1382042951_aaoneabm-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4668.6323	4653.2130	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	105.98172	74.327190	1	762	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.508436	-6.098991
1382042951_aaonekin-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4593.0610	4594.6304	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	129.50813	74.598130	1	406	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.555394	-6.015698
1382042951_aaoofhco-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4529.3900	4348.2266	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	136.37189	72.475784	1	948	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.518584	-5.763877
1382042951_aaopgffc-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	4524.6533	4382.6323	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	129.41782	78.391250	1	941	...	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1	3.536286	-5.793757

5 rows × 40 columns

def extract_value_set(
    cell_metadata_df: pd.DataFrame,
    input_value_set: pd.DataFrame,
    input_value_set_label: str,
    value_set_label: str=None
):
    """Add color and order columns to the cell metadata dataframe based on the input
    value set.

    Columns are added as {input_value_set_label}_color and {input_value_set_label}_order.

    Parameters
    ----------
    cell_metadata_df : pd.DataFrame
        DataFrame containing cell metadata.
    input_value_set : pd.DataFrame
        DataFrame containing the value set information.
    input_value_set_label : str
        The the column name to extract color and order information for. will be added to the cell metadata.
    """
    if value_set_label is None:
        value_set_label = input_value_set_label
    cell_metadata_df[f'{input_value_set_label}_color'] = input_value_set[
        input_value_set['field'] == value_set_label
    ].loc[cell_metadata_df[input_value_set_label]]['color_hex_triplet'].values
    cell_metadata_df[f'{input_value_set_label}_order'] = input_value_set[
        input_value_set['field'] == value_set_label
    ].loc[cell_metadata_df[input_value_set_label]]['order'].values

extract_value_set(cell_metadata_with_genes, value_sets, 'Cluster_neurotransmitter', 'neurotransmitter')
extract_value_set(cell_metadata_with_genes, value_sets, 'Group_neurotransmitter', 'neurotransmitter')
extract_value_set(cell_metadata_with_genes, value_sets, 'Cluster_dominant_region', 'dominant_region')
extract_value_set(cell_metadata_with_genes, value_sets, 'Group_dominant_region', 'dominant_region')

Plotting heatmap expression data.

We define a helper functions aggregate_by_metadata to compute the average expression for a given category. Note that this is a simple wrapper of the Pandas GroupBy functionality. Other summary statistics beyond just the mean are listed here: https://pandas.pydata.org/docs/reference/groupby.html#dataframegroupby-computations-descriptive-stats

We also define a function to plot the resultant averaged data in a heatmap.

def aggregate_by_metadata(df: pd.DataFrame, gnames: List[str], value: str, sort: bool = False) -> pd.DataFrame:
    """Aggregate gene expression data by metadata.

    Parameters
    ----------
    df: pandas.DataFrame
        DataFrame containing gene expression and metadata.
    gnames: List[str]
        List of gene names to aggregate.
    value: str
    sort: bool
        Metadata column to group by.
        If True, sort the output by the first gene in gnames.

    Returns
    -------
    pandas.DataFrame
        DataFrame with aggregated gene expression values.
    """
    grouped = df.groupby(value)[gnames].mean()
    if sort:
        grouped = grouped.sort_values(by=gnames[0], ascending=False)
    return grouped


def plot_heatmap(
        df: pd.DataFrame,
        fig_width: int = 8,
        fig_height: int = 4,
        cmap: mpl.colormaps = plt.cm.magma,
        vmax: Optional[float] = None,
        unit: str = '[ln(batch-normalized + 1)]'
    ) -> Tuple[plt.Figure, plt.Axes]:
    """Plot a heatmap from a DataFrame.

    Parameters
    ----------
    df: pandas.DataFrame
        DataFrame to plot as a heatmap.
    fig_width: int
        Width of the figure.
    fig_height: int
        Height of the figure.
    cmap: matplotlib.colors.Colormap
        Colormap to use for the heatmap.

    Returns
    -------
    fig: matplotlib.pyplot.Figure
        Figure object containing the heatmap.
    ax: matplotlib.pyplot.Axes
        Axes object containing the heatmap.
    """

    arr = df.to_numpy().astype('float')

    vmin = arr.min()
    vmax = arr.max()
    norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    cmap = sm.get_cmap()

    fig, ax = plt.subplots()
    fig.set_size_inches(fig_width, fig_height)

    ax.imshow(arr, cmap=cmap, aspect='auto')
    xlabs = df.columns.values
    ylabs = df.index.values

    ax.set_xticks(range(len(xlabs)))
    ax.set_xticklabels(xlabs)

    ax.set_yticks(range(len(ylabs)))
    ax.set_yticklabels(ylabs)

    cbar = fig.colorbar(sm, ax=ax, orientation='vertical', fraction=0.1, pad=0.01)
    cbar.set_label(unit)

    return fig, ax

Visualize gene expression for selected genes by the HMBA-Subcortex taxonomy.

Below we average the expression of our selected genes across taxons from the HMBA-Subcortex taxonomy. We show the expression for taxons at the class level for each species. For a side by side comparison of gene expression heatmaps for each species, see the HMBA 10X gene expression tutorial.

gene_names = cell_metadata_with_genes.columns[
    cell_metadata_with_genes.columns.isin(gene_names)
].to_numpy()
agg = aggregate_by_metadata(cell_metadata_with_genes, gene_names, 'Class', sort=True)
plot_heatmap(agg, 8, 5)
plt.title('Class')
plt.show()

Visualize gene expression for selected genes by dominant region.

Below we average the expression of our selected genes across the different dominant regions that each cell is uniquely assigned to. We compute these expressions for each species individually.

gene_names = cell_metadata_with_genes.columns[
    cell_metadata_with_genes.columns.isin(gene_names)
].to_numpy()
for level in ['Group', 'Cluster']:
    agg = aggregate_by_metadata(cell_metadata_with_genes, gene_names, f'{level}_dominant_region', sort=True)
    plot_heatmap(agg, 8, 5)
    plt.title(f'{level} Dominant Region')
    plt.show()

Visualize gene expression for selected genes by the HMBA-Subcortex taxonomy.

Below we average the expression of our selected genes across the neurotransmitter at the Group and Cluster level.

gene_names = cell_metadata_with_genes.columns[
    cell_metadata_with_genes.columns.isin(gene_names)
].to_numpy()
for level in ['Group', 'Cluster']:
    agg = aggregate_by_metadata(cell_metadata_with_genes, gene_names, f'{level}_neurotransmitter', sort=True)
    plot_heatmap(agg, 8, 5)
    plt.title(f'{level} Neurotransmitter')
    plt.show()

To learn more about the spatial data see the HBMA-Subcortex Spatial Slabs and Taxonomy notebook.

To learn more about the taxonomy used in this notebook and the data it is derived from, see the HMBA-Subcortex Clustering Analysis and Taxonomy notebook.