Tutorial 4: BMMCs profiled with ASAP-seq (ATAC-seq + Protein Epitopes)#

Tutorials 4 and 5 study a human bone marrow dataset by Mimitou et al., 2021 showing the differentiation of hematopoietic stem cells to multiple developmental fates.

Loading and processing data#

The dataset has two modalities: chromatin accessibility (ATAC-seq) and protein epitopes. Raw data is available on GEO. For analyses, we use count matrices filtered by the dataset authors from here. We converted count matrices to the annotated h5ad format, available for download on figshare.

Download annotated data and import necessary packages.

[1]:

import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'

import ocelli as oci
import anndata as ad
import scanpy as sc
import numpy as np
from scipy.sparse import csr_matrix
from matplotlib.colors import LinearSegmentedColormap

! wget --content-disposition https://figshare.com/ndownloader/articles/28303928/versions/1
! unzip -o 28303928.zip
! rm 28303928.zip

--2025-01-29 18:39:10--  https://figshare.com/ndownloader/articles/28303928/versions/1
Resolving figshare.com (figshare.com)... 52.17.87.113, 54.170.28.5, 2a05:d018:1f4:d000:2b52:700c:3d58:74bf, ...
connected. to figshare.com (figshare.com)|52.17.87.113|:443...
HTTP request sent, awaiting response... 200 OK
Length: 37581543 (36M) [application/zip]
Saving to: ‘28303928.zip’

28303928.zip        100%[===================>]  35.84M  16.4MB/s    in 2.2s

2025-01-29 18:39:13 (16.4 MB/s) - ‘28303928.zip’ saved [37581543/37581543]

Archive:  28303928.zip
 extracting: bmmc_asapseq_atac_hscs.h5ad
 extracting: bmmc_asapseq_atac.h5ad
 extracting: bmmc_asapseq_prot_hscs.h5ad
 extracting: bmmc_asapseq_prot.h5ad

Load the chromatin accessibility annotated data and compute its LDA embedding. The ATAC-seq count matrix comprises 10,927 cells and 3,000 genes (achieved by mapping chromatin accessibility fragments to gene regions with Signac). Loaded data includes 15 precomputed LDA topics.

[2]:

atac = oci.read.h5ad('bmmc_asapseq_atac.h5ad')

# compute LDA topics
# oci.pp.lda(atac, n_components=15, max_iter=50, random_state=17, n_jobs=50, verbose=1)

atac

[2]:

AnnData object with n_obs × n_vars = 10927 × 3000
    obs: 'celltype'
    var: 'modality'
    uns: 'X_lda_params'
    obsm: 'X_lda'
    varm: 'X_lda'

Repeat for protein epitopes data consisting of 238 protein tags. Loaded data includes 15 precomputed LDA topics.

[3]:

prot = oci.read.h5ad('bmmc_asapseq_prot.h5ad')

# compute LDA topics
# oci.pp.lda(prot, n_components=15, max_iter=50, random_state=17, n_jobs=50, verbose=1)

prot

[3]:

AnnData object with n_obs × n_vars = 10927 × 238
    obs: 'celltype'
    var: 'modality'
    uns: 'X_lda_params'
    obsm: 'X_lda'
    varm: 'X_lda'

Multimodal Diffusion Maps#

Create a joint AnnData object storing LDA embeddings of both modalities. Annotated data must include processed modality matrices (in bmmc.obsm) and a list of their keys (in bmmc.uns).

[4]:

# create empty AnnData object
bmmc = ad.AnnData(csr_matrix(([],([],[])), shape=(10927, 1)))

# save modality arrays
bmmc.obsm['atac'] = atac.obsm['X_lda']
bmmc.obsm['prot'] = prot.obsm['X_lda']
bmmc.uns['modalities'] = ['atac', 'prot']

# save cell metadata
bmmc.obs.index = list(atac.obs.index)
bmmc.obs['celltype'] = list(atac.obs['celltype'])

bmmc

[4]:

AnnData object with n_obs × n_vars = 10927 × 1
    obs: 'celltype'
    uns: 'modalities'
    obsm: 'atac', 'prot'

Run MDM to find a multimodal representation of the dataset. Follow a three-step process from previous tutorials:

find nearest neighbors across all modalities,
calculate multimodal weights,
and compute the MDM embedding.

[5]:

oci.pp.neighbors(bmmc, n_neighbors=70, n_jobs=50, verbose=True)
oci.tl.modality_weights(bmmc, random_state=17, n_jobs=50, verbose=True)
oci.tl.MDM(bmmc, n_components=10, random_state=17, n_jobs=50, verbose=True)

[atac]  70 nearest neighbors calculated.
[prot]  70 nearest neighbors calculated.

2025-01-29 18:39:19,224 INFO worker.py:1518 -- Started a local Ray instance.

Multimodal weights estimated.

2025-01-29 18:39:25,289 INFO worker.py:1518 -- Started a local Ray instance.

[atac]  Unimodal Markov chain calculated.
[prot]  Unimodal Markov chain calculated.
Multimodal Markov chain calculated.
Eigendecomposition finished.
10 Multimodal Diffusion Maps components calculated.

Visualizing MDM components#

Create a nearest-neighbors graph in the MDM latent space and visualize it with ForceAtlas2.

[6]:

oci.pp.neighbors(bmmc, x=['X_mdm'], n_neighbors=20, n_jobs=50, verbose=True)
oci.tl.neighbors_graph(bmmc, x='X_mdm', n_edges=20, verbose=True)
oci.tl.fa2(bmmc, n_components=2, n_jobs=50, random_state=142, n_iter=20000,
           linlogmode=True, flags='--jitterTolerance 0.1 --scalingRatio 0.01 '\
           '--outboundAttractionDistribution true')

[X_mdm] 20 nearest neighbors calculated.
Nearest neighbors-based graph constructed.

Jan 29, 2025 6:39:38 PM org.netbeans.modules.masterfs.watcher.Watcher getNotifierForPlatform
INFO: Native file watcher is disabled
Jan 29, 2025 6:39:39 PM org.gephi.io.processor.plugin.DefaultProcessor process
INFO: # Nodes loaded: 10,927 (10,927 added)
Jan 29, 2025 6:39:39 PM org.gephi.io.processor.plugin.DefaultProcessor process
INFO: # Edges loaded: 218,540 (143,249 added)

*************************25%
*************************50%
*************************75%
*************************100%
Time = 176.437s

Plot the embedding.

[7]:

cdict = {'CD14+ Monocytes-2': '#1f77b4',
         'CD14+ Monocytes-3': '#ff7f0e',
         'CD8+ naive T cells-1': '#2ca02c',
         'CD8+ naive T cells-2': '#d62728',
         'CD4+ naive T cells': '#9467bd',
         'B cells': '#8c564b',
         'Neutrophil progenitors': '#e377c2',
         'CD14+ Monocytes-1': '#bcbd22',
         'NK cells': '#17becf',
         'Erythroid cells': '#b3c6e5',
         'cDCs': '#f5be82',
         'HSCs': '#a7dd93',
         'CD16+ Monocytes': '#f19d99',
         'Pre-B cells': '#c1b1d2',
         'pDCs': '#be9d96'}

oci.pl.scatter(bmmc, x='X_fa2', c='celltype', cdict=cdict, s=2.5, title='',
               fontsize=10, markerscale=2, figsize=(8, 6))

[7]:

Exploring modality weights#

Multimodal weight distribution shows that chromatin accessibility is more informative for:

hematopoietic stem cells (HSCs),
neutrophil progenitors,
erythroid cells,
plasmacytoid dendritic cells (pDCs).

Proteins are more informative for monocytes, T cells, or B cells.

[8]:

oci.pl.violin(bmmc, groups='celltype', values='weights', cdict=cdict, fontsize=8,
              figsize=(10, 4))

[8]:

[9]:

# define a colormap
cmap = LinearSegmentedColormap.from_list(
    'custom', ['#eaeaea', '#eaeaea', '#efb0b0', '#ff0000', '#ff0000'], N=256)

oci.pl.scatter(bmmc, x='X_fa2', c='weights', cmap=cmap, s=2,
               fontsize=10, figsize=(10, 6))

[9]:

Investigating protein marker levels#

Protein markers are also an excellent tool for annotating cells. In this section, you will plot the normalized expression of selected markers.

Normalize the protein count matrix.

[10]:

# define a colormap
cmap = LinearSegmentedColormap.from_list(
    'custom', ['#eaeaea', '#eaeaea', '#eaeaea', '#EFB0B0', '#FF0000', '#000000'], N=256)

# normalize the count matrix
sc.pp.normalize_total(prot, target_sum=10000)

CD14 (CD14+ Monocytes)#

[11]:

# specify a marker
marker = 'CD14'
# save normalized expression do bmmc.obs
bmmc.obs[marker] = prot[:, marker].X.toarray()
# plot normalized expression
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[11]:

CD71 (Erythroid cells)#

[12]:

marker = 'CD71'
bmmc.obs[marker] = prot[:, marker].X.toarray()
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[12]:

CD34 (Hematopoietic stem cells)#

[13]:

marker = 'CD34'
bmmc.obs[marker] = prot[:, marker].X.toarray()
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[13]:

CD123 (Plasmacytoid dendritic cells)#

[14]:

marker = 'CD123'
bmmc.obs[marker] = prot[:, marker].X.toarray()
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[14]:

CD19 (B cells)#

[15]:

marker = 'CD19'
bmmc.obs[marker] = prot[:, marker].X.toarray()
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[15]:

CD4 (CD4+ naive T cells)#

[16]:

marker = 'CD4-1'
bmmc.obs[marker] = prot[:, marker].X.toarray()
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[16]:

CD8 (CD8+ naive T cells)#

[17]:

marker = 'CD8a'
bmmc.obs[marker] = prot[:, marker].X.toarray()
oci.pl.scatter(bmmc, x='X_fa2', c=marker, cmap=cmap, s=2,
               fontsize=10, vmax=np.percentile(bmmc.obs[marker], 98),
               figsize=(8, 5))

[17]:

In Tutorial 5, you will use Ocelli’s multimodal diffusion functionality to explore developmental subtrajectories originating in hematopoietic stem cells.

Tutorial 4: BMMCs profiled with ASAP-seq (ATAC-seq + Protein Epitopes)

Contents