作者，Evil Genius

今天我们复习。

需要实现的目标

1、目标细胞的邻域组成

2、不同组的邻域差异

3、空间特征分数

来吧， 大家一定要好好学习，我经历过了，深知一篇好的文章有多么重要,现在努力一分产出2分，随着越来越卷，以后努力1分有没有产出都是问题。

第一部分

import anndata as ad
import pandas as pd
import cellcharter as cc
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import scanpy as sc
import squidpy as sq

adata = ad.read_h5ad('cosmx.h5ad')

adata.uns['spatial'] = {s: {} for s in adata.obs['sample'].unique()}

adata_sample = adata[adata.obs['patient'] == 'Lung5'] ###选择其中一个样本
cc.gr.nhood_enrichment(
    adata_sample,
    cluster_key = 'cluster_20'
)

cc.pl.nhood_enrichment(
    adata_sample,
    cluster_key='cluster_20',
    row_groups=[4, 7],
    col_groups=[1, 8, 10, 11, 16, 17],
    vmin=-0.5, vmax=0.5,
    min_freq=0.001,
    transpose=True
)

signatures = pd.read_excel('data/tumor_cell_state_signatures_gavish.xlsx')
signatures.columns = [c.replace('/', '').replace('  ', '_').replace(' - ', '-').replace(' ', '_') for c in signatures.columns]
signatures_dict = signatures.to_dict(orient='list')

signatures_dict['Cell_Proliferation'] = [x for x in pd.concat((
        signatures['MP1_Cell_Cycle-G2M'],
        signatures['MP2_Cell_Cycle-G1S'],
        signatures['MP3_Cell_Cylce_HMG-rich'],
        signatures['MP4_Cell_Cylce_Chromatin_']
    ))
               if x in adata.var_names]

signatures_dict['GO_Cytokine_Cytokine_Receptor_Interaction']= ['VEGFA', 'CXCL1', 'CXCL3', 'CXCR4', 'PDGFA', 'CXCL8', 'CXCL2', 'TNF', 'CCL20']

def score_cells(adata, signature_name, genes):
    adata.obs[f'score_{signature_name}'] = np.nan
    for patient in adata.obs['patient'].cat.categories:
        adata_patient = adata[adata.obs['patient'] == patient].copy()
        sc.tl.score_genes(adata_patient, gene_list=[x for x in genes if x in adata_patient.var_names], score_name=f'score_{signature_name}')
        adata.obs[f'score_{signature_name}'][adata.obs['patient'] == patient] = adata_patient.obs[f'score_{signature_name}']

def smooth_signature(adata, signature_name, n_neighbors, group_key=None, groups=None):
    if groups:
        adata_groups = adata[adata.obs[group_key].isin(groups)]
    sq.gr.spatial_neighbors(adata_groups, library_key='sample', n_neighs=n_neighbors, coord_type='generic')
    adata_groups.obsp['spatial_connectivities'].setdiag(1)
    
    adata.obs[f'score_{signature_name}_smoothed'] = np.nan
    for patient in adata_groups.obs['patient'].cat.categories:
        adata_groups_patient = adata_groups[adata_groups.obs['patient'] == patient]
        adata_groups_patient.obs[f'score_{signature_name}_smoothed'] = np.nan*np.ones(adata_groups_patient.shape[0])
        adj = adata_groups_patient.obsp['spatial_connectivities']
        score_smoothed = (adj @ adata_groups_patient.obs[f'score_{signature_name}']) / np.array(np.sum(adj, axis=1)).squeeze()
        score_smoothed[score_smoothed == float('inf')] = 0

        low = np.nanpercentile(score_smoothed, 5)
        score_smoothed[score_smoothed < low] = low

        high = np.nanpercentile(score_smoothed, 95)
        score_smoothed[score_smoothed > high] = high

        adata.obs[f'score_{signature_name}_smoothed'][(adata.obs['patient'] == patient) & (adata.obs[group_key].isin(groups))] = score_smoothed

signatures_names  = ['GO_Cytokine_Cytokine_Receptor_Interaction', 'MP13_EMT-II', 'MP6_Hypoxia', 'Cell_Proliferation']

for SIGNATURE_NAME in signatures_names:
    score_cells(adata, SIGNATURE_NAME, signatures_dict[SIGNATURE_NAME])
    smooth_signature(adata, SIGNATURE_NAME, 50, 'cell_type_reduced', ['tumor'])

adata_sample = adata[adata.obs['sample'] == 'LUAD-9 R1'].copy()
adata_sample.obsm['spatial'][:, 1] = np.max(adata_sample.obsm['spatial'][:, 1]) - adata_sample.obsm['spatial'][:, 1]

sq.pl.spatial_scatter(
    adata_sample,
    color=[f'score_{signature}_smoothed' for signature in signatures_names],
    cmap='Spectral_r',
    library_key='sample',
    library_id='LUAD-9 R1',
    img=None,
    size=30,
    na_color="#E5E5E5",
    ncols=4,
    title=['Cytoking-Cytokine interation', 'EMT', 'Response to Hypoxia', 'Cell Proliferation'],
    legend_loc=False,
)
plt.show()

第二部分

from matplotlib import cm
import anndata as ad
import squidpy as sq
import cellcharter as cc
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

adata = ad.read_h5ad('data/codex.h5ad')

adata.obs['condition'] = adata.obs['sample'].str.split('-').str[0].astype('category')

adata.obs_names_make_unique()

adata.obs['spatial_cluster'] = adata.obs['cluster_11'].map(cluster2region)
adata.obs['spatial_cluster'] = adata.obs['spatial_cluster'].astype('category')
adata.obs['spatial_cluster'] = adata.obs['spatial_cluster'].cat.reorder_categories(list(cluster2region.values()))

adata_balbc = adata[adata.obs['condition'] == 'BALBc']
cc.gr.nhood_enrichment(
    adata_balbc,
    cluster_key='spatial_cluster',
)

cc.pl.nhood_enrichment(
    adata_balbc,
    cluster_key='spatial_cluster',
    vmin=-1,
    vmax=1,
    title='Healthy spleen',
    figsize=(5, 5),
)

adata_balbc = adata[adata.obs['condition'] == 'MRL']
cc.gr.nhood_enrichment(
    adata_balbc,
    cluster_key='spatial_cluster',
)

cc.pl.nhood_enrichment(
    adata_balbc,
    cluster_key='spatial_cluster',
    vmin=-1,
    vmax=1,
    title='Systemic lupus',
    figsize=(5, 5),
)

cc.gr.diff_nhood_enrichment(
    adata,
    cluster_key='spatial_cluster',
    condition_key='condition',
    condition_groups=['MRL', 'BALBc'],
    library_key='sample',
)

cc.pl.diff_nhood_enrichment(
    adata,
    cluster_key='spatial_cluster',
    condition_key='condition',
    condition_groups=['MRL', 'BALBc'],
    annotate=True,
    figsize=(7,7),
)

生活很好，有你更好