Human and Mammalian Brain Atlas (HMBA) - Marmoset Subcortex 10X snRNA-seq: Gene Expression

Human and Mammalian Brain Atlas (HMBA) - Marmoset Subcortex 10X snRNA-seq: 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).

The dataset includes the gene expression data from the snRNA-seq component of ~680k high-quality nuclei collected via 10x Multiome (paired snRNA-seq + ATAC-seq). The entire subcortex of four adult marmoset donors (2M, 2F) was sampled via contiguous tiles and resulted in ~680k cells passing rigorous quality control (QC). The data were clustered to resolve 1,737 transcriptomic Clusters. Coupled with the anatomical anchoring provided by the Xenium spatial transcriptomics dataset, these Clusters were hierarchically organized 212 Groups and 17 Classes to form a reference cell type taxonomy of the non-human primate subcortex.

For more information on this dataset and the associated Xenium spatial transcriptomics 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’ll explore gene expressions and combine them with the cell metadata we showed in the previous clustering analysis tutorial.

You need to be connected to the internet to run this notebook or connected to a cache that has the HMBA-Subcortex-10xMultiome-Marmoset data downloaded already.

The notebook presented here shows quick visualizations from precomputed metadata in the atlas. For examples on accessing the expression matrices, specifically selecting genes from expression matrices, see the general Accessing expression data tutorial.

%matplotlib inline

import pandas as pd
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
from typing import List, Optional, Tuple

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 up a Pandas DataFrame. See the getting_started notebook for more details on using the cache including installing it if it has not already been.

Change the download_base variable to where you would like to download the data to your system or a location where a cache is already available.

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'

Below we create the expanded cell metadata as was done previously in the cluster annotation tutorial.

cell_to_cluster_membership = abc_cache.get_metadata_dataframe(
    'HMBA-Subcortex-taxonomy-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')

cell_metadata = abc_cache.get_metadata_dataframe('HMBA-Subcortex-10xMultiome-Marmoset', 'cell_metadata').set_index('cell_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']]

cell_extended = cell_metadata.join(cell_to_cluster_membership, how='inner')
cell_extended = cell_extended.join(cluster_details, on='cluster_alias')
cell_extended = cell_extended.join(cluster_colors, on='cluster_alias', rsuffix='_color')
cell_extended = cell_extended.join(cluster_order, on='cluster_alias', rsuffix='_order')

barcoded_cell_samples = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-10xMultiome-Marmoset',
    file_name='barcoded_cell_sample'
).set_index('barcoded_cell_sample_label')
donor = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-10xMultiome-Marmoset',
    file_name='donor'
).set_index('donor_label')

cell_extended = cell_extended.join(barcoded_cell_samples, on='barcoded_cell_sample_label', rsuffix='_bc')
cell_extended = cell_extended.join(donor, on='donor_label', rsuffix='_donor')

value_sets = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-10xMultiome-Marmoset',
    file_name='value_sets'
).set_index('label')

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_extended = cell_extended.join(cluster_to_dominant_region, on='cluster_alias')
cell_extended = cell_extended.join(cluster_to_neurotransmitter, on='cluster_alias')

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

# Add region of interest color and order
extract_value_set(cell_extended, value_sets, 'region_of_interest_label')
# Add species common name color and order
extract_value_set(cell_extended, value_sets, 'species_genus')
# Add species scientific name color and order
extract_value_set(cell_extended, value_sets, 'species_scientific_name')
# Add donor sex color and order
extract_value_set(cell_extended, value_sets, 'donor_sex')
extract_value_set(cell_extended, value_sets, 'Cluster_neurotransmitter', 'neurotransmitter')
extract_value_set(cell_extended, value_sets, 'Group_neurotransmitter', 'neurotransmitter')
extract_value_set(cell_extended, value_sets, 'Cluster_dominant_region', 'dominant_region')
extract_value_set(cell_extended, value_sets, 'Group_dominant_region', 'dominant_region')
cell_extended.head()
barcoded_cell_sample_label nCount_RNA nFeature_RNA doublet_score donor_label alignment_job_id In_BG_taxonomy_dataset feature_matrix_label dataset_label cluster_alias ... donor_sex_color donor_sex_order Cluster_neurotransmitter_color Cluster_neurotransmitter_order Group_neurotransmitter_color Group_neurotransmitter_order Cluster_dominant_region_color Cluster_dominant_region_order Group_dominant_region_color Group_dominant_region_order
cell_label
GTTCCTGGTTAGCATG-P0005_1 P0005_1 33732 6363 0.118006 CJ23.56.002 240822 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_1 ... #ADC4C3 2 #E0E0E0 10 #E0E0E0 10 #E0E0E0 17 #E0E0E0 17
TACTGAGGTTAGGATT-P0005_1 P0005_1 29309 6123 0.520672 CJ23.56.002 240822 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_17 ... #ADC4C3 2 #ef3a5b 2 #ef3a5b 2 #f06597 10 #f06597 10
TCTTCAAGTTTATGGG-P0005_1 P0005_1 28284 5758 0.125381 CJ23.56.002 240822 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_42 ... #ADC4C3 2 #ef3a5b 2 #ef3a5b 2 #f06597 10 #f06597 10
CCCTCAGTCAGCACCA-P0005_1 P0005_1 25495 5729 0.225403 CJ23.56.002 240822 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_41 ... #ADC4C3 2 #ef3a5b 2 #ef3a5b 2 #f06597 10 #f06597 10
AACTACTCACCTACTT-P0005_1 P0005_1 23770 5132 0.557602 CJ23.56.002 240822 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_18 ... #ADC4C3 2 #ef3a5b 2 #ef3a5b 2 #f06597 10 #f06597 10

5 rows × 62 columns

Single cell/nucleus transcriptomes#

Below we use the convenience function get_gene_data to download and extract specific genes from the gene expression h5ad files. This function can be used to pull expression for the full set of cells or any subset from the set of cell metadata for specific genes. See accessing gene expression data tutorial for more information.

We first load the set of genes for these data.

gene = abc_cache.get_metadata_dataframe(
    directory='HMBA-10xMultiome-BG',
    file_name='marmoset_gene'
).set_index('gene_identifier')
gene.head(5)

marmoset_gene.csv: 100%|██████████| 2.61M/2.61M [00:00<00:00, 3.15MMB/s]
gene_symbol description molecular_type
gene_identifier
NCBIGene:118152095 LOC118152095 glucose-6-phosphate 1-dehydrogenase-like transcribed_pseudogene
NCBIGene:118147983 LOC118147983 U6 spliceosomal RNA snRNA
NCBIGene:100395011 LOC100395011 alpha-internexin-like pseudogene
NCBIGene:118147921 LOC118147921 U4 spliceosomal RNA snRNA
NCBIGene:100385748 LOC100385748 heterogeneous nuclear ribonucleoprotein A1-like protein_coding

Below we list the genes we will use in this notebook and the example method used to load the expression for these specific genes from the h5ad file.

To process and extract the gene expressions for yourself, uncomment the code block below. Warning that this is a large amount of data and may take a significant fraction of time to download. This download action will only need to be performed once, however.

For more details on how to extract specific genes from the data see our accessing gene expression data tutorial

gene_names = ['APOE', 'DRD1', 'DRD2', 'AQP4', 'FLT1', 'MME', 'RELN', 'LRP8']

"""
gene_data = get_gene_data(
    abc_atlas_cache=abc_cache,
    all_cells=cell_metadata,
    all_genes=gene,
    selected_genes=gene_names
)"""

'\ngene_data = get_gene_data(\n    abc_atlas_cache=abc_cache,\n    all_cells=cell_metadata,\n    all_genes=gene,\n    selected_genes=gene_names\n)'

Instead of processing the gene expressions, we load pre-processed files containing the log gene expression for the genes listed above. The expression in the table below are presented as log2(CPM + 1).

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

example_gene_expression.csv: 100%|██████████| 43.7M/43.7M [00:07<00:00, 5.60MMB/s]
DRD1 FLT1 LRP8 RELN DRD2 AQP4 MME APOE
cell_label
GTTCCTGGTTAGCATG-P0005_1 0.000000 0.0 4.937601 0.000000 0.0 0.0 0.000000 0.0
TACTGAGGTTAGGATT-P0005_1 0.000000 0.0 6.691501 9.264785 0.0 0.0 0.000000 0.0
TCTTCAAGTTTATGGG-P0005_1 0.000000 0.0 0.000000 0.000000 0.0 0.0 0.000000 0.0
CCCTCAGTCAGCACCA-P0005_1 0.000000 0.0 6.890813 0.000000 0.0 0.0 5.329962 0.0
AACTACTCACCTACTT-P0005_1 5.428606 0.0 5.428606 8.982532 0.0 0.0 6.411760 0.0

We load and merge the expression into each of our cell metadata.

cell_extended = cell_extended.join(gene_data)

Example use cases#

In this section, we show a use case plotting our sets of genes. First we’ll show the gene expression in heatmaps plotted disease metadata and the taxonomy. We’ll then plot the genes in a UMAP.

Heatmap of Average gene expression#

The helper function below creates a heatmap showing the relation between various parameters in the combined cell metadata and the genes we loaded.

import matplotlib as mpl

def plot_heatmap(
    df: pd.DataFrame,
    gnames: List[str],
    value: str,
    sort: bool = False,
    fig_width: float = 8,
    fig_height: float = 4,
    vmax: float = None,
    cmap: plt.cm = plt.cm.magma,
):
    """Plot a heatmap of gene expression values for a list of genes across species.

    Parameters
    ----------
    df : pd.DataFrame
        DataFrame containing cell metadata and gene expression values.
    gnames : list
        List of gene names to plot.
    value : str
        Column name in df to group by (e.g., 'species_genus').
    sort : bool, optional
        Whether to sort the gene expression values within each species.
    fig_width : float, optional
        Width of the figure in inches.
    fig_height : float, optional
        Height of the figure in inches.
    vmax : float, optional
        Maximum value for the color scale. If None, it is set to the maximum value in the data.
    cmap : matplotlib colormap, optional
        Colormap to use for the heatmap.

    Returns
    -------
    fig : matplotlib.figure.Figure
        The figure object containing the heatmap.
    ax : array of matplotlib.axes.Axes
        The axes objects for each species.
    """

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

    grouped = df.groupby(value)[gnames].mean()
    vmin = grouped.min().min()
    vmax = grouped.max().max()
    norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    cmap = sm.get_cmap()

    if sort:
        grouped = grouped.sort_values(by=gnames[0], ascending=False)
        grouped = grouped.loc[sorted(grouped.index)]

    value_order = df.groupby(value)[[f'{value}_order']].first()
    grouped = grouped.loc[
        value_order.sort_values(f'{value}_order').index
    ]

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

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


    ax.set_yticks(range(len(ylabs)))
    ax.set_yticklabels(ylabs)
    ax.set_xticks(range(len(xlabs)))
    ax.set_xticklabels(xlabs)

    cbar = fig.colorbar(sm, ax=ax, orientation='vertical', fraction=0.01, pad=0.01)
    cbar.set_label('Mean Expression [log2(CPM + 1)]')
    plt.subplots_adjust(wspace=0.00, hspace=0.00)
    
    return fig, ax

We now use the above function to show the average expression across Class, Group, and Group.

fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='Class',
    fig_width=15,
    fig_height=5
)
fig.suptitle('Class')
plt.show()
../_images/c521793545a9c61cef2cbb441de24c6f1c6895f63a68dbc67a246f1d78a9e04f.png 
fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='Group',
    fig_width=15,
    fig_height=30,
    sort=True
)
fig.suptitle('Group')
plt.show()
../_images/5a2133093e916ff19e5f64b3834747aebce86b95497ca41283db462671c0815c.png

Next wee plot the expression against the dominant Group region.

fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='Group_dominant_region',
    fig_width=15,
    fig_height=5,
    sort=True
)
fig.suptitle('Group Dominant Region')
plt.show()
../_images/d1f94df038933f4abee46f0a7bca8f40596f0d56a9cce5a89c72c9ccce7c5231.png

Expression in the UMAP#

We load the UMAP coordinates for our cells and plot the expression in the UMAP for each of our selected genes. The UMAP is the same as the previous notebook.

cell_2d_embedding_coordinates = abc_cache.get_metadata_dataframe(
    directory='HMBA-Subcortex-taxonomy-Marmoset',
    file_name='cell_2d_embedding_coordinates'
).set_index('cell_label')
cell_2d_embedding_coordinates.head()
x y
cell_label
GTTCCTGGTTAGCATG-P0005_1 14.662888 10.134843
TACTGAGGTTAGGATT-P0005_1 -8.280969 18.335058
TCTTCAAGTTTATGGG-P0005_1 7.183354 13.461928
CCCTCAGTCAGCACCA-P0005_1 7.140094 13.440581
AACTACTCACCTACTT-P0005_1 -8.869915 18.777521

Join the coordinates into the cell metadata.

cell_extended = cell_extended.join(
    cell_2d_embedding_coordinates,
    how='inner'
).sample(frac=1)
cell_extended.head()
barcoded_cell_sample_label nCount_RNA nFeature_RNA doublet_score donor_label alignment_job_id In_BG_taxonomy_dataset feature_matrix_label dataset_label cluster_alias ... DRD1 FLT1 LRP8 RELN DRD2 AQP4 MME APOE x y
cell_label
TTACAGCAGCCTGATG-P0042_2 P0042_2 18368 5296 0.064197 CJ23.56.003 240807 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_37 ... 0.0 0.0 6.779851 8.577797 0.000000 0.0 0.00000 0.0 -7.590099 18.316439
TATTTGCTCGCACACA-P0076_2 P0076_2 33623 6449 0.251674 CJ24.56.004 250529 0 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_54 ... 0.0 0.0 5.918460 0.000000 8.358242 0.0 6.49545 0.0 5.548254 16.739714
GTGCTGATCGTGCTTC-P0019_3 P0019_3 11071 3819 0.070829 CJ23.56.003 240110 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_1741 ... 0.0 0.0 7.505034 0.000000 0.000000 0.0 0.00000 0.0 1.108928 -7.322788
GACCTAAGTCCTCCAA-P0024_1 P0024_1 3057 1588 0.125816 CJ23.56.003 240323 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_7 ... 0.0 0.0 0.000000 0.000000 0.000000 0.0 0.00000 0.0 11.842602 11.607631
ATTGCAGGTGCTCCGT-P0052_2 P0052_2 9522 3654 0.052085 CJ24.56.001 241212 1 HMBA-Subcortex-10xMultiome-Marmoset HMBA-Subcortex-10xMultiome-Marmoset cluster_1689 ... 0.0 0.0 8.304054 0.000000 0.000000 0.0 0.00000 0.0 3.154555 11.765785

5 rows × 72 columns

We use the same UMAP plotting convenience function as previously used in the HMBA-Subcortex clustering and annotation tutorial.

def plot_umap(
    xx: np.ndarray,
    yy: np.ndarray,
    cc: np.ndarray = None,
    val: np.ndarray = None,
    fig_width: float = 8,
    fig_height: float = 8,
    cmap: Optional[plt.Colormap] = None,
    labels: np.ndarray = None,
    term_orders: np.ndarray = None,
    colorbar: bool = False,
    sizes: np.ndarray = None,
    fig: plt.Figure = None,
    ax: plt.Axes = None,
 ) -> Tuple[plt.Figure, plt.Axes]:
    """
    Plot a scatter plot of the UMAP coordinates.

    Parameters
    ----------
    xx : array-like
        x-coordinates of the points to plot.
    yy : array-like
        y-coordinates of the points to plot.
    cc : array-like, optional
        colors of the points to plot. If None, the points will be colored by the values in `val`.
    val : array-like, optional
        values of the points to plot. If None, the points will be colored by the values in `cc`.
    fig_width : float, optional
        width of the figure in inches. Default is 8.
    fig_height : float, optional
        height of the figure in inches. Default is 8.
    cmap : str, optional
        colormap to use for coloring the points. If None, the points will be colored by the values in `cc`.
    labels : array-like, optional
        labels for the points to plot. If None, no labels will be added to the plot.
    term_orders : array-like, optional
        order of the labels for the legend. If None, the labels will be ordered by their appearance in `labels`.
    colorbar : bool, optional
        whether to add a colorbar to the plot. Default is False.
    sizes : array-like, optional
        sizes of the points to plot. If None, all points will have the same size.
    fig : matplotlib.figure.Figure, optional
        figure to plot on. If None, a new figure will be created with 1 subplot.
    ax : matplotlib.axes.Axes, optional
        axes to plot on. If None, a new figure will be created with 1 subplot.
    """

    if sizes is None:
        sizes = 1
    if ax is None or fig is None:
        fig, ax = plt.subplots()
        fig.set_size_inches(fig_width, fig_height)

    if cmap is not None:
        scatt = ax.scatter(xx, yy, c=val, s=0.5, marker='.', cmap=cmap, alpha=sizes)
    elif cc is not None:
        scatt = ax.scatter(xx, yy, c=cc, s=0.5, marker='.', alpha=sizes)

    if labels is not None:
        from matplotlib.patches import Rectangle
        unique_label_colors = (labels + ',' + cc).unique()
        unique_labels = np.array([label_color.split(',')[0] for label_color in unique_label_colors])
        unique_colors = np.array([label_color.split(',')[1] for label_color in unique_label_colors])

        if term_orders is not None:
            unique_order = term_orders.unique()
            term_order = np.argsort(unique_order)
            unique_labels = unique_labels[term_order]
            unique_colors = unique_colors[term_order]
            
        rects = []
        for color in unique_colors:
            rects.append(Rectangle((0, 0), 1, 1, fc=color))

        legend = ax.legend(rects, unique_labels, loc=0)
        # ax.add_artist(legend)
    
    ax.set_xticks([])
    ax.set_yticks([])

    if colorbar:
        fig.colorbar(scatt, ax=ax)
    
    return fig, ax

Below we plot the genes in our set next to the UMAP colored by one of our features.

term_to_plot = 'Class'

# Plot UMAPs for the first two genes in gene_names
for gene_name in gene_names:
    fig, ax = plt.subplots(1, 2)
    fig.set_size_inches(18, 9)
    ax = ax.flatten()
    plot_umap(
        cell_extended['x'],
        cell_extended['y'],
        cc=cell_extended[term_to_plot + '_color'],
        labels=cell_extended[term_to_plot],
        term_orders=cell_extended[term_to_plot + '_order'],
        fig=fig,
        ax=ax[0],
    )
    ax[0].set_title(term_to_plot)
    plot_umap(
        cell_extended['x'],
        cell_extended['y'],
        val=cell_extended[gene_name],
        cmap=plt.cm.magma,
        fig=fig,
        ax=ax[1],
        colorbar=False,
    )
    ax[1].set_title(f'Gene Expression: {gene_name}')
    plt.tight_layout()
    plt.show()
../_images/a6934fc8be2d2f4b2e1ae85652e5f4c1ad43d12e424725b186c5a05b9865ac71.png ../_images/464ced525533c5d4dec5266f81fe6c0056ad821e1a14007e2dc3ad8e7a8dee39.png ../_images/3b9954580b0fff5f36ecdbbcd8ab7722e5ba14d4c7540119bf3e434dd95f4db5.png ../_images/a3cf0a2d67b68e107d66b2dec8ebc1f1f99223de73dc937b03d0bd6d77c31ceb.png ../_images/a5adabf3542bb1d10753b266953eec3a99e47055bce6d1d2eeed4f0f5aff769c.png ../_images/4180c7e18d3252c1ebf7ea0023395ec790dbf8f1f36876e60dc0803406047475.png ../_images/7fcfaf76b7ad1565fd887c929985a3bbe41a49b4de4adebf47f90136db461056.png ../_images/ae91ab98ba1b0aa7b69dba9325eaf4687145ea34399fa16b069222799047edaf.png

Also available is the HMBA-Subcortex taxonomy notebook.