SEA-AD Mutlitregion taxonomy and snRNASeq analysis: 10X gene expression

In this notebook we’ll explore some 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 SEA-AD Multiregion 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/allen-brain-cell-atlas-staging/')
abc_cache = AbcProjectCache.from_cache_dir(
   download_base
)

abc_cache.current_manifest

'releases/20260711/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(
    'SEA-AD-Multiregion-taxonomy', 'cell_to_cluster_membership'
).set_index('cell_label')
cluster = abc_cache.get_metadata_dataframe(
    'SEA-AD-Multiregion-taxonomy', 'cluster'
).set_index('label')
cluster_annotation_term = abc_cache.get_metadata_dataframe(
    'SEA-AD-Multiregion-taxonomy', 'cluster_annotation_term'
).set_index('label')

cluster_annotation_term_set = abc_cache.get_metadata_dataframe(
    directory='SEA-AD-Multiregion-taxonomy',
    file_name='cluster_annotation_term_set'
).set_index('label')

cluster_to_cluster_annotation_membership = abc_cache.get_metadata_dataframe(
    directory='SEA-AD-Multiregion-taxonomy',
    file_name='cluster_to_cluster_annotation_membership'
).set_index('cluster_annotation_term_label')

cell_metadata = abc_cache.get_metadata_dataframe('SEA-AD-Multiregion-10X', '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')

library = abc_cache.get_metadata_dataframe(
    directory='SEA-AD-Multiregion-10X',
    file_name='library'
).set_index('library_label')
donor = abc_cache.get_metadata_dataframe(
    directory='SEA-AD-Multiregion-10X',
    file_name='donor'
).set_index('donor_label')
disease = abc_cache.get_metadata_dataframe(
    directory='SEA-AD-Multiregion-10X',
    file_name='disease'
).set_index('donor_label')

cell_extended = cell_extended.join(library, on='library_label', rsuffix='_lib')
cell_extended = cell_extended.join(donor, on='donor_label', rsuffix='_donor')
cell_extended = cell_extended.join(disease, on='donor_label', rsuffix='_disease')

value_sets = abc_cache.get_metadata_dataframe(
    directory='SEA-AD-Multiregion-10X',
    file_name='value_sets'
).set_index('label')

def extract_value_set(
        cell_metadata_df: pd.DataFrame,
        input_value_set: pd.DataFrame,
        input_value_set_label: str,
        dataframe_column: Optional[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 dataframe_column is None:
        dataframe_column = input_value_set_label
    cell_metadata_df[f'{dataframe_column}_color'] = input_value_set[
        input_value_set['field'] == input_value_set_label
    ].loc[cell_metadata_df[dataframe_column]]['color_hex_triplet'].values
    cell_metadata_df[f'{dataframe_column}_order'] = input_value_set[
        input_value_set['field'] == input_value_set_label
    ].loc[cell_metadata_df[dataframe_column]]['order'].values

# Add region of interest color and order
extract_value_set(cell_extended, value_sets, 'region_of_interest_label')
# Add disease color and order
extract_value_set(cell_extended, value_sets, 'APOE Genotype')
extract_value_set(cell_extended, value_sets, 'Arteriolosclerosis')
extract_value_set(cell_extended, value_sets, 'Atherosclerosis')
extract_value_set(cell_extended, value_sets, 'Braak')
extract_value_set(cell_extended, value_sets, 'CERAD score')
extract_value_set(cell_extended, value_sets, 'Cognitive Status')
extract_value_set(cell_extended, value_sets, 'LATE')
extract_value_set(cell_extended, value_sets, 'Highest Lewy Body Disease')
extract_value_set(cell_extended, value_sets, 'Overall AD neuropathological Change')
extract_value_set(cell_extended, value_sets, 'Severely Affected Donor')
extract_value_set(cell_extended, value_sets, 'Thal')
# Add donor sex/gender color and order
extract_value_set(cell_extended, value_sets, 'donor_sex')
extract_value_set(cell_extended, value_sets, 'donor_gender')
# Add race
extract_value_set(cell_extended, value_sets, 'donor_race')
cell_extended.head()

/Users/chris.morrison/src/abc_atlas_access/src/abc_atlas_access/abc_atlas_cache/abc_project_cache.py:643: DtypeWarning: Columns (4) have mixed types. Specify dtype option on import or set low_memory=False.
  return pd.read_csv(path, **kwargs)

	cell_barcode	barcoded_cell_sample_label	library_label	alignment_job_id	doublet_score	umi_count	donor_label	exp_component_name	feature_matrix_label	dataset_label	...	Severely Affected Donor_color	Severely Affected Donor_order	Thal_color	Thal_order	donor_sex_color	donor_sex_order	donor_gender_color	donor_gender_order	donor_race_color	donor_race_order
cell_label
AAACAGCCACTGGCTG-2001_A08	AAACAGCCACTGGCTG	2001_A08	L8XR_231221_02_D02	1322484698	0.184615	54020.0	H21.33.001	AAACAGCCACTGGCTG-L8XR_231221_02_D02-1322484698	AnG-10X	SEA-AD-Multiregion-10X	...	#fb6a4a	1	#fc8161	3	#ADC4C3	2	#ADC4C3	2	#f7f184	125
AAACAGCCAGGTTATT-2001_A08	AAACAGCCAGGTTATT	2001_A08	L8XR_231221_02_D02	1322484698	0.060606	4438.0	H21.33.001	AAACAGCCAGGTTATT-L8XR_231221_02_D02-1322484698	AnG-10X	SEA-AD-Multiregion-10X	...	#fb6a4a	1	#fc8161	3	#ADC4C3	2	#ADC4C3	2	#f7f184	125
AAACAGCCATTCCTCG-2001_A08	AAACAGCCATTCCTCG	2001_A08	L8XR_231221_02_D02	1322484698	0.276923	66285.0	H21.33.001	AAACAGCCATTCCTCG-L8XR_231221_02_D02-1322484698	AnG-10X	SEA-AD-Multiregion-10X	...	#fb6a4a	1	#fc8161	3	#ADC4C3	2	#ADC4C3	2	#f7f184	125
AAACATGCAATGCCTA-2001_A08	AAACATGCAATGCCTA	2001_A08	L8XR_231221_02_D02	1322484698	0.060606	6019.0	H21.33.001	AAACATGCAATGCCTA-L8XR_231221_02_D02-1322484698	AnG-10X	SEA-AD-Multiregion-10X	...	#fb6a4a	1	#fc8161	3	#ADC4C3	2	#ADC4C3	2	#f7f184	125
AAACATGCACCTCACC-2001_A08	AAACATGCACCTCACC	2001_A08	L8XR_231221_02_D02	1322484698	0.181818	48653.0	H21.33.001	AAACATGCACCTCACC-L8XR_231221_02_D02-1322484698	AnG-10X	SEA-AD-Multiregion-10X	...	#fb6a4a	1	#fc8161	3	#ADC4C3	2	#ADC4C3	2	#f7f184	125

5 rows × 105 columns

Single cell and nuclei 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='SEA-AD-Multiregion-10X',
    file_name='gene'
).set_index('gene_identifier')
gene.head(5)

gene.csv: 100%|██████████| 1.74M/1.74M [00:00<00:00, 3.91MMB/s]

	gene_symbol	feature_types	genome
gene_identifier
ENSG00000243485	MIR1302-2HG	Gene Expression	GRCh38
ENSG00000237613	FAM138A	Gene Expression	GRCh38
ENSG00000186092	OR4F5	Gene Expression	GRCh38
ENSG00000238009	AL627309.1	Gene Expression	GRCh38
ENSG00000239945	AL627309.3	Gene Expression	GRCh38

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', 'SLC17A7', 'GAD1', 'SST', '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='SEA-AD-Multiregion-10X',
    file_name='example_gene_expression'
).set_index('cell_label')
gene_data.head()

example_gene_expression.csv: 100%|██████████| 462M/462M [00:46<00:00, 9.99MMB/s]    

	LRP8	GAD1	MME	SST	RELN	APOE	SLC17A7
cell_label
AAACAGCCACTGGCTG-2001_A08	5.248813	4.286265	0.0	0.0	0.000000	0.000000	0.000000
AAACAGCCAGGTTATT-2001_A08	0.000000	0.000000	0.0	0.0	7.822263	7.822263	0.000000
AAACAGCCATTCCTCG-2001_A08	0.000000	0.000000	0.0	0.0	0.000000	0.000000	6.927078
AAACATGCAATGCCTA-2001_A08	8.380596	8.380596	0.0	0.0	0.000000	0.000000	0.000000
AAACATGCACCTCACC-2001_A08	7.539031	10.470869	0.0	0.0	0.000000	0.000000	0.000000

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 against ROI, disease metadata, and 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 and subclass.

fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='Class',
    fig_width=15,
    fig_height=5
)
fig.suptitle('Class')
plt.show()

fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='Subclass',
    fig_width=15,
    fig_height=10,
    sort=True
)
fig.suptitle('Subclass')
plt.show()

Nex, we break down the expression by region of interest.

fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='region_of_interest_label',
    fig_width=15,
    fig_height=5,
    sort=True,
)
fig.suptitle('Brain Region')
plt.show()

Finally we show the average expression in the heatmap vs one of the disease markers. Here we use Thal Phase

fig, ax = plot_heatmap(
    df=cell_extended,
    gnames=gene_names,
    value='Thal',
    fig_width=15,
    fig_height=5,
    sort=True,
)
fig.suptitle('Thal Phase')
plt.show()

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='SEA-AD-Multiregion-taxonomy',
    file_name='cell_2d_embedding_coordinates'
).set_index('cell_label')
cell_2d_embedding_coordinates.head()

	x	y
cell_label
AAACAGCCACTGGCTG-2001_A08	7.009472	14.111731
AAACAGCCAGGTTATT-2001_A08	-1.319671	-1.221017
AAACAGCCATTCCTCG-2001_A08	17.003418	10.816355
AAACATGCAATGCCTA-2001_A08	17.077744	2.317490
AAACATGCACCTCACC-2001_A08	19.195032	1.132652

Join the coordinates into the cell metadata.

cell_extended = cell_extended.join(
    cell_2d_embedding_coordinates,
    how='inner'
).sample(frac=1)
cell_extended.head()

	cell_barcode	barcoded_cell_sample_label	library_label	alignment_job_id	doublet_score	umi_count	donor_label	exp_component_name	feature_matrix_label	dataset_label	...	donor_race_order	LRP8	GAD1	MME	SST	RELN	APOE	SLC17A7	x	y
cell_label
TTGTTTGTCGATTGAC-669_F01	TTGTTTGTCGATTGAC	669_F01	L8TX_210603_01_A10	1109152711	0.090000	39340.0	H20.33.014	TTGTTTGTCGATTGAC-L8TX_210603_01_A10-1109152711	DFC-10X	SEA-AD-Multiregion-10X	...	125	4.723527	0.0	0.0	0.0	0.0	0.000000	7.001094	6.068086	17.445227
ACTAATCCAAGATTCT-1438_A02	ACTAATCCAAGATTCT	1438_A02	L8XR_220929_01_E10	1219250383	0.058824	64995.0	H21.33.018	ACTAATCCAAGATTCT-L8XR_220929_01_E10-1219250383	HIP-10X	SEA-AD-Multiregion-10X	...	125	0.000000	0.0	0.0	0.0	0.0	4.034374	7.411459	5.655291	8.783364
GTCTCACCAAGAATAC-768_G01	GTCTCACCAAGAATAC	768_G01	L8TX_210812_01_A11	1124416554	0.043478	9413.0	H20.33.040	GTCTCACCAAGAATAC-L8TX_210812_01_A11-1124416554	DFC-10X	SEA-AD-Multiregion-10X	...	125	0.000000	0.0	0.0	0.0	0.0	6.744646	0.000000	19.601860	16.896631
GAGAGGTAGCACCTGC-1096_B04	GAGAGGTAGCACCTGC	1096_B04	L8TX_220224_02_G07	1162262296	0.120482	53524.0	H20.33.029	GAGAGGTAGCACCTGC-L8TX_220224_02_G07-1162262296	MEC-10X	SEA-AD-Multiregion-10X	...	65	6.242847	0.0	0.0	0.0	0.0	0.000000	7.042015	18.256517	14.536785
AGTCAACTCCCATTTA-644_D03	AGTCAACTCCCATTTA	644_D03	L8TX_210514_01_G12	1142430439	0.096386	55228.0	H20.33.004	AGTCAACTCCCATTTA-L8TX_210514_01_G12-1142430439	MTG-10X	SEA-AD-Multiregion-10X	...	125	0.000000	0.0	0.0	0.0	0.0	0.000000	6.516233	11.829694	8.897827

5 rows × 114 columns

We use the same UMAP plotting convenience function as previously used in the SEA-AD Multiregion 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. We select every 10th cell via [::10] for ease of plotting. These can be removed to plot the full dataset.

term_to_plot = 'region_of_interest_label'

# 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'][::10],
        cell_extended['y'][::10],
        cc=cell_extended[term_to_plot + '_color'][::10],
        labels=cell_extended[term_to_plot][::10],
        term_orders=cell_extended[term_to_plot + '_order'][::10],
        fig=fig,
        ax=ax[0],
    )
    ax[0].set_title(term_to_plot)
    plot_umap(
        cell_extended['x'][::10],
        cell_extended['y'][::10],
        val=cell_extended[gene_name][::10],
        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()

Also available is the SEA-AD Multiregion taxonomy notebook. You can also investigate the SEA-AD Caudate data and taxonomy here.