Tutorial 3: Diffusion-based multimodal imputation#
Tutorials 2 and 3 investigate a multimodal single-cell hair follicle dataset by Ma et al., 2020. Data comes from female mouse dorsal skin and consists of two single-cell modalities: ATAC-seq (chromatin accessibility, raw data available on GEO) and RNA-seq (gene expression, raw data available on GEO).
The single-cell profiling of cells inevitably generates data of high sparsity levels, resulting in 0-inflated feature values that obstruct the inference of relationships between features. This section shows how to use Ocelli to strengthen a weak single-cell signal using diffusion-based multimodal imputation. This process restores the missing expression values by exponentiation of the multimodal Markov matrix to a power \(t\) and multiplying the result by the original expression matrix.
Below, you will impute values of RNA-seq gene expression and ATAC-seq promoter activity of two hair follicle’s Henle’s lineage markers, Krt71 and Krt73.
Loading and processing data#
Download filtered data, available on figshare, and import necessary packages.
[1]:
import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
import ocelli as oci
import anndata as ad
import scvelo as scv
from scipy.sparse import csr_matrix
from matplotlib.colors import LinearSegmentedColormap
! wget --content-disposition https://figshare.com/ndownloader/articles/28303598/versions/1
! unzip -o 28303598.zip
! rm 28303598.zip
--2025-01-29 18:37:18-- https://figshare.com/ndownloader/articles/28303598/versions/1
52.17.87.113, 54.170.28.5, 2a05:d018:1f4:d003:6dec:dfbb:5bc1:fcf0, ...
Connecting to figshare.com (figshare.com)|52.17.87.113|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 167301182 (160M) [application/zip]
Saving to: ‘28303598.zip’
28303598.zip 100%[===================>] 159.55M 16.7MB/s in 11s
2025-01-29 18:37:29 (14.8 MB/s) - ‘28303598.zip’ saved [167301182/167301182]
Archive: 28303598.zip
extracting: hf_shareseq_atac_peaks.h5ad
extracting: hf_shareseq_atac.h5ad
extracting: hf_shareseq_rna.h5ad
extracting: hf_shareseq_signature_HS_Co.csv
extracting: hf_shareseq_signature_HS_Me.csv
extracting: hf_shareseq_signature_IRS_He.csv
extracting: hf_shareseq_signature_IRS_Hu.csv
Preprocess data as in tutorial 2.
[2]:
# prepare RNA-seq data
rna = oci.read.h5ad('hf_shareseq_rna.h5ad')
scv.pp.normalize_per_cell(rna, counts_per_cell_after=10000)
scv.pp.filter_genes_dispersion(rna, n_top_genes=1000)
scv.pp.log1p(rna)
scv.tl.velocity(rna, mode='stochastic')
scv.tl.velocity_graph(rna)
# prepare ATAC-seq data
atac = oci.read.h5ad('hf_shareseq_atac.h5ad')
# create a joint AnnData object
hf = ad.AnnData(csr_matrix(([],([],[])), shape=(7160, 1000)))
hf.obsm['rna'] = rna.obsm['X_lda']
hf.obsm['atac'] = atac.obsm['X_lda']
hf.uns['modalities'] = ['rna', 'atac']
hf.obs.index = list(rna.obs.index)
hf.obs['celltype'] = list(rna.obs['celltype'])
hf.uns['velocity_graph'] = rna.uns['velocity_graph']
hf.uns['velocity_graph_neg'] = rna.uns['velocity_graph_neg']
hf.layers['velocity'] = rna.layers['velocity']
Normalized count data: X, spliced, unspliced.
Extracted 1000 highly variable genes.
computing neighbors
/tmp/ipykernel_516574/913483603.py:6: DeprecationWarning: `log1p` is deprecated since scVelo v0.3.0 and will be removed in a future version. Please use `log1p` from `scanpy.pp` instead.
scv.pp.log1p(rna)
finished (0:00:03) --> added
'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
finished (0:00:00) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
computing velocities
finished (0:00:00) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 1/256 cores)
finished (0:00:03) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
Multimodal Diffusion Maps#
Repeat the MDM training procedure from tutorial 2. For imputation, Ocelli needs eigenvectors and eigenvalues of the multimodal Markov matrix. To save them in the hf object, set save_eigvec=True and save_eigval=True in ocelli.tl.MDM (this is done by default).
[3]:
oci.pp.neighbors(hf, n_neighbors=20, n_jobs=50, verbose=True)
oci.tl.modality_weights(hf, random_state=17, n_jobs=50, verbose=True)
oci.tl.MDM(hf, n_components=20, random_state=17, save_eigvec=True, save_eigval=True,
n_jobs=50, verbose=True)
oci.pp.neighbors(hf, x=['X_mdm'], n_neighbors=100, n_jobs=50, verbose=True)
oci.tl.transitions_graph(hf, x='X_mdm', transitions='velocity_graph', n_edges=3,
n_jobs=50, verbose=True)
oci.tl.fa2(hf, n_components=3, n_jobs=50, random_state=17)
oci.tl.projection(hf, x='X_fa2', phi=135, theta=135, random_state=17)
[rna] 20 nearest neighbors calculated.
[atac] 20 nearest neighbors calculated.
2025-01-29 18:37:43,131 INFO worker.py:1518 -- Started a local Ray instance.
Multimodal weights estimated.
2025-01-29 18:37:49,081 INFO worker.py:1518 -- Started a local Ray instance.
[rna] Unimodal Markov chain calculated.
[atac] Unimodal Markov chain calculated.
Multimodal Markov chain calculated.
Eigendecomposition finished.
20 Multimodal Diffusion Maps components calculated.
[X_mdm] 100 nearest neighbors calculated.
Transitions-based graph constructed.
Jan 29, 2025 6:38:02 PM org.netbeans.modules.masterfs.watcher.Watcher getNotifierForPlatform
INFO: Native file watcher is disabled
Jan 29, 2025 6:38:03 PM org.gephi.io.processor.plugin.DefaultProcessor process
INFO: # Nodes loaded: 7,160 (7,160 added)
Jan 29, 2025 6:38:03 PM org.gephi.io.processor.plugin.DefaultProcessor process
INFO: # Edges loaded: 21,480 (20,575 added)
*************************25%
*************************50%
*************************75%
*************************100%
Time = 38.478s
Below, you will impute values of RNA-seq gene expression and ATAC-seq promoter activity of two hair follicle’s Henle’s lineage markers, Krt71 and Krt73.
Imputation of RNA-seq gene expression#
Load raw RNA-seq gene expression of Krt71 and Krt73.
[4]:
rna = oci.read.h5ad('hf_shareseq_rna.h5ad')[:, ['Krt71', 'Krt73']]
# save FLE for plotting purposes
rna.obsm['X'] = hf.obsm['X_proj']
Plot the raw gene expression of Krt71.
[5]:
cmap = LinearSegmentedColormap.from_list('custom', ['#e3e3e3', '#ff0000', '#000000'], N=256)
rna.obs['Krt71_t0'] = rna[:, 'Krt71'].X.toarray().flatten()
oci.pl.scatter(rna, x='X', c='Krt71_t0', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[5]:
Plot the raw gene expression of Krt73.
[6]:
rna.obs['Krt73_t0'] = rna[:, 'Krt73'].X.toarray().flatten()
oci.pl.scatter(rna, x='X', c='Krt73_t0', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[6]:
Diffusion-based multimodal imputation with Ocelli requires multimodal eigenvectors and their corresponding eigenvalues. Add them to the rna object and impute rna.X using MDM-generated multimodal information.
[7]:
rna.uns['eigenvectors'] = hf.uns['eigenvectors']
rna.uns['eigenvalues'] = hf.uns['eigenvalues']
oci.tl.imputation(rna, t=1)
Plot the imputed gene expression of Krt71.
[8]:
rna.obs['Krt71_t1'] = rna[:, 'Krt71'].X.toarray().flatten()
oci.pl.scatter(rna, x='X', c='Krt71_t1', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[8]:
Plot the imputed gene expression of Krt73.
[9]:
rna.obs['Krt73_t1'] = rna[:, 'Krt73'].X.toarray().flatten()
oci.pl.scatter(rna, x='X', c='Krt73_t1', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[9]:
Ocelli’s multimodal imputation successfully strengthened the signal of both lowly expressed markers Krt71 and Krt73.
Imputation of ATAC-seq peak activity#
Repeat the above steps to impute the activity of Krt71 and Krt73 promoters. We found promoters of both markers using GREAT.
[10]:
promoters = {'Krt71': ['P_175393'], 'Krt73': ['P_337703']}
Load the raw promoter activity.
[11]:
peaks = oci.read.h5ad('hf_shareseq_atac_peaks.h5ad')[:, ['P_175393', 'P_337703']]
peaks.obsm['X'] = hf.obsm['X_proj']
Plot the raw activity of Krt71’s promoter.
[12]:
peaks.obs['Krt71_t0'] = peaks[:, promoters['Krt71']].X.toarray().flatten()
oci.pl.scatter(peaks, x='X', c='Krt71_t0', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[12]:
Plot the raw activity of Krt73’s promoter.
[13]:
peaks.obs['Krt73_t0'] = peaks[:, promoters['Krt73']].X.toarray().flatten()
oci.pl.scatter(peaks, x='X', c='Krt73_t0', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[13]:
Add eigenvectors and eigenvalues to the peaks obejct.
[14]:
peaks.uns['eigenvalues'] = hf.uns['eigenvalues']
peaks.uns['eigenvectors'] = hf.uns['eigenvectors']
oci.tl.imputation(peaks, t=1)
Plot the imputed activity of Krt71’s promoter.
[15]:
peaks.obs['Krt71_t1'] = peaks[:, promoters['Krt71']].X.toarray().flatten()
oci.pl.scatter(peaks, x='X', c='Krt71_t1', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[15]:
Plot the imputed activity of Krt73’s promoter.
[16]:
peaks.obs['Krt73_t1'] = peaks[:, promoters['Krt73']].X.toarray().flatten()
oci.pl.scatter(peaks, x='X', c='Krt73_t1', cmap=cmap, s=3, markerscale=2,
fontsize=15, figsize=(9, 6), title='')
[16]:
Ocelli successfully imputed chromatin accessibility levels at the promoters of Krt71 and Krt73 genes.
In Tutorial 4, you will explore a bimodal single-cell dataset of BMMCs profiled with ASAP-seq.