Human-Mammalian Brain - Basal Ganglia spatial transcriptomic slab coordinates and annotation

The basal ganglia are a set of subcortical structures critical for motor control, particularly in the context of action selection, motor learning and emotional state, whose coarse functional organization is well-described in the literature. Although single-cell and single-nucleus RNA-seq studies have greatly enhanced our understanding of cellular diversity in the brain, the exact spatial distribution of cell types within the brain has been difficult to determine. Using spatial transcriptomics platforms such as MERSCOPE or Xenium enables spatial profiling of hundreds of genes for each cell, enabling mapping of the cells to the consensus basal ganglia cell type taxonomy.

In this notebook, we walk through the data available in this release and show how they are joined together to enable analyses at the cell level.

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

from abc_atlas_access.abc_atlas_cache.abc_project_cache import AbcProjectCache

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 in your system.

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

abc_cache.current_manifest

'releases/20261101/manifest.json'

Data Overview

Below we list the available metadata and expression matrix files available for each species.

The various species in this dataset are divided up by their donor_id and transcriptomic platform. MERSCOPE for Human and Macaque, Xenium for Marmoset.

The mapping between donor to species is:

• Human: H22.30.001

• Macaque: QM23.50.001

• Marmoset: CJ23.56.004

for dataset_name in ['HMBA-BG-MERSCOPE-Human', 'HMBA-BG-MERSCOPE-Macaque', 'HMBA-BG-Xenium-Marmoset']:
    print(f"Dataset: {dataset_name}")
    metadata_files = abc_cache.list_metadata_files(dataset_name)
    print(f"Metadata files: {metadata_files}")
    expression_matrix_files = abc_cache.list_expression_matrix_files(dataset_name)
    print(f"Expression matrix files: {expression_matrix_files}\n")

Dataset: HMBA-BG-MERSCOPE-Human
Metadata files: ['additional_cell_metadata', 'cell_anatomical_annotations', 'cell_metadata', 'cirro_spatial_grid', 'donor', 'example_gene_expression', 'gene', 'mmc_results', 'slab_plane_coordinates', 'specimen_metadata', 'value_sets', 'visualization']
Expression matrix files: ['MERSCOPE-H22.30.001-BG-batch-normalized/ln', 'MERSCOPE-H22.30.001-BG-batch-normalized/raw', 'MERSCOPE-H22.30.001-BG/raw']

Dataset: HMBA-BG-MERSCOPE-Macaque
Metadata files: ['additional_cell_metadata', 'cell_anatomical_annotations', 'cell_metadata', 'cirro_spatial_grid', 'donor', 'example_gene_expression', 'gene', 'mmc_results', 'slab_plane_coordinates', 'specimen_metadata', 'value_sets', 'visualization']
Expression matrix files: ['MERSCOPE-QM23.50.001-BG-batch-normalized/ln', 'MERSCOPE-QM23.50.001-BG-batch-normalized/raw', 'MERSCOPE-QM23.50.001-BG/raw']

Dataset: HMBA-BG-Xenium-Marmoset
Metadata files: ['additional_cell_metadata', 'ccf_coordinates', 'cell_anatomical_annotations', 'cell_metadata', 'cirro_spatial_grid', 'donor', 'donor_coordinates', 'example_gene_expression', 'gene', 'mmc_results', 'rotated_coordinates', 'slab_plane_coordinates', 'specimen_metadata', 'value_sets', 'visualization']
Expression matrix files: ['Xenium-CJ23.56.004-BG-batch-normalized/ln', 'Xenium-CJ23.56.004-BG-batch-normalized/raw', 'Xenium-CJ23.56.004-BG/raw']

Cell metadata

Each species in the HMBA spatial dataset has it’s own h5ad files and associated cell metadata tables. We load each of these in turn starting from Human to Macaque to Marmoset. 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. Below we load the cell metadata for each species.

Load the Human cells.

human_cell_metadata = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Human',
    file_name='cell_metadata',
    dtype={"cell_label": str},
).set_index('cell_label')
human_cell_metadata

cell_metadata.csv: 100%|██████████| 1.16G/1.16G [00:50<00:00, 22.8MMB/s]   

	brain_section_label	brain_section_barcode	x_experiment	y_experiment	segmentation_job_id	cell_volume	total_counts	total_counts_genes	total_counts_genes_batch_normalized	n_unique_genes	qc_pass	abc_atlas_visualization	dataset_label	feature_matrix_label
cell_label
1315611012_1-20250227	H22.30.001.CX.16.02.03.03	1315611012	4687.0664	2342.8290	sis_cellpose_v1	78.062865	26	26	48.419523	13	True	True	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1315611012_2-20250227	H22.30.001.CX.16.02.03.03	1315611012	4939.2110	2343.9695	sis_cellpose_v1	NaN	5	5	9.311447	5	False	False	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1315611012_3-20250227	H22.30.001.CX.16.02.03.03	1315611012	4788.4240	2343.7327	sis_cellpose_v1	8.110529	11	11	20.485183	7	True	True	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1315611012_4-20250227	H22.30.001.CX.16.02.03.03	1315611012	4974.0500	2352.2866	sis_cellpose_v1	NaN	4	4	7.449157	4	False	False	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1315611012_5-20250227	H22.30.001.CX.16.02.03.03	1315611012	4681.3800	2346.6519	sis_cellpose_v1	NaN	5	5	9.311447	5	False	False	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
1384732733_104386-20250227	H22.30.001.CX.20.02.03.01	1384732733	8903.6710	2691.2349	sis_cellpose_v1	38.255657	46	46	65.389384	26	True	True	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1384732733_104387-20250227	H22.30.001.CX.20.02.03.01	1384732733	8999.2290	2791.0234	sis_cellpose_v1	NaN	1	1	1.421508	1	False	False	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1384732733_104388-20250227	H22.30.001.CX.20.02.03.01	1384732733	8899.7460	2821.3567	sis_cellpose_v1	NaN	3	3	4.264525	3	False	False	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1384732733_104389-20250227	H22.30.001.CX.20.02.03.01	1384732733	9038.0340	2646.0881	sis_cellpose_v1	37.946513	20	20	28.430167	16	True	True	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized
1384732733_104390-20250227	H22.30.001.CX.20.02.03.01	1384732733	9122.9000	2747.7932	sis_cellpose_v1	26.485758	18	18	25.587150	8	True	True	HMBA-BG-MERSCOPE-Human	MERSCOPE-H22.30.001-BG-batch-normalized

5429191 rows × 14 columns

Load the Macaque cells.

macaque_cell_metadata = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Macaque',
    file_name='cell_metadata',
    dtype={"cell_label": str},
).set_index('cell_label')
macaque_cell_metadata.head()

cell_metadata.csv: 100%|██████████| 725M/725M [00:30<00:00, 23.4MMB/s]    

	brain_section_label	brain_section_barcode	x_experiment	y_experiment	segmentation_job_id	cell_volume	total_counts	total_counts_genes	total_counts_genes_batch_normalized	n_unique_genes	qc_pass	abc_atlas_visualization	dataset_label	feature_matrix_label
cell_label
1291814091_147-20250227	QM23.50.001.CX.44.01.05.01	1291814091	4432.6885	3979.6113	sis_cellpose_v1	756.243171	330	330	636.778305	83	True	True	HMBA-BG-MERSCOPE-Macaque	MERSCOPE-QM23.50.001-BG-batch-normalized
1291814091_157-20250227	QM23.50.001.CX.44.01.05.01	1291814091	4430.7207	3986.3145	sis_cellpose_v1	564.354610	326	326	629.059780	69	True	True	HMBA-BG-MERSCOPE-Macaque	MERSCOPE-QM23.50.001-BG-batch-normalized
1291814091_169-20250227	QM23.50.001.CX.44.01.05.01	1291814091	4377.2495	4001.7097	sis_cellpose_v1	167.832120	62	62	119.637136	34	True	True	HMBA-BG-MERSCOPE-Macaque	MERSCOPE-QM23.50.001-BG-batch-normalized
1291814091_180-20250227	QM23.50.001.CX.44.01.05.01	1291814091	4399.9390	4017.4290	sis_cellpose_v1	400.366153	231	231	445.744813	77	True	True	HMBA-BG-MERSCOPE-Macaque	MERSCOPE-QM23.50.001-BG-batch-normalized
1291814091_181-20250227	QM23.50.001.CX.44.01.05.01	1291814091	4392.6035	4016.5780	sis_cellpose_v1	308.688416	203	203	391.715139	70	True	True	HMBA-BG-MERSCOPE-Macaque	MERSCOPE-QM23.50.001-BG-batch-normalized

Load the Marmoset cells.

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

cell_metadata.csv: 100%|██████████| 378M/378M [00:22<00:00, 16.5MMB/s]    

	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	total_counts_genes	total_counts_genes_batch_normalized	n_unique_genes	qc_pass	abc_atlas_visualization	dataset_label	feature_matrix_label
cell_label
1382042951_abmkfncc-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	5923.4927	6160.6543	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	123.005630	72.566090	1	354	354	270.073610	111	True	True	HMBA-BG-Xenium-Marmoset	Xenium-CJ23.56.004-BG-batch-normalized
1382042951_abmmdgea-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	6137.5723	6217.0757	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	46.194847	NaN	0	193	193	147.243522	93	True	True	HMBA-BG-Xenium-Marmoset	Xenium-CJ23.56.004-BG-batch-normalized
1382042951_abmojfce-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	5802.6750	5614.0080	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	52.652190	19.462345	2	313	313	238.793898	122	True	True	HMBA-BG-Xenium-Marmoset	Xenium-CJ23.56.004-BG-batch-normalized
1382042951_abmomnim-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	5805.9424	5599.3257	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	35.492813	14.766094	1	161	161	122.830088	73	True	True	HMBA-BG-Xenium-Marmoset	Xenium-CJ23.56.004-BG-batch-normalized
1382042951_aeiamgga-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	6291.2680	6068.7200	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	167.800630	66.921560	2	603	602	459.277721	137	True	True	HMBA-BG-Xenium-Marmoset	Xenium-CJ23.56.004-BG-batch-normalized

Specimen metadata

In order to associate each cell to the associated specimen in comes from, we load the specimen table for each species. This provides a mapping for each cell to the specimen/tissue/brain section it originates from. It also allows us to associate the cell to its associated slab_plane through the specimen id. We’ll explain more regarding the slab_plane when we load the our cells coordinates in the plane.

Below we load the Human specimen data.

human_specimen = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Human',
    file_name='specimen_metadata',
    dtype={"cell_label": str}
).set_index('brain_section_label')
human_specimen.head()

specimen_metadata.csv: 100%|██████████| 21.4k/21.4k [00:00<00:00, 322kMB/s]

	brain_section_barcode	brain_section_thickness	brain_section_plane_of_section	brain_section_ordinal	brain_section_depth	section_set_label	section_set_ordinal	section_set_base_ordinal	brain_block_label	brain_block_number	...	brain_slab_ordinal	brain_slab_thickness	hemisphere	brain_division_label	brain_division_ordinal	donor_label	slab_plane_ordinal	slab_plane_label	slab_plane_order	slab_z_order
brain_section_label
H22.30.001.BS.3.01.01.02	1327888285	14	transverse (rostral to caudal)	74	1036	H22.30.001.BS.3.01.01	1	72	H22.30.001.BS.3.01	1	...	2	4.0	right	H22.30.001.BS	2	H22.30.001	1	H22.30.001.BS.3.Plane.01	20201	0
H22.30.001.BS.3.01.03.02	1327888566	14	transverse (rostral to caudal)	109	1526	H22.30.001.BS.3.01.03	3	107	H22.30.001.BS.3.01	1	...	2	4.0	right	H22.30.001.BS	2	H22.30.001	3	H22.30.001.BS.3.Plane.03	20203	1
H22.30.001.BS.3.02.01.02	1327890979	14	transverse (rostral to caudal)	37	518	H22.30.001.BS.3.02.01	1	35	H22.30.001.BS.3.02	2	...	2	4.0	right	H22.30.001.BS	2	H22.30.001	1	H22.30.001.BS.3.Plane.01	20201	0
H22.30.001.BS.3.02.03.03	1327891526	14	transverse (rostral to caudal)	146	2044	H22.30.001.BS.3.02.03	3	143	H22.30.001.BS.3.02	2	...	2	4.0	right	H22.30.001.BS	2	H22.30.001	3	H22.30.001.BS.3.Plane.03	20203	1
H22.30.001.CX.13.01.03.02	1304618482	10	coronal (caudal to rostral)	143	1430	H22.30.001.CX.13.01.03	3	141	H22.30.001.CX.13.01	1	...	13	4.0	right	H22.30.001.CX	1	H22.30.001	5	H22.30.001.CX.13.Plane.05	11305	4

5 rows × 21 columns

The Macaque specimen table.

macaque_specimen = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Macaque',
    file_name='specimen_metadata',
    dtype={"cell_label": str}
).set_index('brain_section_label')
macaque_specimen.head()

specimen_metadata.csv: 100%|██████████| 20.1k/20.1k [00:00<00:00, 384kMB/s]

	brain_section_barcode	brain_section_thickness	brain_section_plane_of_section	brain_section_ordinal	brain_section_depth	section_set_label	section_set_ordinal	section_set_base_ordinal	brain_block_label	brain_block_number	...	brain_slab_ordinal	brain_slab_thickness	hemisphere	brain_division_label	brain_division_ordinal	donor_label	slab_plane_ordinal	slab_plane_label	slab_plane_order	slab_z_order
brain_section_label
QM23.50.001.CX.42.01.03.02	1291809823	10	coronal (caudal to rostral)	133	1330	QM23.50.001.CX.42.01.03	3	131	QM23.50.001.CX.42.01	1	...	42	4.0	right	QM23.50.001.CX	1	QM23.50.001	5	QM23.50.001.CX.42.Plane.05	14205	2
QM23.50.001.CX.42.01.05.03	1291809799	10	coronal (caudal to rostral)	234	2340	QM23.50.001.CX.42.01.05	5	231	QM23.50.001.CX.42.01	1	...	42	4.0	right	QM23.50.001.CX	1	QM23.50.001	3	QM23.50.001.CX.42.Plane.03	14203	1
QM23.50.001.CX.42.02.03.02	1291810606	10	coronal (caudal to rostral)	134	1340	QM23.50.001.CX.42.02.03	3	132	QM23.50.001.CX.42.02	2	...	42	4.0	right	QM23.50.001.CX	1	QM23.50.001	5	QM23.50.001.CX.42.Plane.05	14205	2
QM23.50.001.CX.42.02.05.03	1291810584	10	coronal (caudal to rostral)	234	2340	QM23.50.001.CX.42.02.05	5	231	QM23.50.001.CX.42.02	2	...	42	4.0	right	QM23.50.001.CX	1	QM23.50.001	3	QM23.50.001.CX.42.Plane.03	14203	1
QM23.50.001.CX.42.02.07.03	1291810555	10	coronal (caudal to rostral)	334	3340	QM23.50.001.CX.42.02.07	7	331	QM23.50.001.CX.42.02	2	...	42	4.0	right	QM23.50.001.CX	1	QM23.50.001	1	QM23.50.001.CX.42.Plane.01	14201	0

5 rows × 21 columns

And finally the Marmoset specimen table.

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

specimen_metadata.csv: 100%|██████████| 10.6k/10.6k [00:00<00:00, 124kMB/s]

	brain_section_barcode	brain_section_thickness	brain_section_plane_of_section	brain_section_ordinal	brain_section_depth	section_set_label	section_set_ordinal	section_set_base_ordinal	brain_block_label	brain_block_number	...	brain_slab_ordinal	brain_slab_thickness	hemisphere	brain_division_label	brain_division_ordinal	donor_label	slab_plane_ordinal	slab_plane_label	slab_plane_order	slab_z_order
brain_section_label
CJ23.56.004.CX.42.01.06	1382042951	10	coronal (rostral to caudal)	7	70	CJ23.56.004.CX.42.01	1	1	CJ23.56.004.CX.42	42	...	42	5.0	right	CJ23.56.004.CX	1	CJ23.56.004	1	CJ23.56.004.CX.42.Plane.01	14201	1
CJ23.56.004.CX.42.02.06	1382042748	10	coronal (rostral to caudal)	27	270	CJ23.56.004.CX.42.02	2	21	CJ23.56.004.CX.42	42	...	42	5.0	right	CJ23.56.004.CX	1	CJ23.56.004	2	CJ23.56.004.CX.42.Plane.02	14202	2
CJ23.56.004.CX.42.03.07	1382042618	10	coronal (rostral to caudal)	48	480	CJ23.56.004.CX.42.03	3	41	CJ23.56.004.CX.42	42	...	42	5.0	right	CJ23.56.004.CX	1	CJ23.56.004	3	CJ23.56.004.CX.42.Plane.03	14203	3
CJ23.56.004.CX.42.04.06	1382042399	10	coronal (rostral to caudal)	67	670	CJ23.56.004.CX.42.04	4	61	CJ23.56.004.CX.42	42	...	42	5.0	right	CJ23.56.004.CX	1	CJ23.56.004	4	CJ23.56.004.CX.42.Plane.04	14204	4
CJ23.56.004.CX.42.05.06	1382042114	10	coronal (rostral to caudal)	87	870	CJ23.56.004.CX.42.05	5	81	CJ23.56.004.CX.42	42	...	42	5.0	right	CJ23.56.004.CX	1	CJ23.56.004	5	CJ23.56.004.CX.42.Plane.05	14205	5

5 rows × 21 columns

Join the specimen tables with the cell metadata for each species.

macaque_cell_metadata = macaque_cell_metadata.join(macaque_specimen, on='brain_section_label', rsuffix='_specimen')
marmoset_cell_metadata = marmoset_cell_metadata.join(marmoset_specimen, on='brain_section_label', rsuffix='_specimen')
human_cell_metadata = human_cell_metadata.join(human_specimen, on='brain_section_label', rsuffix='_specimen')

# clean up memory
del macaque_specimen
del marmoset_specimen
del human_specimen

Spatial, slab-plane coordinates

We load the slab-plane coordinates for all of our cells. Each cell comes from a brain section that is a slice of a block cut from a cross-sectional slab of brain. The brain sections are spatially registered with sections from the neighboring blocks within the brain slab. These associated brain sections are referred to a slab-plane. This slab-plane coordinate system in millimeters are the coordinates we will plot all of our cells in this notebook.

Again, we load the coordinates for each of the 3 different species.

human_cell_slab_plane_coordinates = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Human',
    file_name='slab_plane_coordinates',
    dtype={"cell_label": str}
).set_index('cell_label')
human_cell_slab_plane_coordinates.head()

slab_plane_coordinates.csv: 100%|██████████| 680M/680M [00:36<00:00, 18.5MMB/s]    

	x_slab_pixel	y_slab_pixel	x_slab_mm	y_slab_mm	slab_plane_label
cell_label
1315611012_1-20250227	1148.461025	3098.882230	22.969220	61.977645	H22.30.001.CX.16.Plane.03
1315611012_2-20250227	1156.057947	3109.691747	23.121159	62.193835	H22.30.001.CX.16.Plane.03
1315611012_3-20250227	1151.497108	3103.241802	23.029942	62.064836	H22.30.001.CX.16.Plane.03
1315611012_4-20250227	1156.782462	3111.448269	23.135649	62.228965	H22.30.001.CX.16.Plane.03
1315611012_5-20250227	1148.136335	3098.762477	22.962727	61.975250	H22.30.001.CX.16.Plane.03

macaque_cell_slab_plane_coordinates = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Macaque',
    file_name='slab_plane_coordinates',
    dtype={"cell_label": str},
).set_index('cell_label')
macaque_cell_slab_plane_coordinates.head()

slab_plane_coordinates.csv: 100%|██████████| 438M/438M [00:17<00:00, 24.5MMB/s]    

	x_slab_pixel	y_slab_pixel	x_slab_mm	y_slab_mm	registration_subset_label	slab_plane_label
cell_label
1291814091_147-20250227	2031.897028	1110.876577	40.637941	22.217532	None/NA	QM23.50.001.CX.44.Plane.05
1291814091_157-20250227	2031.616827	1111.085137	40.632337	22.221703	None/NA	QM23.50.001.CX.44.Plane.05
1291814091_169-20250227	2029.018004	1110.091872	40.580360	22.201837	None/NA	QM23.50.001.CX.44.Plane.05
1291814091_180-20250227	2029.451516	1111.402151	40.589030	22.228043	None/NA	QM23.50.001.CX.44.Plane.05
1291814091_181-20250227	2029.184113	1111.147536	40.583682	22.222951	None/NA	QM23.50.001.CX.44.Plane.05

marmoset_cell_slab_plane_coordinates = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-Xenium-Marmoset',
    file_name='slab_plane_coordinates',
    dtype={"cell_label": str}
).set_index('cell_label')
marmoset_cell_slab_plane_coordinates.head()

slab_plane_coordinates.csv: 100%|██████████| 131M/131M [00:06<00:00, 20.3MMB/s]   

	x_slab_pixel	y_slab_pixel	x_slab_mm	y_slab_mm	slab_plane_label
cell_label
1382042951_abmkfncc-1-20250227	1117.5521	659.85284	22.351042	13.197058	CJ23.56.004.CX.42.Plane.01
1382042951_abmmdgea-1-20250227	1108.0673	665.55994	22.161345	13.311200	CJ23.56.004.CX.42.Plane.01
1382042951_abmojfce-1-20250227	1115.6930	631.92270	22.313862	12.638454	CJ23.56.004.CX.42.Plane.01
1382042951_abmomnim-1-20250227	1115.3304	631.26380	22.306610	12.625276	CJ23.56.004.CX.42.Plane.01
1382042951_aeiamgga-1-20250227	1098.6119	660.59240	21.972239	13.211849	CJ23.56.004.CX.42.Plane.01

And finally we join each species into their respective cell metadata.

human_cell_metadata = human_cell_metadata.join(
    human_cell_slab_plane_coordinates[['x_slab_mm', 'y_slab_mm']], how='inner'
)
macaque_cell_metadata = macaque_cell_metadata.join(
    macaque_cell_slab_plane_coordinates[['x_slab_mm', 'y_slab_mm']], how='inner'
)
marmoset_cell_metadata = marmoset_cell_metadata.join(
    marmoset_cell_slab_plane_coordinates[['x_slab_mm', 'y_slab_mm']], how='inner'
)

# clean up memory
del human_cell_slab_plane_coordinates
del macaque_cell_slab_plane_coordinates
del marmoset_cell_slab_plane_coordinates

Geojson polygons

The dataset includes polygons outlining select regions of the Basal Ganglia, which can be plotted over slab coordinates. These polygon annotations were manually drawn using a combination of taxonomy cell types, select gene expression from the spatial transcriptomics gene panel, and cell segmentation stains with best attempts made to be consistent across species and planes. Manual annotations include the striatum (STR; encompasses the caudate, putamen, nucleus accumbens, and olfactory tubercle), the globus pallidus (GP) and its nested subregions (the external (GPe) and internal (GPi) segments), the ventral pallidum (VeP), the substantia nigra (SN), and the subthalamic nucleus (STH).. However, please note that these polygons are not derived from registration to the Harmonized Ontology of Mammalian Brain Anatomy (HOMBA). They should instead be used as cellular groupings to orient the user in the anatomical context of the spatial transcriptomics data and subset the data into more manageable, anatomically-related chunks. 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.

First though, we load the term set for these data. While these geojsons are manually drawn and not derived from the HOMBA, we do use the abbreviations and colors associated with these structures through their DHBA identifier when such an identifier exists. These colorings match those of the previous HMBA-BG 10X tutorial. We load these colors and identifiers below.

value_sets = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Macaque',
    file_name='value_sets',
).set_index('label')
pred = value_sets['field'] == 'parcellation_term_symbol'
parcellation_terms = value_sets[pred].copy()
parcellation_terms

value_sets.csv: 100%|██████████| 3.98k/3.98k [00:00<00:00, 38.3kMB/s]

	field	table	description	order	external_identifier	parent_label	color_hex_triplet
label
Basal ganglia and adjacent	parcellation_term_symbol	cell_anatomical_annotations	Basal ganglia and adjacent	1	basal_ganglia_and_adjacent	NaN	#BEBEBE
STR	parcellation_term_symbol	cell_anatomical_annotations	striatum	2	HOMBA:10333	Basal ganglia and adjacent	#cd86cc
GP	parcellation_term_symbol	cell_anatomical_annotations	globus pallidus	3	HOMBA:10342	Basal ganglia and adjacent	#7AE361
GPe	parcellation_term_symbol	cell_anatomical_annotations	external segment of globus pallidus	4	HOMBA:10343	Basal ganglia and adjacent	#53DB33
GPi	parcellation_term_symbol	cell_anatomical_annotations	internal segment of globus pallidus	5	HOMBA:10344	Basal ganglia and adjacent	#A7DB33
VeP	parcellation_term_symbol	cell_anatomical_annotations	ventral pallidus	6	HOMBA:10345	Basal ganglia and adjacent	#9b19f5
STH	parcellation_term_symbol	cell_anatomical_annotations	subthalamic nucleus	7	HOMBA:10466	Basal ganglia and adjacent	#907991
SN	parcellation_term_symbol	cell_anatomical_annotations	substantia nigra	8	HOMBA:12251	Basal ganglia and adjacent	#89522A

We define two functions that will be used in our plotting of the geojsons. The first overplots the geojson over the specified plot, handling the different kinds of geojson shapes/polygons. 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.
        # 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, species: str) -> dict:
    """Load the geojson polygons for the specified slab and species.

    Parameters
    ----------
    slab_name: str
        Name of the slab to load.
    species: str
        Species of the slab, one of 'macaque', 'marmoset', or 'human'.

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

    if species == 'macaque':
        base_dir = 'HMBA-BG-MERSCOPE-Macaque'
    elif species == 'marmoset':
        base_dir = 'HMBA-BG-Xenium-Marmoset'
    elif species == 'human':
        base_dir = 'HMBA-BG-MERSCOPE-Human'
    else:
        raise ValueError('Invalid species')

    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

Manual anatomical annotations

Using the manually drawn anatomical polygons described above, each cell was assigned membership in the annotated anatomical regions. The tables we load below contain boolean columns denoting which anatomical regions the cell could be considered a member of. Note that a given cell can belong to multiple different anatomical annotations as the polygons overlap and can be subsets of each other. For instance, all cells in this release are within the “Basal Ganglia and adjacent” annotation, and most cells in either “GPe” or “GPi” are within the parent annotation “GP”. The final column, parcellation_term_symbol is a single value representing a single unique assignment to a structure based on selecting from STH in reverse order of the columns, keeping the first True assignment. For this notebook, we will use only the parcellation_term_symbol to plot our data.

human_cell_anatomical_annotations = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Human',
    file_name='cell_anatomical_annotations'
).set_index('cell_label')
human_cell_anatomical_annotations.head()

cell_anatomical_annotations.csv: 100%|██████████| 448M/448M [00:16<00:00, 26.5MMB/s]    

	Basal ganglia and adjacent	STR	GP	GPe	GPi	VeP	SN	STH	parcellation_term_symbol
cell_label
1315611012_1-20250227	True	True	False	False	False	False	False	False	STR
1315611012_2-20250227	True	True	False	False	False	False	False	False	STR
1315611012_3-20250227	True	True	False	False	False	False	False	False	STR
1315611012_4-20250227	True	True	False	False	False	False	False	False	STR
1315611012_5-20250227	True	True	False	False	False	False	False	False	STR

macaque_cell_anatomical_annotations = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Macaque',
    file_name='cell_anatomical_annotations'
).set_index('cell_label')
macaque_cell_anatomical_annotations.head()

cell_anatomical_annotations.csv: 100%|██████████| 275M/275M [00:10<00:00, 25.4MMB/s]   

	Basal ganglia and adjacent	STR	GP	GPe	GPi	VeP	SN	STH	parcellation_term_symbol
cell_label
1291814091_147-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1291814091_157-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1291814091_169-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1291814091_180-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1291814091_181-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent

marmoset_cell_anatomical_annotations = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-Xenium-Marmoset',
    file_name='cell_anatomical_annotations'
).set_index('cell_label')
marmoset_cell_anatomical_annotations.head()

cell_anatomical_annotations.csv: 100%|██████████| 116M/116M [00:04<00:00, 24.6MMB/s]   

	Basal ganglia and adjacent	STR	GP	GPe	GPi	VeP	SN	STH	parcellation_term_symbol
cell_label
1382042951_abmkfncc-1-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1382042951_abmmdgea-1-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1382042951_abmojfce-1-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1382042951_abmomnim-1-20250227	True	False	False	False	False	False	False	False	Basal ganglia and adjacent
1382042951_aeiamgga-1-20250227	True	True	False	False	False	False	False	False	STR

Below, we define a function to extract color and order information from the value_set table for each anatomical parcellation term. This will allow us to plot the color of the anatomical region we assigned each cell to using the manually drawn, geojson polygons.

def extract_value_set(cell_metadata_df: pd.DataFrame, input_value_set: pd.DataFrame, input_value_set_label: str):
    """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.
    """
    cell_metadata_df[f'{input_value_set_label}_color'] = input_value_set[
        input_value_set['field'] == input_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'] == input_value_set_label
    ].loc[cell_metadata_df[input_value_set_label]]['order'].values

Merge the manual annotation information into the cell meatadata tables and add in the color and ordering information.

human_cell_metadata = human_cell_metadata.join(
    human_cell_anatomical_annotations[['parcellation_term_symbol']]
)
macaque_cell_metadata = macaque_cell_metadata.join(
    macaque_cell_anatomical_annotations[['parcellation_term_symbol']]
)
marmoset_cell_metadata = marmoset_cell_metadata.join(
    marmoset_cell_anatomical_annotations[['parcellation_term_symbol']]
)
extract_value_set(human_cell_metadata, parcellation_terms, 'parcellation_term_symbol')
extract_value_set(macaque_cell_metadata, parcellation_terms, 'parcellation_term_symbol')
extract_value_set(marmoset_cell_metadata, parcellation_terms, 'parcellation_term_symbol')

human_cell_metadata = human_cell_metadata.sort_values(by='parcellation_term_symbol_order')
macaque_cell_metadata = macaque_cell_metadata.sort_values(by='parcellation_term_symbol_order')
marmoset_cell_metadata = marmoset_cell_metadata.sort_values(by='parcellation_term_symbol_order')

# clean up memory
del human_cell_anatomical_annotations
del macaque_cell_anatomical_annotations
del marmoset_cell_anatomical_annotations

Plotting the Slabs

We define a helper function plot_slabs to visualize the cells for a specified set of slabs either by colorized metadata or gene expression. The code also allows for plotting the annotation polygons over the data if requested. This function can be used to plot individual frames as well.

def plot_slabs(
    df: pd.DataFrame,
    feature: str,
    slab_list: List[str],
    species: 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.
    species: str
        Species to plot (human, macaque, marmoset).
    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_slab_mm' in df.columns:
        x_col = 'x_slab_mm'
    if 'y' in df.columns:
        y_col = 'y'
    elif 'y_slab_mm' in df.columns:
        y_col = 'y_slab_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, species=species), 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_ylim(ylim[1], ylim[0])
        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

Displaying the manual anatomical annotations.

We now have a set of cells in x,y slab coordinates and metadata on the anatomical regions these cells have been manually annotated as members of. Below we plot the cells using the colors as defined in our term set, displaying the cells by manually annotated regions of the Basal Ganglia.

human_slab_list = [
    'H22.30.001.CX.16.Plane.07', 'H22.30.001.CX.17.Plane.07',
    'H22.30.001.CX.19.Plane.07', 'H22.30.001.CX.20.Plane.03', 'H22.30.001.CX.20.Plane.05'
]
color_type = 'parcellation_term_symbol_color'
fig, ax = plot_slabs(
    human_cell_metadata,
    color_type,
    human_slab_list,
    species='human',
    plot_polygons=True,
)
fig.suptitle(f'Human: Manual Anatomical Annotations', size=20)
plt.show()

H22.30.001.CX.16.Plane.07.geojson: 100%|██████████| 5.47k/5.47k [00:00<00:00, 39.6kMB/s]
H22.30.001.CX.17.Plane.07.geojson: 100%|██████████| 8.08k/8.08k [00:00<00:00, 131kMB/s]
H22.30.001.CX.19.Plane.07.geojson: 100%|██████████| 5.97k/5.97k [00:00<00:00, 130kMB/s]
H22.30.001.CX.20.Plane.03.geojson: 100%|██████████| 3.00k/3.00k [00:00<00:00, 23.6kMB/s]
H22.30.001.CX.20.Plane.05.geojson: 100%|██████████| 2.35k/2.35k [00:00<00:00, 25.0kMB/s]

macaque_slab_list = [
    'QM23.50.001.CX.42.Plane.03', 'QM23.50.001.CX.43.Plane.07', 'QM23.50.001.CX.44.Plane.05',
    'QM23.50.001.CX.45.Plane.05', 'QM23.50.001.CX.46.Plane.05'
]
color_type = 'parcellation_term_symbol_color'
fig, ax = plot_slabs(
    macaque_cell_metadata,
    color_type,
    macaque_slab_list,
    species='macaque',
    plot_polygons=True
)
fig.suptitle(f'Macaque: Manual Anatomical Annotations', size=20)
plt.show()

QM23.50.001.CX.42.Plane.03.geojson: 100%|██████████| 2.33k/2.33k [00:00<00:00, 40.5kMB/s]
QM23.50.001.CX.43.Plane.07.geojson: 100%|██████████| 2.99k/2.99k [00:00<00:00, 53.5kMB/s]
QM23.50.001.CX.44.Plane.05.geojson: 100%|██████████| 5.53k/5.53k [00:00<00:00, 91.0kMB/s]
QM23.50.001.CX.45.Plane.05.geojson: 100%|██████████| 6.83k/6.83k [00:00<00:00, 58.3kMB/s]
QM23.50.001.CX.46.Plane.05.geojson: 100%|██████████| 5.00k/5.00k [00:00<00:00, 95.0kMB/s]

marmoset_slab_list = [
    'CJ23.56.004.CX.42.Plane.05', 'CJ23.56.004.CX.42.Plane.22', 'CJ23.56.004.CX.43.Plane.05',
    'CJ23.56.004.CX.43.Plane.15', 'CJ23.56.004.CX.44.Plane.03'
]
color_type = 'parcellation_term_symbol_color'
fig, ax = plot_slabs(
    marmoset_cell_metadata,
    color_type,
    marmoset_slab_list,
    species='marmoset',
    plot_polygons=True
)
fig.suptitle(f'Marmoset: Manual Anatomical Annotations', size=20)
plt.show()

CJ23.56.004.CX.42.Plane.05.geojson: 100%|██████████| 6.43k/6.43k [00:00<00:00, 53.2kMB/s]
CJ23.56.004.CX.42.Plane.22.geojson: 100%|██████████| 7.33k/7.33k [00:00<00:00, 68.0kMB/s]
CJ23.56.004.CX.43.Plane.05.geojson: 100%|██████████| 12.2k/12.2k [00:00<00:00, 308kMB/s]
CJ23.56.004.CX.43.Plane.15.geojson: 100%|██████████| 9.98k/9.98k [00:00<00:00, 127kMB/s]
CJ23.56.004.CX.44.Plane.03.geojson: 100%|██████████| 5.76k/5.76k [00:00<00:00, 85.0kMB/s]

Taxonomy mapping

Now that we have our base functions for plotting set up and we’ve exercised them using the parcellation colors, let’s load the taxonomy information so that we can plot the spatial locations of various cell types. Below, we load and merge all the information we need to plot cell colors and taxonomy as done in the previous HMBA-BG tutorials.

cluster = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-taxonomy',
    file_name='cluster',
    dtype={'number_of_cells': 'Int64'}
).rename(columns={'label': 'cluster_annotation_term_label'}).set_index('cluster_annotation_term_label')
cluster_annotation_term_set = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-taxonomy',
    file_name='cluster_annotation_term_set'
).rename(columns={'label': 'cluster_annotation_term_label'})

cluster_annotation_term = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-taxonomy',
    file_name='cluster_annotation_term',
).rename(columns={'label': 'cluster_annotation_term_label'}).set_index('cluster_annotation_term_label')

cluster_to_cluster_annotation_membership = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-taxonomy',
    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_annotation_term_set = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-taxonomy',
    file_name='cluster_annotation_term_set'
).rename(columns={'label': 'cluster_annotation_term_label'})

# term_sets = abc_cache.get_metadata_dataframe(directory='WHB-taxonomy', file_name='cluster_annotation_term_set').set_index('label')
cluster_details = membership_groupby['cluster_annotation_term_name'].first().unstack()
cluster_order = membership_groupby['term_order'].first().unstack()
cluster_order.sort_values(['Neighborhood', 'Class', 'Subclass', 'Group', 'Cluster'], inplace=True)
cluster_order.rename(
    columns={'Neighborhood': 'Neighborhood_order',
             'Class': 'Class_order',
             'Subclass': 'Subclass_order',
             'Group': 'Group_order',
             'Cluster': 'Cluster_order'},
    inplace=True
)

cluster_colors = membership_groupby['color_hex_triplet'].first().unstack()
cluster_colors = cluster_colors[cluster_annotation_term_set['name']]
cluster_colors.sort_values(
    ['Neighborhood', 'Class', 'Subclass', 'Group', 'Cluster'],
    inplace=True
)

cluster.csv: 100%|██████████| 79.9k/79.9k [00:00<00:00, 878kMB/s]
cluster_annotation_term_set.csv: 100%|██████████| 223/223 [00:00<00:00, 3.22kMB/s]
cluster_annotation_term.csv: 100%|██████████| 207k/207k [00:00<00:00, 1.16MMB/s] 
cluster_to_cluster_annotation_membership.csv: 100%|██████████| 569k/569k [00:00<00:00, 3.94MMB/s] 

Each of the three species has been mapped to the HMBA 10X taxonomy using MapMyCells. Below, we load the mapping of each cell into the taxonomy. The mmc_results contain quality metrics for the full mapping however, for this tutorial we will only utilize the mapping of cell -> Cluster.

Again, we load each species.

human_mmc_results = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Human',
    file_name='mmc_results'
).set_index('cell_label')

macaque_mmc_results = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-MERSCOPE-Macaque',
    file_name='mmc_results',
).set_index('cell_label')

marmoset_mmc_results = abc_cache.get_metadata_dataframe(
    directory='HMBA-BG-Xenium-Marmoset',
    file_name='mmc_results',
).set_index('cell_label')

mmc_results.csv: 100%|██████████| 898M/898M [00:40<00:00, 22.3MMB/s]    
mmc_results.csv: 100%|██████████| 514M/514M [00:25<00:00, 20.3MMB/s]    
mmc_results.csv: 100%|██████████| 215M/215M [00:09<00:00, 21.9MMB/s]    

The spatial data are registered to an extended version of the taxonomy that contains clusters that are considered Adjacent to the Basal Ganglia. These are cell types that are present in the data but are not considered part of the basal ganglia. These clusters are mapped into a Group, Subclass, Class, and Neighborhood in a single path of the taxonomy called Adjacent. The release consensus taxonomy does not explicitly contain this taxon. Any cell in a cluster that would be annotated into Adjacent, will end up being NaN in the joins between the clusters and the taxonomy. We thus fill in these NaNs, assigning any cell that was not mapped into to a group in the taxonomy an Adjacent cell type at all levels with a grey color. We give these Adjacent taxons an ordering at the very end of the set of taxons at each level.”

Although the BG Consensus Taxonomy is intended to be specific to BG cell types, when the data were sampled for the snRNA-Seq reference, some cells from adjacent non-BG regions were included. Because these cells originated from non-BG regions, we did not assign them a cell type identity (i.e. a label at the Group, Subclass, etc. levels). However, because 1) our spatial transcriptomics datasets also contain cells from adjacent non-BG regions, and 2) MapMyCells requires every cell to be assigned a label from the reference, including Adjacent cells in the reference decreases spurious mapping results in BG-adjacent regions. Adjacent cells are included in the release data for completeness, where they are labeled as “Adjacent” at all taxonomy levels. When visualizing cells in BG regions, we recommend filtering them out.

def fill_adjacent_data(input_df: pd.DataFrame) -> pd.DataFrame:
    output_df = input_df.fillna(
        {'cluster_alias': 'Adjacent',
         'cluster_label': 'Adjacent',
         'Class': 'Adjacent',
         'Cluster': 'Adjacent',
         'Group': 'Adjacent',
         'Neighborhood': 'Adjacent',
         'Subclass': 'Adjacent',
         'Neighborhood_color': '#808080',
         'Class_color': '#808080',
         'Subclass_color': '#808080',
         'Group_color': '#808080',
         'Cluster_color': '#808080',
         'Neighborhood_order': cluster_order.max()['Neighborhood_order'] + 1,
         'Class_order': cluster_order.max()['Class_order'] + 1,
         'Subclass_order': cluster_order.max()['Subclass_order'] + 1,
         'Group_order': cluster_order.max()['Group_order'] + 1,
         'Cluster_order': cluster_order.max()['Cluster_order'] + 1,
        }
    )
    return output_df.astype(
        {
            'Neighborhood_order': 'int',
            'Class_order': 'int',
            'Subclass_order': 'int',
            'Group_order': 'int',
            'Cluster_order': 'int',
        }
    )

Now we merge the cell type classification into the cell metadata. Note that we use the cluster_alias only here and load the full taxonomy colors and orders from the taxonomy itself.

human_cell_metadata = human_cell_metadata.join(human_mmc_results[['Cluster']]).rename(columns={'Cluster': 'cluster_alias'})
human_cell_metadata = human_cell_metadata.join(cluster_details, on='cluster_alias')
human_cell_metadata = human_cell_metadata.join(cluster_colors, on='cluster_alias', rsuffix='_color')
human_cell_metadata = human_cell_metadata.join(cluster_order, on='cluster_alias')
human_cell_metadata = fill_adjacent_data(human_cell_metadata)

human_cell_metadata.head()

	brain_section_label	brain_section_barcode	x_experiment	y_experiment	segmentation_job_id	cell_volume	total_counts	total_counts_genes	total_counts_genes_batch_normalized	n_unique_genes	...	Neighborhood_color	Class_color	Subclass_color	Group_color	Cluster_color	Class_order	Cluster_order	Group_order	Neighborhood_order	Subclass_order
cell_label
1326508078_47253-20250227	H22.30.001.CX.17.09.07.02	1326508078	2537.3235	7008.8370	sis_cellpose_v1	153.942654	73	73	66.153050	36	...	#f2ca7d	#6ec0da	#7f207a	#c96e0a	#f36c66	1	74	3	1	2
1326469674_61070-20250227	H22.30.001.CX.17.06.05.02	1326469674	1731.3826	7811.5690	sis_cellpose_v1	494.349680	207	207	223.686463	44	...	#f2ca7d	#8605d4	#def476	#401e66	#6af4cf	3	161	10	1	7
1326469674_61071-20250227	H22.30.001.CX.17.06.05.02	1326469674	1529.1754	7808.3413	sis_cellpose_v1	46.788776	40	40	43.224437	25	...	#f2ca7d	#8605d4	#e3cba4	#922047	#6520f6	3	221	13	1	9
1326469674_61072-20250227	H22.30.001.CX.17.06.05.02	1326469674	1677.0059	7816.4243	sis_cellpose_v1	226.519122	108	108	116.705980	31	...	#f2ca7d	#8605d4	#def476	#401e66	#6af4cf	3	161	10	1	7
1326469674_61073-20250227	H22.30.001.CX.17.06.05.02	1326469674	1569.8965	7821.2876	sis_cellpose_v1	616.443236	265	265	286.361897	53	...	#f2ca7d	#8605d4	#def476	#401e66	#6af4cf	3	161	10	1	7

5 rows × 56 columns

Load the Macaque and Marmoset data in the same way.

macaque_cell_metadata = macaque_cell_metadata.join(macaque_mmc_results[['Cluster']]).rename(columns={'Cluster': 'cluster_alias'})
macaque_cell_metadata = macaque_cell_metadata.join(cluster_details, on='cluster_alias')
macaque_cell_metadata = macaque_cell_metadata.join(cluster_colors, on='cluster_alias', rsuffix='_color')
macaque_cell_metadata = macaque_cell_metadata.join(cluster_order, on='cluster_alias')
macaque_cell_metadata = fill_adjacent_data(macaque_cell_metadata)

macaque_cell_metadata.head()

	brain_section_label	brain_section_barcode	x_experiment	y_experiment	segmentation_job_id	cell_volume	total_counts	total_counts_genes	total_counts_genes_batch_normalized	n_unique_genes	...	Neighborhood_color	Class_color	Subclass_color	Group_color	Cluster_color	Class_order	Cluster_order	Group_order	Neighborhood_order	Subclass_order
cell_label
1291814091_147-20250227	QM23.50.001.CX.44.01.05.01	1291814091	4432.6885	3979.6113	sis_cellpose_v1	756.243171	330	330	636.778305	83	...	#f2ca7d	#8605d4	#def476	#488edc	#a5e9f2	3	149	9	1	7
1291812419_27709-20250227	QM23.50.001.CX.43.01.05.02	1291812419	4813.3880	3526.3945	sis_cellpose_v1	570.255755	129	128	142.601574	8	...	#f2ca7d	#6ec0da	#604f15	#e16c95	#83d8d6	1	26	1	1	1
1291812419_27708-20250227	QM23.50.001.CX.43.01.05.02	1291812419	4757.3520	3525.8088	sis_cellpose_v1	927.675333	432	430	479.052163	113	...	#f2ca7d	#6ec0da	#7f207a	#c96e0a	#1e975e	1	75	3	1	2
1291812419_27707-20250227	QM23.50.001.CX.43.01.05.02	1291812419	5032.3677	3521.7780	sis_cellpose_v1	309.492010	270	270	300.800195	64	...	#f2ca7d	#8605d4	#def476	#401e66	#8ae6d9	3	167	10	1	7
1291812419_27705-20250227	QM23.50.001.CX.43.01.05.02	1291812419	4802.1180	3520.3070	sis_cellpose_v1	710.057491	536	536	597.144091	96	...	#f2ca7d	#6ec0da	#604f15	#e16c95	#83d8d6	1	26	1	1	1

5 rows × 56 columns

marmoset_cell_metadata = marmoset_cell_metadata.join(marmoset_mmc_results[['Cluster']]).rename(columns={'Cluster': 'cluster_alias'})
marmoset_cell_metadata = marmoset_cell_metadata.join(cluster_details, on='cluster_alias')
marmoset_cell_metadata = marmoset_cell_metadata.join(cluster_colors, on='cluster_alias', rsuffix='_color')
marmoset_cell_metadata = marmoset_cell_metadata.join(cluster_order, on='cluster_alias')
marmoset_cell_metadata = fill_adjacent_data(marmoset_cell_metadata)

marmoset_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	...	Neighborhood_color	Class_color	Subclass_color	Group_color	Cluster_color	Class_order	Cluster_order	Group_order	Neighborhood_order	Subclass_order
cell_label
1382042951_abmkfncc-1-20250227	CJ23.56.004.CX.42.01.06	1382042951	5923.4927	6160.6543	10X cell segmentation kit - default	Segmented by boundary stain (ATP1A1+CD45+E-Cad...	123.00563	72.566090	1	354	...	#f2ca7d	#8605d4	#def476	#488edc	#5a9ad6	3	158	9	1	7
1386591529_oglkpmal-1-20250227	CJ23.56.004.CX.43.01.07	1386591529	10066.0640	9964.9790	10X cell segmentation kit - default	Segmented by nucleus expansion of 5.0µm	289.45157	75.546410	1	95	...	#f2ca7d	#aa6dec	#761646	#cbfc1e	#54ca78	2	113	4	1	3
1386591529_oglkphgf-1-20250227	CJ23.56.004.CX.43.01.07	1386591529	10032.3370	9951.0440	10X cell segmentation kit - default	Segmented by nucleus expansion of 5.0µm	250.75267	57.032345	1	190	...	#f2ca7d	#8605d4	#c68ae4	#f5e689	#f222e1	3	219	12	1	8
1386591529_oglkpgff-1-20250227	CJ23.56.004.CX.43.01.07	1386591529	10078.6850	9936.4230	10X cell segmentation kit - default	Segmented by nucleus expansion of 5.0µm	237.34126	48.091408	1	121	...	#f2ca7d	#cba5ec	#f6390c	#361d61	#1bf6d7	4	253	14	1	10
1386591529_oglikmic-1-20250227	CJ23.56.004.CX.43.01.07	1386591529	9677.9650	11089.3920	10X cell segmentation kit - default	Segmented by nucleus expansion of 5.0µm	233.14172	48.859066	1	99	...	#f2ca7d	#aa6dec	#566c21	#195f8d	#80f3ea	2	130	5	1	4

5 rows × 59 columns

Merging in the taxonomy allows one to select cells by the specific taxons they are associated with, as we will see in a function we create below.

Plotting the taxonomy

Now that the taxonomy is loaded, we plot the the cells showing the spatial locations of the cell types in the Basal Ganglia. Below we plot each species in the dataset for each of the 5 levels in the taxonomy selection using the same set of slabs we used for the above to plot the manual anatomical annotations. As we step through each slab plane, we can see hints of the various regions of the basal ganglia, helped by the manual annotation polygons that we plot of the BG regions.

For clarity, we’ll filter out the Adjacent taxonomy classes as well as the Noneuronal cells. We’ll also leave an example for filtering out just the Adjacent data as well.

# Filter out Adjacent
# filtered_human_metadata = human_cell_metadata[~(human_cell_metadata['Neighborhood'] == 'Adjacent')]
# Filter out both adjacent and non-neuronal cells.
filtered_human_metadata = human_cell_metadata[
    np.logical_and(human_cell_metadata['Neighborhood'] != 'Adjacent',
                   human_cell_metadata['Neighborhood'] != 'Nonneuron')
]

for color_type in ['Neighborhood_color', 'Class_color', 'Subclass_color', 'Group_color', 'Cluster_color']:
    fig, ax = plot_slabs(
        filtered_human_metadata,
        color_type,
        human_slab_list,
        species='human',
        plot_polygons=True
    )
    fig.suptitle(f'Human: {color_type.split("_")[0]}')
    plt.show()

# Filter out Adjacent
# filtered_macaque_metadata = macaque_cell_metadata[~(macaque_cell_metadata['Neighborhood'] == 'Adjacent')]
# Filter out both adjacent and non-neuronal cells.
filtered_macaque_metadata = macaque_cell_metadata[
    np.logical_and(macaque_cell_metadata['Neighborhood'] != 'Adjacent',
                   macaque_cell_metadata['Neighborhood'] != 'Nonneuron')
]
for color_type in ['Neighborhood_color', 'Class_color', 'Subclass_color', 'Group_color', 'Cluster_color']:
    fig, ax = plot_slabs(
        filtered_macaque_metadata,
        color_type,
        macaque_slab_list,
        species='macaque',
        plot_polygons=True
    )
    fig.suptitle(f'Macaque: {color_type.split("_")[0]}')
    plt.show()

# Filter out Adjacent
# filtered_marmoset_metadata = marmoset_cell_metadata[~(marmoset_cell_metadata['Neighborhood'] == 'Adjacent')]
# Filter out both adjacent and non-neuronal cells.
filtered_marmoset_metadata = marmoset_cell_metadata[
    np.logical_and(marmoset_cell_metadata['Neighborhood'] != 'Adjacent',
                   marmoset_cell_metadata['Neighborhood'] != 'Nonneuron')
]
for color_type in ['Neighborhood_color', 'Class_color', 'Subclass_color', 'Group_color', 'Cluster_color']:
    fig, ax = plot_slabs(
        filtered_marmoset_metadata,
        color_type,
        marmoset_slab_list,
        species='marmoset',
        plot_polygons=True
    )
    fig.suptitle(f'Marmoset: {color_type.split("_")[0]}')
    plt.show()

Selecting cell types

We define a slightly modified version of our slab plotting function that allows us to plot specific cell types, coloring by cell type level that is nominally lower in the taxonomy. This allows for us to investigate the spatial distribution of cell types more clearly by, for instance, plotting the spatial distributions of Groups within the a Medium Spiny Neuron subclass.

In general, cells associated with individual cell types and taxons can be selected by masking the DataFrame using the desired column and value in that column, as we showed in our examples above. For instance to select the Class, CN LGE GABA, we can create a mask thusly human_cell_metadata['Class']=='CN LGE GABA' and similarly for the other species. Again, this function can be used to plot individual slab planes and selections as well.

def plot_selected_cells(
    df: pd.DataFrame,
    feature: str,
    selection_list: List[Tuple[str, str]],
    slab: str,
    species: str,
    cmap: Optional[plt.Colormap] = None,
    fig_width: float = 40,
    fig_height: float = 8,
    plot_polygons: bool = True,
    point_size: float = 1.0,
    fix_bounds: bool = False
):
    """Plot spatial transcriptomics data selected a cell type and colored by
    a feature of interest.

    Parameters
    ----------
    df: pandas.DataFrame
        cell metadata data frame with x,y locations to plot.
    feature: str
        Expression value or color column to plot.
    selection_list: List[str]
        Set of cell types to select by. Specify as ('taxonomy_level_name', 'cell type')
    slab: str
        Specific slab to plot.
    species: str
        Species to plot (human, macaque, marmoset).
    cmap: matplotlib.colors.Colormap
        Colormap for expression features.
    fid_width: int
        Width of the plotted figure.
    fig_height: iont
        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.
    """
    
    fig, ax = plt.subplots(1, len(selection_list))
    fig.set_size_inches(fig_width, fig_height)

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

    if 'slab_plane_label' in df.columns:
        slab_filtered = df[df['slab_plane_label']==slab]
    elif 'brain_section_label' in df.columns:
        slab_filtered = df[df['brain_section_label']==slab]

    if 'x' in df.columns:
        x_col = 'x'
    elif 'x_slab_mm' in df.columns:
        x_col = 'x_slab_mm'
    if 'y' in df.columns:
        y_col = 'y'
    elif 'y_slab_mm' in df.columns:
        y_col = 'y_slab_mm'

    if fix_bounds:
        x_min = slab_filtered[x_col].min()
        x_max = slab_filtered[x_col].max()
        y_min = slab_filtered[y_col].min()
        y_max = slab_filtered[y_col].max()
    
    for idx, (col, name) in enumerate(selection_list):

        if 'slab_plane_label' in df.columns:
            filtered = slab_filtered[slab_filtered[col]==name]
        elif 'brain_section_label' in df.columns:
            filtered = slab_filtered[slab_filtered[col]==name]
        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, species=species), ax[idx])
            
        ax[idx].axis('equal')
        ylim = ax[idx].get_ylim()
        if fix_bounds:
            ylim = (y_min, y_max)
            ax[idx].set_xlim(x_min, x_max)
        ax[idx].set_ylim(ylim[1], ylim[0])
        ax[idx].set_xticks([])
        ax[idx].set_yticks([])
        
        ax[idx].set_title(name)
        ax[idx].legend(loc=0)
        
    plt.subplots_adjust(wspace=0.01, hspace=0.01)
    fig.suptitle(f'{species}, {slab}')
    return fig, ax

Using the above function, we plot the cells colored by their Group color and select by specific classes and subclasses. We can see commonalities in the distribution across the different species and sections,of the cell types exceptionally clearly.

for slab in human_slab_list:
    fig, ax = plot_selected_cells(
        human_cell_metadata,
        feature='Group_color',
        selection_list=[
            ('Class', 'CN LGE GABA'),
            ('Subclass', 'STR D1 SPN'),
            ('Subclass', 'STR D2 SPN'),
            ('Subclass', 'STR Eccentric SPN'),
            ('Subclass', 'OT Granule GABA')
        ],
        slab=slab,
        species='human',
        fig_width=30,
        fig_height=6,
        plot_polygons=True,
    )
    plt.show()

for slab in macaque_slab_list:
    fig, ax = plot_selected_cells(
        macaque_cell_metadata,
        feature='Group_color',
        selection_list=[
            ('Class', 'CN LGE GABA'),
            ('Subclass', 'STR D1 SPN'),
            ('Subclass', 'STR D2 SPN'),
            ('Subclass', 'STR Eccentric SPN'),
            ('Subclass', 'OT Granule GABA')
        ],
        slab=slab,
        species='macaque',
        fig_width=30,
        fig_height=6,
        plot_polygons=True,
    )
    plt.show()

for slab in marmoset_slab_list:
    fig, ax = plot_selected_cells(
        marmoset_cell_metadata,
        feature='Group_color',
        selection_list=[
            ('Class', 'CN LGE GABA'),
            ('Subclass', 'STR D1 SPN'),
            ('Subclass', 'STR D2 SPN'),
            ('Subclass', 'STR Eccentric SPN'),
            ('Subclass', 'OT Granule GABA')
        ],
        slab=slab,
        species='marmoset',
        fig_width=30,
        fig_height=6,
        plot_polygons=True,
    )
    plt.show()

In the next tutorial, we will explore gene expression patterns across these spatial transcriptomics datasets.

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