10x RNA-seq gene expression data (part 2b)#

We can continue to explore our examples looking at the expression of canonical neurotransmitter transporter genes and gene Tac2 over the whole brain.

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']['WMB-10X']['metadata']
view_directory = os.path.join( download_base, 
                               manifest['directory_listing']['WMB-10X']['directories']['metadata']['relative_path'], 
                              'views')
cache_views = False
if cache_views :
    os.makedirs( view_directory, exist_ok=True )

Read in the expanded cell metadata table we created in part 1.

rpath = metadata['cell_metadata_with_cluster_annotation']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath )
cell = pd.read_csv(file, keep_default_na=False)
cell.set_index('cell_label',inplace=True)
cell.columns
Index(['cell_barcode', 'barcoded_cell_sample_label', 'library_label',
       'feature_matrix_label', 'entity', 'brain_section_label',
       'library_method', 'region_of_interest_acronym', 'donor_label',
       'donor_genotype', 'donor_sex', 'dataset_label', 'x', 'y',
       'cluster_alias', 'neurotransmitter', 'class', 'subclass', 'supertype',
       'cluster', 'neurotransmitter_color', 'class_color', 'subclass_color',
       'supertype_color', 'cluster_color', 'region_of_interest_order',
       'region_of_interest_color'],
      dtype='object')

Read in the gene expression dataframe we previously cached in the metadata directory from running part 2a tutorial.

rpath = metadata['example_genes_all_cells_expression']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath)
exp = pd.read_csv(file)
exp.set_index('cell_label',inplace=True)

We define a helper functions aggregate_by_metadata to compute the average expression for a given catergory and plot_umap to plot cells in a UMAP colorized by metadata or expression values.

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
def plot_umap( 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.axis('equal')
    ax.set_xlim(-18,27)
    ax.set_ylim(-18,27)
    ax.set_xticks([])
    ax.set_yticks([])
    
    return fig, ax

Expression of canonical neurotransmitter transporter genes in the whole brain#

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 expression of the 9 solute carrier family genes for all the cells in the dataset. We then group the cells by the assigned neurotransmitter class and compute the mean expression for each group. We define a simple helper function to visualized the expression as a heatmap.

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 = 6)
    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
ntgenes = ['Slc17a7','Slc17a6','Slc17a8','Slc32a1','Slc6a5','Slc18a3','Slc6a3','Slc6a4','Slc6a2']
filtered = exp[ntgenes]
joined = cell.join( filtered )

The results are similar that in part 1. Using data from the whole brain, gene Slc17a7 is now most enriched in glutamatergic assigned cells with some expression in cholinergic types. Gene Slc17a6 is most enriched in noradrenergic, then cholinergic types. Genes Slc6a5, Slc18a3, Slc6a3, Slc6a4 and Slc6a2 shows high specificity to glycinergic, cholinergic, dopaminergic, serotonergic and noradrenergic respectively.

agg = aggregate_by_metadata( joined, ntgenes, 'neurotransmitter' )
agg = agg[ntgenes]
res = plot_heatmap(agg,8,3)
../_images/a1b5b87867729747beec651461cad231edd8873e86b61ee819c31b5c3f34d34c.png

We can also visualize the relationship between transport genes and neurotransmitter identities on the UMAP.

subsampled = joined.loc[::10]
fig, ax = plot_umap(subsampled['x'], subsampled['y'],cc=subsampled['neurotransmitter_color'])
res = ax.set_title("Neurotransmitter identity")
../_images/82a3b1da5b3e515aba6e4ca74a574cb21f0f00d0f783e77254375456107e19b9.png
fig, ax = plot_umap(subsampled['x'], subsampled['y'],val=subsampled['Slc17a7'],cmap=plt.cm.magma_r)
res = ax.set_title("Gene Slc17a7")
../_images/cb14b7ede1a68fabaff9ca93f369010db6bf26b5ac658509a990fb180fc0443e.png
fig, ax = plot_umap(subsampled['x'], subsampled['y'],val=subsampled['Slc32a1'],cmap=plt.cm.magma_r)
res = ax.set_title("Gene Slc32a1")
../_images/1271d0ad6ddcbaeed4304ff9f9e532cbf877c270cf6565abdfe49e814b7568ee.png
fig, ax = plot_umap(subsampled['x'], subsampled['y'],val=subsampled['Slc18a3'],cmap=plt.cm.magma_r)
res = ax.set_title("Gene Slc18a3")
../_images/91b29c71ff20629efb6eed15cc7320df722805100f1d901ce285d11f29b9f5cd.png
fig, ax = plot_umap(subsampled['x'], subsampled['y'],val=subsampled['Slc6a5'],cmap=plt.cm.magma_r)
res = ax.set_title("Gene Slc6a5")
../_images/8d2843d8f601419b30bd8fcaa73afbfe16a35af5197fa581a26982249d438ac1.png

Grouping expression by dissection region of interest shows that each of these genes have distinct spatial patterns. Genes Slc17a7 is enriched in the cortex, Slc32a1 in the striatum, Slc17a6 in the thalamus, Slc6a3 in the midbrain, and Slc6a5 in the medulla.

rpath = metadata['region_of_interest_metadata']['files']['csv']['relative_path']
file = os.path.join( download_base, rpath)
roi = pd.read_csv(file)
roi.set_index('acronym',inplace=True)
agg = aggregate_by_metadata( joined, ntgenes, 'region_of_interest_acronym' )
agg = agg.loc[roi.index] # order roi by anatomical order
res = plot_heatmap(agg,8,8)
../_images/6de11fe2d8bcc043c8ef5a171956232f99a1473b7ace9da080ca8e7d7b856732.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 expression values of Tac2 for all cells across the whole brain. As with the thalamus only example, grouping expression by neurotransmitter show that Tac2 gene is highly enriched in cholinergic cell types. With the rest of brain included, we can observe that Tac2 is also enriched in Glut-GABA cell types as well.

exgenes = ['Tac2']
filtered = exp[exgenes]
joined = cell.join( filtered )
subsampled = joined.loc[::10]
agg = aggregate_by_metadata( joined, exgenes, 'neurotransmitter', True )
plot_heatmap(agg,1,3)
<matplotlib.image.AxesImage at 0x7f8e059d7be0>
../_images/d923624be011012316c8f6cccff049b5972e973488cc7c715aa12df530849411.png

Grouping by 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 type and “06 CTX-CGE GABA” GABAergic cells originating from the caudal ganglionic eminence (CGE).

agg = aggregate_by_metadata( joined, exgenes, 'class', True ).head(8)
plot_heatmap(agg,1,3)
<matplotlib.image.AxesImage at 0x7f8e05921c40>
../_images/874da37fe81a6b4303845dccf32bf2a833d6fbed60b70171ee12d9473d361745.png

At the next level, grouping by subclass reveals enrichment is highly anatomically localized cell types such as the medial habenula (MH), bed nuclei of the stria terminalis (BST), spinal nucleus of the trigeminal (SPVC), main olfactory blub (MOB), central amygdalar nucleus (CEA) and arcuate hypothalamic nucleus (ARH).

agg = aggregate_by_metadata( joined, exgenes, 'subclass', True ).head(15)
subclass_list = agg.index[0:10]
plot_heatmap(agg,1,3)
<matplotlib.image.AxesImage at 0x7f8e05aad1f0>
../_images/3d0786917d74361156d54472d1c64ea12bf3d074b47ffb797d97f4f51692de9d.png

Plotting Tac2 expression on the UMAP enables you see whole brain expression in one view showing Tac2 expression is highly restricted to a small number of clusters. Further we can compare the expression profile with the top 10 enriched subclasses.

fig, ax = plot_umap(subsampled['x'], subsampled['y'],val=subsampled['Tac2'],cmap=plt.cm.magma_r)
res = ax.set_title("Gene Tac2")
../_images/1acdb875ba5e154c6cdc59fc4b9bc78fbbbc3e183d5a715cbb69e42e432bfb35.png
pred = [x in subclass_list for x in joined['subclass']]
filtered = joined[pred]

fig, ax = plot_umap(filtered['x'], filtered['y'],cc=filtered['subclass_color'])
res = ax.set_title("Top 10 enriched subclass for Tac2")
../_images/25f6080f832b7a36324f045e09a4eb5dc877ed525dd5d111707a90f78288dc16.png