Intracellular Communication Modeling with scMagnify

Intracellular Communication Modeling with scMagnify#

Preliminaries#

In this tutorial, you will learn how to:

Cell-cell Communication analysis with LIANA+
Infer Dynamic Signaling-to-Transcription Axes. Correlate receptor expression with intracellular TF activity along a pseudotime trajectory.

Rationale#

Intercellular communication via ligand-receptor (L-R) interactions is fundamental to coordinating cellular responses in development and disease. However, a key challenge is understanding how these extracellular signals are dynamically translated into specific intracellular transcriptional programs to orchestrate cell state transitions [].

To address this, scMagnify provides a dynamic communication module. This analysis moves beyond simple L-R pairing and aims to connect intercellular signaling (receptors) directly to the intracellular TF activity along a defined differentiation trajectory.

Briefly, the workflow first establishes potential ligand-receptor-TF links from a knowledgebase. It then correlates receptor expression with TF activity across metacells ordered by pseudotime. Finally, a permutation test is employed to identify statistically significant and robust signaling-to-transcription axes that change dynamically during the process.

Import packages#

%load_ext autoreload
%autoreload 2

import warnings
from numba.core.errors import NumbaDeprecationWarning

warnings.simplefilter("ignore", category=NumbaDeprecationWarning)
warnings.simplefilter("ignore", FutureWarning)
warnings.simplefilter("ignore", UserWarning)
warnings.simplefilter("ignore", RuntimeWarning)

import os
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns

import scanpy as sc

import liana as li

import scmagnify as scm
from scmagnify.settings import settings

scm.info()

Installed version:	v0.0.0
Key dependencies:	scanpy v1.10.3, mudata v0.2.3, cellrank v2.0.7, decoupler v2.1.1, SEACells v0.3.3
PyTorch version:	v2.0.0+cu117
CUDA available:	True
scmagnify data cached:	True
Repository:	https://github.com/your-username/your-repo

Configurations#

scm.settings.verbosity = 2

%matplotlib inline

scm.settings.set_figure_params(
    dpi=100,
    facecolor="white",
    frameon=False,
)

scm.load_fonts(["Arial"])

plt.rcParams["font.family"] = "Arial"
plt.rcParams["grid.alpha"] = 0

# Setting a workspace
dirPjtHome = "/mnt/TrueNas/project/chenxufeng/Data/PMID38199997_NatCommun2024"
workDir = os.path.join(dirPjtHome, "scmagnify_wd")
scm.set_workspace(workDir)

workspace: /mnt/TrueNas/project/chenxufeng/Data/PMID38199997_NatCommun2024/scmagnify_wd/
├── data
├── models
├── tmpfiles
└── figures

scm.set_genome(version="hg38", genomes_dir="/home/chenxufeng/picb_cxf/Ref/human/hg38/")

                        Genome Information                        
┏━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Version ┃ Provider ┃ Directory                                 ┃
┡━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ hg38    │ UCSC     │ /home/chenxufeng/picb_cxf/Ref/human/hg38/ │
└─────────┴──────────┴───────────────────────────────────────────┘

Load the Data#

gdata = scm.read(os.path.join(settings.data_dir, "kidney-injury-tal_H11CORE.h5mu"))

gdata

Gene Regulatory Network (GRN) with 30523 edges.
MuData object with n_obs × n_vars = 8080 × 359416
  uns:	'attention_weights', 'filtered_network', 'motif_scan', 'network', 'peak_gene_corrs', 'regfactors', 'regfactors_colors'
  4 modalities
    RNA:	8080 x 22857
      obs:	'nCount_RNA', 'nFeature_RNA', 'library', 'percent.er', 'percent.mt', 'experiment', 'subclass.l3', 'subclass.l2', 'subclass.l1', 'nCount_ATAC', 'nFeature_ATAC', 'nucleosome_signal', 'nucleosome_percentile', 'TSS.enrichment', 'TSS.percentile', 'Total_fragments', 'FRiP', 'RNA.weight', 'ATAC.weight', 'dpt_pseudotime', 'celltype', 'leiden_res_0.50', 'n_counts', 'SEACell'
      var:	'name', 'n_cells', 'significant_genes', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
      uns:	'celltype_colors', 'celltype_sizes', 'diffmap_evals', 'draw_graph', 'hvg', 'iroot', 'leiden_res_0.50', 'leiden_res_0.50_colors', 'leiden_res_0.50_sizes', 'library_colors', 'log1p', 'neighbors', 'paga', 'subclass.l1_colors', 'subclass.l2_colors', 'subclass.l2_sizes', 'subclass.l3_colors', 'test_assoc', 'umap'
      obsm:	'X_diffmap', 'X_draw_graph_fa', 'X_lsi', 'X_pca', 'X_phate', 'X_umap', 'padj_mlm', 'score_mlm'
      varm:	'test_assoc_res'
      layers:	'counts', 'log1p_norm'
      obsp:	'connectivities', 'distances'
    ATAC:	8080 x 336500
      obs:	'nCount_RNA', 'nFeature_RNA', 'library', 'percent.er', 'percent.mt', 'experiment', 'subclass.l3', 'subclass.l2', 'subclass.l1', 'nCount_ATAC', 'nFeature_ATAC', 'nucleosome_signal', 'nucleosome_percentile', 'TSS.enrichment', 'TSS.percentile', 'Total_fragments', 'FRiP', 'RNA.weight', 'ATAC.weight', 'SEACell'
      var:	'count', 'percentile', 'AA', 'AC', 'AG', 'AT', 'CA', 'CC', 'CG', 'CT', 'GA', 'GC', 'GG', 'GT', 'TA', 'TC', 'TG', 'TT', 'GC.percent', 'sequence.length'
      uns:	'library_colors', 'neighbors', 'peak_seq', 'subclass.l2_colors', 'subclass.l3_colors', 'umap'
      obsm:	'X_lsi', 'X_pca', 'X_svd', 'X_umap'
      layers:	'counts'
      obsp:	'connectivities', 'distances'
    GRN:	8080 x 54
      obs:	'nCount_RNA', 'nFeature_RNA', 'library', 'percent.er', 'percent.mt', 'experiment', 'subclass.l3', 'subclass.l2', 'subclass.l1', 'nCount_ATAC', 'nFeature_ATAC', 'nucleosome_signal', 'nucleosome_percentile', 'TSS.enrichment', 'TSS.percentile', 'Total_fragments', 'FRiP', 'RNA.weight', 'ATAC.weight', 'dpt_pseudotime', 'celltype', 'leiden_res_0.50', 'n_counts', 'SEACell'
      var:	'mean_activity'
      uns:	'basal_grn', 'celltype_colors', 'celltype_sizes', 'diffmap_evals', 'draw_graph', 'hvg', 'iroot', 'leiden_res_0.50', 'leiden_res_0.50_colors', 'leiden_res_0.50_sizes', 'library_colors', 'log1p', 'neighbors', 'paga', 'subclass.l1_colors', 'subclass.l2_colors', 'subclass.l2_sizes', 'subclass.l3_colors', 'test_assoc', 'umap'
      obsm:	'X_diffmap', 'X_draw_graph_fa', 'X_lsi', 'X_pca', 'X_phate', 'X_umap', 'padj_mlm', 'score_mlm'
      varm:	'network_score'
    RegFactor:	8080 x 5
      obs:	'nCount_RNA', 'nFeature_RNA', 'library', 'percent.er', 'percent.mt', 'experiment', 'subclass.l3', 'subclass.l2', 'subclass.l1', 'nCount_ATAC', 'nFeature_ATAC', 'nucleosome_signal', 'nucleosome_percentile', 'TSS.enrichment', 'TSS.percentile', 'Total_fragments', 'FRiP', 'RNA.weight', 'ATAC.weight', 'dpt_pseudotime', 'celltype', 'leiden_res_0.50', 'n_counts', 'SEACell'
      uns:	'celltype_colors', 'celltype_sizes', 'diffmap_evals', 'draw_graph', 'hvg', 'iroot', 'leiden_res_0.50', 'leiden_res_0.50_colors', 'leiden_res_0.50_sizes', 'library_colors', 'log1p', 'neighbors', 'paga', 'subclass.l1_colors', 'subclass.l2_colors', 'subclass.l2_sizes', 'subclass.l3_colors', 'test_assoc', 'umap'
      obsm:	'X_diffmap', 'X_draw_graph_fa', 'X_lsi', 'X_pca', 'X_phate', 'X_umap', 'padj_mlm', 'score_mlm'
      varm:	'Lag_loadings', 'TF_loadings', 'TG_loadings'

sc.pl.umap(
    gdata["RNA"],
    color=["celltype"],
    size=10,
    frameon=False,
)

../_images/067a1c11edbd194e3ca10c9560632b637c3829f03ce82fbbe2f09dc750c11d29.png

scm.pl.stripplot(gdata, sortby="degree_centrality", n_top=30)

../_images/ddbb29c9eb19121b2db5c53c70491190fc0c2158c716f2c325eec06a92d22fdc.png

adata_cci = sc.read(
    "/mnt/TrueNas/project/chenxufeng/Data/PMID38199997_NatCommun2024/1_AnnData/kidney-injury-rna_tal_imm.h5ad"
)

adata_cci

AnnData object with n_obs × n_vars = 12813 × 36554
    obs: 'nCount_RNA', 'nFeature_RNA', 'library', 'percent.er', 'percent.mt', 'experiment', 'subclass.l3', 'subclass.l2', 'subclass.l1', 'nCount_ATAC', 'nFeature_ATAC', 'nucleosome_signal', 'nucleosome_percentile', 'TSS.enrichment', 'TSS.percentile', 'Total_fragments', 'FRiP', 'RNA.weight', 'ATAC.weight', 'celltype_hierarchical', 'celltype'
    var: 'name'
    uns: 'celltype_colors', 'dendrogram_celltype_hierarchical', 'library_colors', 'log1p', 'neighbors', 'subclass.l1_colors', 'subclass.l2_colors', 'subclass.l3_colors', 'umap'
    obsm: 'X_lsi', 'X_pca', 'X_umap'
    layers: 'log1p_norm'
    obsp: 'connectivities', 'distances'

sc.pl.umap(adata_cci, color=["celltype", "library", "subclass.l1"], wspace=0.4, ncols=3)

../_images/8f125f15645f5212af7dc5a9b5f8f8c039c90c5919292db9f787dba0b2f733cf.png

meta_mdata = scm.read(os.path.join(settings.data_dir, "kidney-injury-tal_metacells.h5mu"))
meta_mdata["RNA"].layers["log1p_norm"] = meta_mdata["RNA"].X.copy()

RegFactor Analysis#

scm.pl.barplot(
    gdata,
    modal="RegFactor",
    key="TF_loadings",
    swap_df=True,
    n_top=10,
    ncols=5,
    cmap="Blues_r",
)

../_images/9e53c9d6891fffecadf3078db17607239917ba216287b34006862592303f3b08.png

scm.pl.barplot(
    gdata,
    modal="RegFactor",
    key="TG_loadings",
    swap_df=True,
    n_top=10,
    ncols=5,
    cmap="Blues_r",
)

../_images/0a8acdc8dfcb33797491c8dbcade797b1f4487ac5fd260d493ea5bbb1c656ae9.png

sc.pl.violin(
    gdata["RegFactor"], keys=gdata["RegFactor"].var_names, groupby="celltype", rotation=45, stripplot=False, show=True
)

../_images/ccb1352b8d90894cabc812e41494d2f26e994bd5d40885fc7cec8a4053c90759.png

Cell-cell Communication with Liana+#

li.mt.cellchat(
    adata_cci,
    groupby="celltype",
    resource_name="cellchatdb",
    verbose=True,
    use_raw=False,
    layer="log1p_norm",
    key_added="cellchat_res",
)

cellchat_res = adata_cci.uns["cellchat_res"].copy()

Generating ligand-receptor stats for 12813 samples and 878 features

from scmagnify.external.plotting.liana import LianaVisualizer

lvis = LianaVisualizer(
    adata_cci, res_key="cellchat_res", magnitude_col="lr_probs", pvalue_col="cellchat_pvals", cluster_key="celltype"
)

fig = lvis.plot_chord(
    kind="count", normalize="row", link_kws={"ec": "black", "lw": 0, "direction": 1}, label_kws={"size": 15}
)

../_images/fdcdb88822db5030e14723b07c67c5027a23aaa2f0b56866ef15069d3ac87a0c.png

fig, ax = lvis.plot_interact_heatmap(cmap="Reds")

../_images/ea6e9551fe5418b713687d26083b52c715893e35a12b5cd87c604662f7c87d5b.png

fig, ax = lvis.plot_radar(cell="aTAL", mode="incoming", kind="strength")

../_images/90983efb1ef26e0c4598e7039ef9d14629b20f61f662409b86c603be27bb31d8.png

fig, ax = lvis.plot_radar(cell="MYOF", mode="outgoing", kind="strength")

../_images/d6104caaf8260882de65ab7bd04f5e1c8ed3b47f25bcd9962eb02804f16fb0b1.png

fig, ax = lvis.plot_radar(mode="incoming", kind="strength")