IPython magic command to render matplotlib plots.

%matplotlib inline

Whole Human Brain 10x RNA-seq gene expression data (part 2)#

In this notebook we’ll explore some of the gene expression and combine it with the cell metadata we showed in part 1.

You need to be either connected to the internet or connected to a cache that has the WHB data downloaded already and have already downloaded the example data via the getting started notebook to run this notebook.

import pandas as pd
import numpy as np
import anndata
from pathlib import Path
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 requsted data on disk. For metadata, the cache can also directly serve a 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 have downloaded the data in your system.

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

abc_cache.current_manifest
'releases/20240831/manifest.json'

Create the expanded cell metadata as was done previously in the cluster annotation tutorial and part 1 of the 10X tutorial.

# Load the cell metadata.
cell = abc_cache.get_metadata_dataframe(
    directory='WHB-10Xv3',
    file_name='cell_metadata',
    dtype={'cell_label': str}
)
cell.set_index('cell_label', inplace=True)
print("Number of cells = ", len(cell))

# Load the cluster memembership metadata and combine the data with the cell data.
membership = abc_cache.get_metadata_dataframe(
    directory='WHB-taxonomy',
    file_name='cluster_to_cluster_annotation_membership'
)

term_sets = abc_cache.get_metadata_dataframe(directory='WHB-taxonomy', file_name='cluster_annotation_term_set').set_index('label')
cluster_details = membership.groupby(['cluster_alias', 'cluster_annotation_term_set_name'])['cluster_annotation_term_name'].first().unstack()
cluster_details = cluster_details[term_sets['name']] # order columns
cluster_details.fillna('Other', inplace=True)

cluster_details.sort_values(['supercluster', 'cluster', 'subcluster'], inplace=True)
cluster_colors = membership.groupby(['cluster_alias', 'cluster_annotation_term_set_name'])['color_hex_triplet'].first().unstack()
cluster_colors = cluster_colors[term_sets['name']]
cluster_colors.sort_values(['supercluster', 'cluster', 'subcluster'], inplace=True)
cluster_colors

roi = abc_cache.get_metadata_dataframe(directory='WHB-10Xv3', file_name='region_of_interest_structure_map')
roi.set_index('region_of_interest_label', inplace=True)
roi.rename(columns={'color_hex_triplet': 'region_of_interest_color'},
           inplace=True)

del membership
del term_sets

cell_extended = cell.join(cluster_details, on='cluster_alias')
cell_extended = cell_extended.join(cluster_colors, on='cluster_alias', rsuffix='_color')
cell_extended = cell_extended.join(roi[['region_of_interest_color']], on='region_of_interest_label')

del cluster_details
del cluster_colors
del roi

cell_extended.head(5)
Number of cells =  3369219
cell_barcode barcoded_cell_sample_label library_label feature_matrix_label entity brain_section_label library_method donor_label donor_sex dataset_label ... anatomical_division_label subcluster cluster supercluster neurotransmitter subcluster_color cluster_color supercluster_color neurotransmitter_color region_of_interest_color
cell_label
10X386_2:CATGGATTCTCGACGG CATGGATTCTCGACGG 10X386_2 LKTX_210825_01_B01 WHB-10Xv3-Neurons nuclei H19.30.001.CX.51 10Xv3 H19.30.001 M WHB-10Xv3 ... Myelencephalon URL_312_20 URL_312 Upper rhombic lip VGLUT1 #4CB941 #97B8C8 #80BAED #2BDFD1 #5D6CB2
10X383_5:TCTTGCGGTGAATTGA TCTTGCGGTGAATTGA 10X383_5 LKTX_210818_02_E01 WHB-10Xv3-Neurons nuclei H19.30.002.BS.94 10Xv3 H19.30.002 M WHB-10Xv3 ... Myelencephalon URL_312_20 URL_312 Upper rhombic lip VGLUT1 #4CB941 #97B8C8 #80BAED #2BDFD1 #5D6CB2
10X386_2:CTCATCGGTCGAGCAA CTCATCGGTCGAGCAA 10X386_2 LKTX_210825_01_B01 WHB-10Xv3-Neurons nuclei H19.30.001.CX.51 10Xv3 H19.30.001 M WHB-10Xv3 ... Myelencephalon URL_312_17 URL_312 Upper rhombic lip VGLUT1 #C85E40 #97B8C8 #80BAED #2BDFD1 #5D6CB2
10X378_8:TTGGATGAGACAAGCC TTGGATGAGACAAGCC 10X378_8 LKTX_210809_01_H01 WHB-10Xv3-Neurons nuclei H19.30.002.BS.93 10Xv3 H19.30.002 M WHB-10Xv3 ... Pons URL_312_18 URL_312 Upper rhombic lip VGLUT1 #61C1C2 #97B8C8 #80BAED #2BDFD1 #517DBE
10X387_7:TGAACGTAGTATTCCG TGAACGTAGTATTCCG 10X387_7 LKTX_210825_02_G01 WHB-10Xv3-Neurons nuclei H19.30.001.CX.51 10Xv3 H19.30.001 M WHB-10Xv3 ... Myelencephalon URL_312_16 URL_312 Upper rhombic lip VGLUT1 #45328F #97B8C8 #80BAED #2BDFD1 #5D6CB2

5 rows × 24 columns

Single cell transcriptomes#

The ~3 million cell dataset of WHB has been divided into 2 expression matrices: one with Neurons and one with Nonneurons. Each matrix file is formatted as an annadata, h5ad file with minimal metadata.

Below we show some interactions with data from the 10X expression matricies in the WHB dataset. For a deeper dive into how to access specific gene data from the expression matricies, take a look at general_acessing_10x_snRNASeq_tutorial.ipynb. Below we will use precomputed metadata from these matricies to look at gene expression both in relation to different neurotransmitters and locations across the brain.

First, we list the available metadata in the WHB-10Xv3 dataset again. The two files we will be using in this tutorial are the gene metadata and the example_genes_all_cells_expression table.

abc_cache.list_metadata_files('WHB-10Xv3')
['anatomical_division_structure_map',
 'cell_metadata',
 'donor',
 'example_genes_all_cells_expression',
 'gene',
 'region_of_interest_structure_map']

The table below holds metadata for all genes sequenced in the dataset.

gene = abc_cache.get_metadata_dataframe(directory='WHB-10Xv3', file_name='gene')
gene.set_index('gene_identifier', inplace=True)
print("Number of genes = ", len(gene))
gene.head(5)
Number of genes =  59357
gene_symbol biotype name
gene_identifier
ENSG00000000003 TSPAN6 protein_coding tetraspanin 6
ENSG00000000005 TNMD protein_coding tenomodulin
ENSG00000000419 DPM1 protein_coding dolichyl-phosphate mannosyltransferase subunit...
ENSG00000000457 SCYL3 protein_coding SCY1 like pseudokinase 3
ENSG00000000460 C1orf112 protein_coding chromosome 1 open reading frame 112

We’ll skip accessing these data from the expression matricies specifically for now, however, users can learn how to access specific genes from the released expression matricies in thegeneral_acessing_10x_snRNASeq_tutorial notebook.

The precomputed table below contains expressions for the genes SLC17A6, SLC17A7, SLC32A1, PTPRC, PLP1, AQP4, and TTR for all cells across the the WHB dataset. We then join this with our previously created, cell_extended pandas DataFrame from this tutorial.

example_cells_with_genes = abc_cache.get_metadata_dataframe(
    directory='WHB-10Xv3',
    file_name='example_genes_all_cells_expression'
).set_index('cell_label')
example_cells_with_genes.head()
example_genes_all_cells_expression.csv: 100%|█████████████████████████████████████████████████████████████████████████████████| 233M/233M [00:44<00:00, 5.25MMB/s]
PTPRC SLC17A6 SLC32A1 SLC17A7 TTR PLP1 AQP4
cell_label
10X386_2:CATGGATTCTCGACGG 0.0 0.0 0.0 0.000000 13.652457 0.000000 0.0
10X383_5:TCTTGCGGTGAATTGA 0.0 0.0 0.0 0.000000 10.484976 0.000000 0.0
10X386_2:CTCATCGGTCGAGCAA 0.0 0.0 0.0 0.000000 14.006844 0.000000 0.0
10X378_8:TTGGATGAGACAAGCC 0.0 0.0 0.0 9.348821 0.000000 9.348821 0.0
10X387_7:TGAACGTAGTATTCCG 0.0 0.0 0.0 0.000000 12.046448 0.000000 0.0
cell_extended_with_genes = cell_extended.join(example_cells_with_genes)

Example use cases#

In this section, we show a use case with the example genes SLC17A6, SLC17A7, SLC32A1, PTPRC, PLP1, AQP4, and TTR. These genes were selected because they were presented in Siletti et al. 2023 as markers genes for glutamatergic (SLC17A6, SLC17A7) and GABAergic (SLC32A1) neurons, immune cells (PTPRC), oligodendrocytes (PLP1), astrocytes (AQP4), and choroid plexus (TTR). “Marker genes” have much higher expression in the specified cell type or anatomic structure when compared to all other cells, and in many cases are functionally relevant for those cell types. Additional context about these genes, and additional genes of interest can be found in Siletti et al. 2023.

We define a helper functions aggregate_by_metadata to compute the average expression for a given catergory and later plot_umap to plot cells in a UMAP colorized by metadata or expression values similar to what was used in part 1.

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 selected genes in the whole brain#

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

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

Below, we plot the expression of the genes selected in seven genes. We show the genes vs neurotransmitter type/idendity.

agg = aggregate_by_metadata(cell_extended_with_genes, example_cells_with_genes.columns, 'neurotransmitter')
res = plot_heatmap(agg, 8, 8)
../_images/3f5b000d5d50e20ee963c10488dcab1433531afcdbb140bb73765b33f39a1211.png

Below is by dissection region of interest shows that each of these genes are associated with distinct regions.

agg = aggregate_by_metadata(cell_extended_with_genes, example_cells_with_genes.columns, 'anatomical_division_label')
res = plot_heatmap(agg, 8, 8)
../_images/59e0852b744f62ac2ef451c7643d828aac2ec5e7e73834feba36a28bfa46b2e4.png

And finally by supercluster from the taxonomy.

agg = aggregate_by_metadata(cell_extended_with_genes, example_cells_with_genes.columns, 'supercluster')
res = plot_heatmap(agg, 8, 12)
../_images/1826c05b612b176c521cc548560fc6c5939ef4ccf76d23f94c86f9b0af6bec80.png

We can also visualize the relationship these genes and their location in the UMAP. Note again as in part one that Neurons and Non-neurons should be plotted seperately.

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
neurons_subsampled = cell_extended_with_genes[cell_extended_with_genes['feature_matrix_label'] == 'WHB-10Xv3-Neurons'][::10]
non_neurons_subsampled = cell_extended_with_genes[cell_extended_with_genes['feature_matrix_label'] == 'WHB-10Xv3-Nonneurons'][::10]
fig, ax = plot_umap(neurons_subsampled['x'], neurons_subsampled['y'], val=neurons_subsampled['SLC17A6'], cmap=plt.cm.magma_r)
res = ax.set_title("Neuron: Gene SLC17A6")
fig, ax = plot_umap(non_neurons_subsampled['x'], non_neurons_subsampled['y'], val=non_neurons_subsampled['SLC17A6'], cmap=plt.cm.magma_r)
res = ax.set_title("Non-Neuron: Gene SLC17A6")
../_images/f6bc501ed11f228c37ed54d2d5876cc9a10d23fe719b45679f7bcc24424eeed5.png ../_images/cf43476015629f51221fa1ab4eed704b8f4bc19c04147db441fb324ec3822c01.png
fig, ax = plot_umap(neurons_subsampled['x'], neurons_subsampled['y'], val=neurons_subsampled['SLC32A1'], cmap=plt.cm.magma_r)
res = ax.set_title("Neuron: Gene SLC32A1")
fig, ax = plot_umap(non_neurons_subsampled['x'], non_neurons_subsampled['y'], val=non_neurons_subsampled['SLC32A1'], cmap=plt.cm.magma_r)
res = ax.set_title("Non-Neuron: Gene SLC32A1")
../_images/7e1eda55c959c697eefa77bb396fb8f0574db9d96c47370582e2219afb92a838.png ../_images/b56fe810c1f731fc7a7da06a12cd50e3cc25372305a1fd2799d1aed2b9622f76.png
fig, ax = plot_umap(neurons_subsampled['x'], neurons_subsampled['y'], val=neurons_subsampled['PLP1'], cmap=plt.cm.magma_r)
res = ax.set_title("Neuron: Gene PLP1")
fig, ax = plot_umap(non_neurons_subsampled['x'], non_neurons_subsampled['y'], val=non_neurons_subsampled['PLP1'], cmap=plt.cm.magma_r)
res = ax.set_title("Non-Neuron: Gene PLP1")
../_images/0e50ddd3ddce5dab24687600f548b476ce8f474fea6fdfb8790e6cbebd957cb0.png ../_images/6a3c1062be211dad578dc12e0e14ab74b01f48dabebb3315c855abb640371608.png
fig, ax = plot_umap(neurons_subsampled['x'], neurons_subsampled['y'], val=neurons_subsampled['AQP4'], cmap=plt.cm.magma_r)
res = ax.set_title("Neuron: Gene AQP4")
fig, ax = plot_umap(non_neurons_subsampled['x'], non_neurons_subsampled['y'], val=non_neurons_subsampled['AQP4'], cmap=plt.cm.magma_r)
res = ax.set_title("Non-Neuron: Gene AQP4")
../_images/5e7b391124994f42094c32526ea5347377bdbdfcf6f9df8fe86d31648cd57bf3.png ../_images/ca9195e19c76df50e654a287f6ebcf67ab14f9c2f87744e4d624bd1aee104ece.png

(Plots for SLC17A7, PTPRC, and TTR could be generated the same way, but aren’t shown.)