In command line / or jupyter notebook

pip install liana[common]

pip install adjustText==0.8

import numpy as np

import liana as li

import mudata as mu

import scanpy as sc

from matplotlib import pyplot as plt

from adjustText import adjust_text

save figure parameters

kwargs = {'frameon':False, 'size':1.5, 'img_key':'lowres'}

let's download the data

rna = sc.read("sma_rna.h5ad", backup_url="https://figshare.com/ndownloader/files/44624974?private_link=4744950f8768d5c8f68c")

msi = sc.read("sma_msi.h5ad", backup_url="https://figshare.com/ndownloader/files/44624971?private_link=4744950f8768d5c8f68c")

ct = sc.read("sma_ct.h5ad", backup_url="https://figshare.com/ndownloader/files/44624968?private_link=4744950f8768d5c8f68c")

and create a MuData object

mdata = mu.MuData({'rna':rna, 'msi':msi, 'ct':ct})

# examine the object

mdata

examine data

fig, axes = plt.subplots(1, 3, figsize=(12, 3))

sc.pl.spatial(rna, color='log1p_total_counts', ax=axes[0], **kwargs, show=False)

sc.pl.spatial(msi, color='Dopamine', cmap='magma', ax=axes[1], **kwargs, show=False)

sc.pl.spatial(ct, color='MSN1', cmap='viridis', ax=axes[2], **kwargs, show=False)

fig.subplots_adjust(wspace=0, hspace=0)

fig.tight_layout()

examine experimental design

sc.pl.spatial(rna, color=['lesion', 'region'], **kwargs, wspace=0.25)

sc.pp.highly_variable_genes(rna, flavor='cell_ranger', n_top_genes=5000)

sc.pp.highly_variable_genes(msi, flavor='cell_ranger', n_top_genes=150)

ct.var['cv'] = ct.X.A.var(axis=0) / ct.X.A.mean(axis=0)

ct.var['highly_variable'] = ct.var['cv'] > np.percentile(ct.var['cv'], 20)

msi = msi[:, msi.var['highly_variable']]

rna = rna[:, rna.var['highly_variable']]

ct = ct[:, ct.var['highly_variable']]

# scale msi intensities

sc.pp.scale(msi, max_value=5)

metalinks = li.rs.get_metalinks(tissue_location='Brain',

                            biospecimen_location='Cerebrospinal Fluid (CSF)',

                            source=['CellPhoneDB', 'NeuronChat']

                            )

metalinks.head()

Obtain Human to Murine gene symbols Map

map_df = li.rs.get_hcop_orthologs(columns=['human_symbol', 'mouse_symbol'],

                              min_evidence=3

                              ).rename(columns={'human_symbol':'source',

                                               'mouse_symbol':'target'})

metalinks = li.rs.translate_column(resource=metalinks,

                               map_df=map_df,

                               column='gene_symbol',

                               one_to_many=1)

metalinks.head()

Intersect the RNA modality with the receptors

receptors = np.intersect1d(metalinks['gene_symbol'].unique(), rna.var_names)

rec = rna[:, receptors].copy()

We use the metabolite modality as a reference to which we align the other modalities.

We use the spatial_neighbors function to compute the spatial proximity from cell types and brain-specific receptors to the metabolite intensities.

plot, _ = li.ut.query_bandwidth(coordinates=rna.obsm['spatial'], start=0, end=1500)

plot

choose bandwidth

bandwidth = 500

cutoff = 0.1

# distances of metabolties to RNA

reference = mdata.mod["msi"].obsm["spatial"]

Compute proximities to the reference for each modality

li.ut.spatial_neighbors(ct, bandwidth=bandwidth, cutoff=cutoff, spatial_key="spatial", reference=reference, set_diag=False, standardize=False)

li.ut.spatial_neighbors(rec, bandwidth=bandwidth, cutoff=cutoff, spatial_key="spatial", reference=reference, set_diag=False, standardize=False)

MISTy

mdata.update_obs()

misty = li.mt.MistyData({"intra": msi, "receptor": rec, "ct": ct}, enforce_obs=False, obs=mdata.obs)

misty

Learn the relationships between intra and the extra views

misty(model=li.mt.sp.LinearModel, verbose=True, bypass_intra=True, maskby='lesion')

li.pl.target_metrics(misty, stat='multi_R2', return_fig=True, top_n=20, filter_fun=lambda x: x['intra_group']=='intact')

interactions = misty.uns['interactions']

interactions = interactions[(interactions['intra_group'] == 'intact') & (interactions['target'] == 'Dopamine')]

# Create scatter plot

plt.figure(figsize=(5, 4))

# rank rank by abs importances

interactions['rank'] = interactions['importances'].rank(ascending=False)

plt.scatter(interactions['rank'], interactions['importances'], s=11,

        c=interactions['view'].map({'ct': '\#008B8B', 'receptor': '\#a11838'}))

        

# add for top 10

top_n = interactions[interactions['rank'] <= 10]

texts = []

for i, row in top_n.iterrows():

texts.append(plt.text(row['rank'], row['importances'], row['predictor'], fontsize=10))

adjust_text(texts, arrowprops=dict(arrowstyle="->", color='grey', lw=1.5))

plt.tight_layout()

Focusing on Dopamine, we can next use LIANA+'s local metrics to identify the subregions of interactions with MSN1/2 cells.

# to do so, we need to interpolate one modality to the other, such that they are on the same coordinate system, while this is done internally for the multi-view learning, we need to interpolate it as:

metabs = li.ut.interpolate_adata(target=msi, reference=rna, use_raw=False, spatial_key='spatial')

let's rebuild this with the update modalities:

mdata = mu.MuData({'msi': metabs, 'rna':rna, 'deconv':ct}, obsm=rna.obsm, obs=rna.obs, uns=rna.uns)

re-calculate neighbours

li.ut.spatial_neighbors(mdata, bandwidth=bandwidth, cutoff=cutoff, set_diag=True)

define interactions of interest

interactions = metalinks[['metabolite', 'gene_symbol']].apply(tuple, axis=1).tolist()

Let's calculate the local metrics for the Dopamine intensities with and Drd1/2 receptors.

li.mt.bivariate(mdata,

            local_name='cosine',

            x_mod='msi', 

            y_mod='rna',

            key_added='lr',

            x_use_raw=False, 

            y_use_raw=False,

            verbose=True, 

            mask_negatives=True, 

            n_perms=1000,

            interactions=interactions,

            x_transform=sc.pp.scale,

            y_transform=sc.pp.scale,

        )

sc.pl.spatial(mdata.mod['lr'],

          color=['Dopamine^Drd1', 'Dopamine^Drd2'],

          cmap='cividis_r', vmax=1, layer='pvals',

          **kwargs)