脚本更新---高精度（Xenium、CosMx）的邻域通讯分析

作者，Evil Genius

维权路漫漫，且行且珍惜。

今天我们更新脚本，高精度（Xenium、CosMx）的邻域通讯分析

python代码，python大家尽量学一学，如何搭建python环境以及写封装脚本在linux系统上运行，是空转必修课，大家可以参考2025番外--linux、R、python培训.

import stlearn as st
import pathlib as pathlib
import matplotlib.pyplot as plt
import numpy as np
import random as random
import os as os

st.settings.set_figure_params(dpi=120)

# Make sure all the seeds are set
seed = 0
np.random.seed(seed)
random.seed(seed)
os.environ['PYTHONHASHSEED'] = str(seed)

# Ignore all warnings
import warnings
warnings.filterwarnings("ignore")

###Loading the data
# Setup download directory and get data
st.settings.datasetdir =  pathlib.Path.cwd().parent / "data"
library_id = "Xenium_FFPE_Human_Breast_Cancer_Rep1"
data_dir = st.settings.datasetdir / "Xenium_FFPE_Human_Breast_Cancer_Rep1"

st.datasets.xenium_sge(library_id=library_id, include_hires_tiff=True)

adata = st.ReadXenium(feature_cell_matrix_file=data_dir / "cell_feature_matrix.h5",
                      cell_summary_file=data_dir / "cells.csv.gz",
                      library_id=library_id,
                      image_path=data_dir / "he_image.ome.tif",
                      scale=1,
                      spot_diameter_fullres=15,
                      alignment_matrix_file=data_dir / "he_imagealignment.csv",
                      experiment_xenium_file=data_dir / "experiment.xenium",
                      )

st.pp.filter_genes(adata, min_counts=10)
st.pp.filter_cells(adata, min_counts=10)

adata.raw = adata

# Run PCA, neighbors and clustering.
st.em.run_pca(adata, n_comps=50, random_state=0)
st.pp.neighbors(adata, n_neighbors=25, use_rep='X_pca', random_state=0)
st.tl.clustering.louvain(adata, random_state=0)

st.pl.cluster_plot(adata, use_label="louvain", image_alpha=0, size=4, figsize=(10, 10))

识别hotspot区域

st.pp.normalize_total(adata)

### Gridding.
grid = st.tl.cci.grid(adata, n_row=n_, n_col=n_, use_label='louvain')
print(grid.shape)  # Slightly less than the above calculation, since we filter out spots with 0 cells.

fig, axes = plt.subplots(ncols=2, figsize=(20, 8))
st.pl.cluster_plot(grid, use_label='louvain', size=10, ax=axes[0], show_plot=False)
st.pl.cluster_plot(adata, use_label='louvain', ax=axes[1], show_plot=False)
axes[0].set_title(f'Grid louvain dominant spots')
axes[1].set_title(f'Cell louvain labels')
plt.show()

groups = list(grid.obs['louvain'].cat.categories)
for group in [groups[3], groups[0]]:
    fig, axes = plt.subplots(ncols=3, figsize=(20, 8))
    group_props = grid.uns['louvain'][group].values
    grid.obs['group'] = group_props
    st.pl.feat_plot(grid, feature='group', ax=axes[0], show_plot=False, vmax=1, show_color_bar=False)
    st.pl.cluster_plot(grid, use_label='louvain', list_clusters=[group], ax=axes[1], show_plot=False)
    st.pl.cluster_plot(adata, use_label='louvain', list_clusters=[group], ax=axes[2], show_plot=False)
    axes[0].set_title(f'Grid {group} proportions (max = 1)')
    axes[1].set_title(f'Grid {group} max spots')
    axes[2].set_title(f'Individual cell {group}')
    plt.show()

fig, axes = plt.subplots(ncols=2, figsize=(20, 5))
st.pl.gene_plot(grid, gene_symbols='CXCL12', ax=axes[0], show_color_bar=False, show_plot=False)
st.pl.gene_plot(adata, gene_symbols='CXCL12', ax=axes[1], show_color_bar=False, show_plot=False, vmax=80)
axes[0].set_title(f'Grid CXCL12 expression')
axes[1].set_title(f'Cell CXLC12 expression')
plt.show()

LR Permutation Test

lrs = st.tl.cci.load_lrs(['connectomeDB2020_lit'], species='human')
print(len(lrs))

st.tl.cci.run(grid, lrs,
              min_spots=20,  # Filter out any LR pairs with no scores for less than min_spots
              distance=250,  # None defaults to spot+immediate neighbours; distance=0 for within-spot mode
              n_pairs=1000,  # Number of random pairs to generate; low as example, recommend ~10,000
              n_cpus=None,   # Number of CPUs for parallel. If None, detects & use all available.
              )

lr_info = grid.uns['lr_summary']

best_lr = grid.uns['lr_summary'].index.values[0]

stats = ['lr_scores', '-log10(p_adjs)', 'lr_sig_scores']
fig, axes = plt.subplots(ncols=len(stats), figsize=(12, 6))
for i, stat in enumerate(stats):
    st.pl.lr_result_plot(grid, use_result=stat, use_lr=best_lr, show_color_bar=False, ax=axes[i])
    axes[i].set_title(f'{best_lr} {stat}')

生活很好，有你更好