MERFISH whole brain spatial transcriptomics (part 1)#

The MERFISH data consists of 59 coronal sections (10 micron thick) covering the whole anterior to posterior extent of the brain (~200 micron apart) from a single adult male mouse. Each section is assayed using the MERSCOPE platform with a 500 gene panel designed to enable us to map cells to the whole mouse brain taxonomy.

The MERSCOPE platform uses combinatorial labeling, sequential imaging and error robust barcodes to obtain quantitative (counts) RNA-expression of 100s-1000s of genes. Each probe binds to a specific transcripts resulting in “dots” in the image. Cells are segmented using CellPose software based on DAPI/PolyT stains.

Gene expression data is stored either as individual sections or all sections in a single file. The puprose of this set of notebooks is to provide and overview of the data, the file organization and how to combine the data and metadata through example use cases.

You need to be connected to the internet to run this notebook and that you have downloaded the example data via the getting started notebook.

import os
import pandas as pd
import numpy as np
import anndata
import time
import json
import matplotlib.pyplot as plt
import requests

The prerequisite for running this notebook is that the data have been downloaded to local directory maintaining the organization from the manifest.json. Change the download_base variable to where you have downloaded the data in your system.

version = '20231215'
download_base = '../../abc_download_root'

use_local_cache = False
manifest_path = 'releases/%s/manifest.json' % version

if not use_local_cache :
    url = 'https://allen-brain-cell-atlas.s3-us-west-2.amazonaws.com/' + manifest_path
    manifest = json.loads(requests.get(url).text)
else :
    file = os.path.join(download_base,manifest_path)
    with open(file,'rb') as f:
        manifest = json.load(f)
    
metadata = manifest['file_listing']['MERFISH-C57BL6J-638850']['metadata']
view_directory = os.path.join( download_base, 
                               manifest['directory_listing']['MERFISH-C57BL6J-638850']['directories']['metadata']['relative_path'], 
                              'views')
cache_views = False
if cache_views :
    os.makedirs( view_directory, exist_ok=True )

Data overview#

Cell metadata#

Essential cell metadata is stored as a dataframe. Each row represents one cell indexed by a cell label.

Each cell is associated with a brain section label, donor label, donor genotype and donor sex. Each cell also has a set of x, y, z coordinates generated by rotating each section so that it is upright with midline approximatedly in the middle of the frame. The z coordinate takes into account the phyiscal spacing between each section.

Each cell is mapped to the whole mouse brain taxonomy using a hierachical correlation method, resulting in the assignment of a *cluster alias” and average correlation score.

rpath = metadata['cell_metadata']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath)
cell = pd.read_csv(file, dtype={"cell_label":str})
cell.set_index('cell_label',inplace=True)
print("Number of cells = ", len(cell))
cell.head(5)
Number of cells =  3938808
brain_section_label cluster_alias average_correlation_score feature_matrix_label donor_label donor_genotype donor_sex x y z
cell_label
1019171907102340387-1 C57BL6J-638850.37 1408 0.596276 C57BL6J-638850 C57BL6J-638850 wt/wt M 7.226245 4.148963 6.6
1104095349101460194-1 C57BL6J-638850.26 4218 0.641180 C57BL6J-638850 C57BL6J-638850 wt/wt M 5.064889 7.309543 4.2
1017092617101450577 C57BL6J-638850.25 4218 0.763531 C57BL6J-638850 C57BL6J-638850 wt/wt M 5.792921 8.189973 4.0
1018093344101130233 C57BL6J-638850.13 4218 0.558073 C57BL6J-638850 C57BL6J-638850 wt/wt M 3.195950 5.868655 2.4
1019171912201610094 C57BL6J-638850.27 4218 0.591009 C57BL6J-638850 C57BL6J-638850 wt/wt M 5.635732 7.995842 4.4

We define a small helper fuction that uses pandas groupby function to see how many unique items are associated for each field and list them out if the number of items is small.

def print_column_info(df) :
    
    for c in df.columns :
        grouped = df[[c]].groupby(c).count()
        members = ''
        if len(grouped) < 30 :
            members = str(list(grouped.index))
        print("Number of unique %s = %d %s" % (c,len(grouped),members))
print_column_info(cell)
Number of unique brain_section_label = 59 
Number of unique cluster_alias = 5320 
Number of unique average_correlation_score = 3932069 
Number of unique feature_matrix_label = 1 ['C57BL6J-638850']
Number of unique donor_label = 1 ['C57BL6J-638850']
Number of unique donor_genotype = 1 ['wt/wt']
Number of unique donor_sex = 1 ['M']
Number of unique x = 3938808 
Number of unique y = 3938808 
Number of unique z = 59 

We can also bring in the pivot table from the “cluster annotation tutorial” to associate each cell with terms at each cell type classification level and the corresponding color.

taxonomy_metadata = manifest['file_listing']['WMB-taxonomy']['metadata']
rpath = taxonomy_metadata['cluster_to_cluster_annotation_membership_pivoted']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath)
cluster_details = pd.read_csv(file,keep_default_na=False)
cluster_details.set_index('cluster_alias', inplace=True)
cluster_details.head(5)
neurotransmitter class subclass supertype cluster
cluster_alias
1 Glut 01 IT-ET Glut 018 L2 IT PPP-APr Glut 0082 L2 IT PPP-APr Glut_3 0326 L2 IT PPP-APr Glut_3
2 Glut 01 IT-ET Glut 018 L2 IT PPP-APr Glut 0082 L2 IT PPP-APr Glut_3 0327 L2 IT PPP-APr Glut_3
3 Glut 01 IT-ET Glut 018 L2 IT PPP-APr Glut 0081 L2 IT PPP-APr Glut_2 0322 L2 IT PPP-APr Glut_2
4 Glut 01 IT-ET Glut 018 L2 IT PPP-APr Glut 0081 L2 IT PPP-APr Glut_2 0323 L2 IT PPP-APr Glut_2
5 Glut 01 IT-ET Glut 018 L2 IT PPP-APr Glut 0081 L2 IT PPP-APr Glut_2 0325 L2 IT PPP-APr Glut_2
rpath = taxonomy_metadata['cluster_to_cluster_annotation_membership_color']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath)
cluster_colors = pd.read_csv(file)
cluster_colors.set_index('cluster_alias', inplace=True)
cluster_colors.head(5)
neurotransmitter_color class_color subclass_color supertype_color cluster_color
cluster_alias
1 #2B93DF #FA0087 #0F6632 #266DFF #64661F
2 #2B93DF #FA0087 #0F6632 #266DFF #CCA73D
3 #2B93DF #FA0087 #0F6632 #002BCC #99000D
4 #2B93DF #FA0087 #0F6632 #002BCC #5C8899
5 #2B93DF #FA0087 #0F6632 #002BCC #473D66
cell_extended = cell.join(cluster_details,on='cluster_alias')
cell_extended = cell_extended.join(cluster_colors,on='cluster_alias')
cell_extended.head(5)
brain_section_label cluster_alias average_correlation_score feature_matrix_label donor_label donor_genotype donor_sex x y z neurotransmitter class subclass supertype cluster neurotransmitter_color class_color subclass_color supertype_color cluster_color
cell_label
1019171907102340387-1 C57BL6J-638850.37 1408 0.596276 C57BL6J-638850 C57BL6J-638850 wt/wt M 7.226245 4.148963 6.6 04 DG-IMN Glut 038 DG-PIR Ex IMN 0141 DG-PIR Ex IMN_2 0515 DG-PIR Ex IMN_2 #666666 #16f2f2 #3D53CC #CC7A3D #73FFBF
1104095349101460194-1 C57BL6J-638850.26 4218 0.641180 C57BL6J-638850 C57BL6J-638850 wt/wt M 5.064889 7.309543 4.2 Glut 23 P Glut 235 PG-TRN-LRN Fat2 Glut 0953 PG-TRN-LRN Fat2 Glut_1 4199 PG-TRN-LRN Fat2 Glut_1 #2B93DF #6b5ca5 #9B7ACC #990041 #663D63
1017092617101450577 C57BL6J-638850.25 4218 0.763531 C57BL6J-638850 C57BL6J-638850 wt/wt M 5.792921 8.189973 4.0 Glut 23 P Glut 235 PG-TRN-LRN Fat2 Glut 0953 PG-TRN-LRN Fat2 Glut_1 4199 PG-TRN-LRN Fat2 Glut_1 #2B93DF #6b5ca5 #9B7ACC #990041 #663D63
1018093344101130233 C57BL6J-638850.13 4218 0.558073 C57BL6J-638850 C57BL6J-638850 wt/wt M 3.195950 5.868655 2.4 Glut 23 P Glut 235 PG-TRN-LRN Fat2 Glut 0953 PG-TRN-LRN Fat2 Glut_1 4199 PG-TRN-LRN Fat2 Glut_1 #2B93DF #6b5ca5 #9B7ACC #990041 #663D63
1019171912201610094 C57BL6J-638850.27 4218 0.591009 C57BL6J-638850 C57BL6J-638850 wt/wt M 5.635732 7.995842 4.4 Glut 23 P Glut 235 PG-TRN-LRN Fat2 Glut 0953 PG-TRN-LRN Fat2 Glut_1 4199 PG-TRN-LRN Fat2 Glut_1 #2B93DF #6b5ca5 #9B7ACC #990041 #663D63
print_column_info(cell_extended)
Number of unique brain_section_label = 59 
Number of unique cluster_alias = 5320 
Number of unique average_correlation_score = 3932069 
Number of unique feature_matrix_label = 1 ['C57BL6J-638850']
Number of unique donor_label = 1 ['C57BL6J-638850']
Number of unique donor_genotype = 1 ['wt/wt']
Number of unique donor_sex = 1 ['M']
Number of unique x = 3938808 
Number of unique y = 3938808 
Number of unique z = 59 
Number of unique neurotransmitter = 10 ['', 'Chol', 'Dopa', 'GABA', 'GABA-Glyc', 'Glut', 'Glut-GABA', 'Hist', 'Nora', 'Sero']
Number of unique class = 34 
Number of unique subclass = 338 
Number of unique supertype = 1201 
Number of unique cluster = 5320 
Number of unique neurotransmitter_color = 10 ['#03EDFF', '#0a9964', '#2B93DF', '#533691', '#666666', '#73E785', '#820e57', '#FF3358', '#fcf04b', '#ff7621']
Number of unique class_color = 34 
Number of unique subclass_color = 338 
Number of unique supertype_color = 1201 
Number of unique cluster_color = 5296 

For convenience, we can cache this view for later reuse.

if cache_views :
    file = os.path.join( view_directory, 'cell_metadata_with_cluster_annotation.csv')
    cell_extended.to_csv( file )

Gene panel#

A 500 gene panel was selected to enable faciliate the mapping to the whole mouse brain taxonomy. Each gene is uniquely identifier by an Ensembl ID. It is best practice to gene identifier to for tracking and data interchange as gene symbols are not unique and can change over time. The expression matrices also includes data for 50 “blanks” codewords in the spatial transcriptomic codebook that are not associated with any gene. They are included to measure the false discovery rate of gene detection.

Each row of the gene dataframe has Ensembl gene identifier, the identifier of the transcript the probe targets, a gene symbol and name.

rpath = metadata['gene']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath)
gene = pd.read_csv(file)
gene.set_index('gene_identifier',inplace=True)
print("Number of genes = ", len(gene))
gene.head(5)
Number of genes =  550
gene_symbol transcript_identifier name mapped_ncbi_identifier
gene_identifier
ENSMUSG00000026778 Prkcq ENSMUST00000028118 protein kinase C, theta NCBIGene:18761
ENSMUSG00000026837 Col5a1 ENSMUST00000028280 collagen, type V, alpha 1 NCBIGene:12831
ENSMUSG00000001985 Grik3 ENSMUST00000030676 glutamate receptor, ionotropic, kainate 3 NCBIGene:14807
ENSMUSG00000039323 Igfbp2 ENSMUST00000047328 insulin-like growth factor binding protein 2 NCBIGene:16008
ENSMUSG00000048387 Osr1 ENSMUST00000057021 odd-skipped related transcription factor 1 NCBIGene:23967

Gene expression matrix#

Expression values for all sections are stored in an anndata h5ad format with minimal metadata. There are two h5ad files one storing the raw counts and the other log normalization of it. In this next section, we provide example code on how to open the file and connect with the rich cell level metadata discussed above.

expression_matrices = manifest['file_listing']['MERFISH-C57BL6J-638850']['expression_matrices']
expression_matrices['C57BL6J-638850']
{'log2': {'files': {'h5ad': {'version': '20230830',
    'relative_path': 'expression_matrices/MERFISH-C57BL6J-638850/20230830/C57BL6J-638850-log2.h5ad',
    'url': 'https://allen-brain-cell-atlas.s3.us-west-2.amazonaws.com/expression_matrices/MERFISH-C57BL6J-638850/20230830/C57BL6J-638850-log2.h5ad',
    'size': 7627589574}}},
 'raw': {'files': {'h5ad': {'version': '20230830',
    'relative_path': 'expression_matrices/MERFISH-C57BL6J-638850/20230830/C57BL6J-638850-raw.h5ad',
    'url': 'https://allen-brain-cell-atlas.s3.us-west-2.amazonaws.com/expression_matrices/MERFISH-C57BL6J-638850/20230830/C57BL6J-638850-raw.h5ad',
    'size': 7627589574}}}}
rpath = expression_matrices['C57BL6J-638850']['log2']['files']['h5ad']['relative_path']
file = os.path.join( download_base, rpath)
print(file)
../../abc_download_root/expression_matrices/MERFISH-C57BL6J-638850/20230830/C57BL6J-638850-log2.h5ad

We use the anndata’s read_h5ad function to open the package for the log2 normalization file. The “backed=’r’” makes use of the lazy loading functionality to only load required data. By default, anndata will load the entire expression matrix in memory.

adata = anndata.read_h5ad(file,backed='r')
print(adata)
AnnData object with n_obs × n_vars = 4334174 × 550 backed at '../../abc_download_root/expression_matrices/MERFISH-C57BL6J-638850/20230830/C57BL6J-638850-log2.h5ad'
    obs: 'brain_section_label'
    var: 'gene_symbol', 'transcript_identifier'
    uns: 'accessed_on', 'src'

Genes are represented as “variables”. For this data, the var dataframe is indexed by the Ensembl gene identifier with one metadata column gene symbol

print("Number of genes = ", len(adata.var))
adata.var.index[0:5]
Number of genes =  550
Index(['ENSMUSG00000026778', 'ENSMUSG00000026837', 'ENSMUSG00000001985',
       'ENSMUSG00000039323', 'ENSMUSG00000048387'],
      dtype='object', name='gene_identifier')

We define a small helper function plot section to visualize the cells in anatomical context colorized by: neurotransmitter identity, cell types division, class and subclass.

def plot_section( xx, yy, cc=None, val=None, fig_width = 8, fig_height = 8, cmap=None ) :
    fig, ax = plt.subplots()
    fig.set_size_inches(fig_width, fig_height)
    if cmap is not None :
        plt.scatter(xx,yy,s=0.5,c=val,marker='.',cmap=cmap)
    elif cc is not None :
        plt.scatter(xx,yy,s=0.5,color=cc,marker='.')
    ax.set_ylim(11,0)
    ax.set_xlim(0,11)
    ax.axis('equal')
    ax.set_xticks([])
    ax.set_yticks([])
    return fig, ax

Let pick one section “C57BL6J-638850.38” for visualization.

pred = (cell_extended['brain_section_label'] == 'C57BL6J-638850.38')
section = cell_extended[pred]
print("Number of cells in section:",len(section))
Number of cells in section: 120186
fig, ax = plot_section(section['x'], section['y'], section['neurotransmitter_color'])
res = ax.set_title("Neuortransmitter Identity")
../_images/445cba94da494ee7d1349ee665a1b89af52cd301397a9f2f32f52b7c453b83e5.png
fig, ax = plot_section(section['x'], section['y'],section['class_color'])
res = ax.set_title("Cell Types: Class")
../_images/d998d76519887c1f66dc12cd589ea5d046a04197edc1f71446003ed6d3545205.png
fig, ax = plot_section(section['x'], section['y'],section['subclass_color'])
res = ax.set_title("Cell Types: Subclass")
../_images/4d6ca05485f6a66f2a3fb94283c785cf2431c9b671ac87df309f3c111384c68b.png

Example use cases#

In this section, we explore the same two use cases as with the 10x data. The first example looks at the expression of nine canonical neurotransmitter transporter genes and the second the expression of gene Tac2. To support these use cases, we will create a smaller submatrix (all cells and 10 genes) that read it into memory. Note this operation takes around 10 seconds.

ntgenes = ['Slc17a7','Slc17a6','Slc17a8','Slc32a1','Slc6a5','Slc6a3','Slc6a4']
exgenes = ['Tac2']
gnames = ntgenes + exgenes
pred = [x in gnames for x in adata.var.gene_symbol]
gene_filtered = adata.var[pred]
gene_filtered
gene_symbol transcript_identifier
gene_identifier
ENSMUSG00000030500 Slc17a6 ENSMUST00000032710
ENSMUSG00000037771 Slc32a1 ENSMUST00000045738
ENSMUSG00000025400 Tac2 ENSMUST00000026466
ENSMUSG00000039728 Slc6a5 ENSMUST00000056442
ENSMUSG00000070570 Slc17a7 ENSMUST00000085374
ENSMUSG00000019935 Slc17a8 ENSMUST00000020102
ENSMUSG00000021609 Slc6a3 ENSMUST00000022100
ENSMUSG00000020838 Slc6a4 ENSMUST00000021195
start = time.process_time()
asubset = adata[:,gene_filtered.index].to_memory()
print("time taken: ", time.process_time() - start)
print(asubset)
time taken:  6.811778808999996
AnnData object with n_obs × n_vars = 4334174 × 8
    obs: 'brain_section_label'
    var: 'gene_symbol', 'transcript_identifier'
    uns: 'accessed_on', 'src'

We define two helper functions to (1) create_expression_dataframe: create joined gene expression and cell metadata dataframe based for a set of input genes and (2) aggregate_by_metadata which computes the average expression for each term in a given category.

def create_expression_dataframe( ad, gf ) :
    gdata = ad[:,gf.index].to_df()
    gdata.columns = gf.gene_symbol
    joined = section.join( gdata )
    return joined
def aggregate_by_metadata( df, gnames, value, sort = False ) :
    grouped = df.groupby(value)[gnames].mean()
    if sort :
        grouped = grouped.sort_values(by=gnames[0],ascending=False)
    return grouped

Expression of canonical neurotransmitter transporter genes#

During analysis, clusters were assigned neurotransmitter identities based on the expression of of canonical neurotransmitter transporter genes. In this example, we create a dataframe comprising of cells in the ‘C57BL6J-638850.40’ section and expression of the 9 solute carrier family genes. We then group the cells by the assigned neurotransmitter class and compute the mean expression for each group and visualized as a colorized table. The results are similar of the 10x data restricted to cells from the thalamus.

def plot_heatmap( df, fig_width = 8, fig_height = 4, cmap=plt.cm.magma_r) :

    arr = df.to_numpy()

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

    im = ax.imshow(arr, cmap = cmap, aspect='auto', vmin=0, vmax=5)
    xlabs = df.columns.values
    ylabs = df.index.values

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

    ax.set_yticks(range(len(ylabs)))
    res = ax.set_yticklabels( ylabs )
    
    return im
pred = [x in ntgenes for x in asubset.var.gene_symbol]
gf = asubset.var[pred]

ntexp = create_expression_dataframe( asubset, gf )
agg = aggregate_by_metadata(ntexp, gf.gene_symbol, 'neurotransmitter')
agg = agg[ntgenes]
im = plot_heatmap(agg,8,3)
../_images/79812ca7489a385321dccc6d367be939a641148b9c27e8123cb6b29bde421fee.png

The expression data can also be visualized as scatter point in anatomical context. Gene Slc17a7 is highly enriched in the cortex , hippocampus, medial habenula with some expression in the thalamus. On the other hand gene Slc32a1 is enriched in the striatum with sparse expression on the cortical plate.

fig, ax = plot_section(ntexp['x'], ntexp['y'],val=ntexp['Slc17a7'],cmap=plt.cm.magma_r)
res = ax.set_title('Gene Slc17a7')
../_images/37f365aab18e95ef4de1ff05aa617b83b3d9257deef4b80fab327e879d9ab386.png
fig, ax = plot_section(ntexp['x'], ntexp['y'],val=ntexp['Slc32a1'],cmap=plt.cm.magma_r)
res = ax.set_title('Gene Slc32a1')
../_images/0553d9589ed3c8c42e6c689c5faf94e8b817f827b4ffd86e6b0d3683c38248ee.png

Expression of Tachykinin 2 (Tac2) in the whole brain#

In mice, the tachykinin 2 (Tac2) gene encodes neuropeptide called neurokinin B (NkB). Tac2 is produced by neurons in specific regions of the brain know to be invovled in emotion and social behavior. Based on ISH data from the Allen Mouse Brain Atlas, Tac 2 is sparsely expressed in the mouse isocortex and densely enriched is specific subcortical regions such the medial habenula (MH), the amygdala and hypothalamus.

In this example, we create a dataframe comprising of cells in section ‘C57BL6J-638850.40’ and expression values of Tac2 in those cells.

gf = asubset.var[asubset.var.gene_symbol == 'Tac2']
tac2 = create_expression_dataframe( asubset, gf )

Grouping cells by neurotransmitter identites and computing the mean expression in each group, we can observed that Tac2 gene is highly enriched in cholinergic cell types with some expression in Glut-GABA types.

agg = aggregate_by_metadata(tac2, gf.gene_symbol, 'neurotransmitter', True).head(8)
im = plot_heatmap(agg,1,3)
../_images/08cdfbf8158e7a3fafd71995c5c27959ad9c4ffa0f81de090626839f1e8680c8.png

Grouping cells by cell types class, shows that Tac2 is enriched in class “16 MH-LH Glut” with cells restricted to the medial (MH) and lateral (LH) habenula and a mixture of glutamatergic and cholinergic types.

agg = aggregate_by_metadata(tac2, gf.gene_symbol, 'class', True).head(8)
im = plot_heatmap(agg,1,3)
../_images/8e88ecc4fa7d9f636126c0c247aaf3e61f92e46ed75e5c8d67c16b67f871e217.png

At the next level, grouping by subclass, shows Tac2 is enriched in sublcasses with specificity to the medial habenula (MH), bed muclei of the stria terminalis (BST) and main olfactory blub (MOB).

agg = aggregate_by_metadata(tac2, gf.gene_symbol, 'subclass', True).head(8)
im = plot_heatmap(agg,1,3)
../_images/4aac2ae3abaecffec190cd8838dd267edf47ec36f3e658f0ada33fc14df48617.png

Plotting Tac2 expression in anatomical coordinates results in a image that is comparable to the ISH image data showing enriched expression in the medial habenula, amydala and hypothalamus and sparse expression in the cortex.

fig,ax = plot_section(tac2['x'], tac2['y'],val=tac2['Tac2'],cmap=plt.cm.magma_r)
res = ax.set_title('Gene Tac2')
../_images/e9ddca4d768279de36e1eb1dab1145adf9ced0fade25fbba4e05b74432e238fd.png

Close h5ad file and clean up

adata.file.close()
del adata