Xenium多样本分析实战-代码分享

生信大杂烩

发布于 2025-06-15 11:52:06

21100

代码可运行

运行总次数：0

代码可运行

最近陆陆续续完成了几个Xenium项目的数据分析，主要集中在肿瘤样本研究，每个项目由于实验设计的不同，老师关心的科学问题不同，都做了许多个性化的分析，有时间闲下来后自己做了一些总结，将自己在Xenium, Visium HD 以及CosMx这些高精度空间数据分析过程中用到的一些分析方法及代码(这里分享的主要是基于python实现的各种方法，之前一直使用的是R进行分析，数据量大了之后非常慢，于是乎转向了python

)和大家进行分享，希望和大家相互学习。

最近也经常有同学微信问各种分析怎样实现的，个人觉得生信分析的代码实现应该是最容易的，因为除了搞算法开发的，大家都是调包侠，真正厉害的是那些能设计实验、从分析结果中找到有意思的东西、最后讲个完整生物学故事的人，所以希望和大家相互学习。

一、基础分析

1. Xenium Ranger数据读取

详细内容参考我们之前的相关介绍：

Xenium数据分析 | 下机数据读取

2. 低质量细胞过滤

和单细胞数据分析一样，过滤低质量细胞是至关重要的一步预处理操作，提高分析的准确性和可靠性。过滤掉这些技术噪音和低质量数据点后，剩下的细胞更能代表真实的、状态良好的细胞群体。这使得后续分析的生物学信号更清晰，结果更可信，更容易揭示真实的生物学变异（如细胞类型、状态、发育轨迹、对刺激的反应等）

使用中位数绝对偏差(Median Absolute Deviation, MAD)识别低质量细胞：

MAD是一种鲁棒性统计量，对异常值不敏感（异常值不会过度影响MAD）、适用于非正态分布数据、在存在极端值时更稳健，相比标准差更适合存在异常值的数据。
计算步骤：

1. 计算所有数据点与中位数的绝对偏差：|Xᵢ - median(X)|

2. 取这些绝对偏差的中位数：median(|Xᵢ - median(X)|)

nmads的选择

1. 默认3：较严格（约排除~0.3%数据，假设正态分布）

2. 常用5：较宽松（单细胞分析中常用）

3.需结合数据分布调整

from scipy.stats import median_abs_deviation
def is_outlier(adata, metric, nmads=6):
    M = adata.obs[metric]
    outlier = (M < np.median(M) - nmads * median_abs_deviation(M)) | (M > np.median(M) + nmads * median_abs_deviation(M))
    return outlier
sc.pp.calculate_qc_metrics(adata, percent_top=(10, 50, 100), inplace=True)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts'], groupby='sample', jitter=False, rotation=90, multi_panel=True, show=False)
adata.obs["outlier"] = (
        is_outlier(adata, 'total_counts', 6) |
        is_outlier(adata, 'n_genes_by_counts', 6)
    )
adata = adata[~adata.obs['outlier'], :]
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts'], groupby='sample', jitter=False, rotation=90, multi_panel=True, show=False)

Note: Xenium细胞过滤还要考虑另一个指标，细胞核面积大小（因为Xenium有部分细胞的分割是基于识别出细胞核后，向外扩5um当做细胞边界，有些细胞核面积异常大的，这些就是分割有误的细胞，需要删除），下面展示的是细胞核面积与细胞面积比的统计，需要根据实际情况将这部分考虑进细胞过滤条件。

sns.histplot(
    adata.obs["nucleus_area"] / adata.obs["cell_area"],
    kde=False
)

3. 归一化和特征选择

为什么要归一化

不同细胞捕获的 mRNA 总量不同（称为"测序深度"差异），这会导致捕获细胞基因数量多的细胞所有基因表达量都偏高；
标准化可以使所有细胞的总表达量处于同一水平，解决不同细胞在文库制备、测序过程中的技术差异，标准化后细胞间的差异主要反映生物学差异而非技术差异；
使数据分布更合理，原始UMI计数数据是离散的，标准化为后续分析（如PCA、聚类）奠定基础。

常用方法为对每个细胞中的counts深度缩放，使用sc.pp.normalize_total函数，target_sum=1e4参数意味着将每个细胞的总表达量缩放到10,000，随后进行 log1p 转换。Note: target_sum参数如果不设置，会将每个细胞的总表达量缩放到所有细胞counts的中值

什么是特征选择

为了降低单细胞数据的维度，仅包含信息量最大的基因，这一步骤通常称为特征选择(选择高变基因)，这样可以降低数据复杂度，加快计算；
scanpy通过sc.pp.highly_variable_genes函数选择高变基因。

adata.layers['counts'] = adata.X.copy()
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key="sample")
adata = adata[:, adata.var['highly_variable']]

4. 降维聚类

使用Banksy进行降维聚类，传统的单细胞降维聚类，只考虑了细胞中的基因表达情况，忽略了细胞的空间位置。Banky可以调整lambda参数增加空间的特征权重进行降维聚类，更适用于空间转录组数据进行降维聚类。当数据中有多个样本的时候，我们还需要使用harmony进行批次校正。

详细内容参考我们之前的相关介绍：

Visium HD数据分析之空间聚类算法Banksy

详解Banksy多样本空间聚类分析

5. 细胞类型注释

降维聚类结束后，计算每个Cluster中的marker基因，使用这些marker基因进行细胞类型注释。

sc.tl.rank_genes_groups(
    adata, 
    groupby='leiden', 
    # use_raw=True, 
    method="t-test", # method='t-test','wilcoxon',key_added='rank_genes_groups' 
    pts=True
)
sc.pl.rank_genes_groups_dotplot(
    adata, 
    groupby='leiden', 
    standard_scale="var", 
    n_genes=5,
    cmap='RdYlBu_r',
    show=False,
)

enhance_umap(
    adata, 
    color='cell_type_clean', 
    figsize=(8,6)
)
plot_cellular_composition(
    adata, 
    column='cell_type_clean', 
    group='sample', 
    max_cols=4, 
    label_threshold=2, 
    show_labels=True, 
    adjust_labels=False, 
    plot_type="bar"
)

二、空间高级分析

1. 细胞微环境的特征、结构和功能

细胞微环境，也称为空间域、细胞邻域、细胞niche等，是指围绕细胞的局部环境，在决定各种生物过程中起着关键作用，例如维持组织稳态和塑造疾病进展等。空间组学技术的最新进展，提供了单细胞分辨率的分子特征，允许在组织背景下系统地探索细胞状态、功能和相互作用。

现在已经有多种计算方法，将细胞本身的特征(例如细胞的分子谱)及其微环境的各种特征(例如细胞邻域的细胞组成或分子谱)以统一的方式定义为来自不同“视角”的细胞特征，通过整合细胞的分子特征与空间信息来识别细胞微环境

结合细胞基因表达和周围细胞空间特征：

Bansksy(https://prabhakarlab.github.io/Banksy/)通过结合细胞自身的分子特征与其邻居的分子特征生成邻居增强特征，并提供一个特定的超参数来调整细胞及其局部微环境的贡献;
CellCharter(https://cellcharter.readthedocs.io/en/latest/)通过线性加权和小区域聚合将邻居的分子特征整合到细胞自身的分子特征中(空间生态位差异分析);

只考虑细胞类型和空间位置，之前过单样本的计算方法：

Cellular Neighborhood（Xenium | 细胞邻域(Cellular Neighborhood)分析(fixed k-nearest neighbor) Xenium | 细胞邻域(Cellular Neighborhood)分析(fixed radius)）
Cell Niche(空间转录组 | 细胞niche分析)

以下是多样本CN统计代码更新：

def cell_neighborhood(
        adata, 
        mode='knn', 
        sample='sample', 
        cluster='celltype', 
        radius=80, 
        k_neighbor=20, 
        n_clusters=10,
        out_dir=None):
    
    print(f"开始细胞邻域分析，共 {adata.n_obs} 个细胞，{len(adata.obs[sample].unique())} 个样本")
    start_time = time.time()
    
    # 1. 预计算所有细胞类型
    all_cell_types = np.sort(adata.obs[cluster].unique())
    cell_type_to_idx = {ct: i for i, ct in enumerate(all_cell_types)}
    n_cell_types = len(all_cell_types)
    
    # 2. 准备结果存储
    sample_data = []
    samples = adata.obs[sample].unique()
    
    # 3. 并行处理每个样本
    for s in samples:
        sample_start = time.time()
        # 3.1 提取样本数据
        sample_mask = adata.obs[sample] == s
        tmp = adata[sample_mask]
        x_y_coordinates = tmp.obsm["spatial"]
        n_cells = x_y_coordinates.shape[0]
        cts = tmp.obs[cluster].values
        
        # 3.2 构建邻域矩阵
        if mode == 'radius':
            connectivity_matrix = radius_neighbors_graph(
                x_y_coordinates, radius, mode='connectivity', include_self=False)
        else:  # 默认为 knn
            connectivity_matrix = kneighbors_graph(
                x_y_coordinates, k_neighbor, mode='connectivity', include_self=False)
        
        # 3.3 创建稀疏的one-hot编码矩阵
        row_indices = np.arange(n_cells)
        col_indices = np.array([cell_type_to_idx[ct] for ct in cts])
        data = np.ones(n_cells)
        one_hot_sparse = csr_matrix((data, (row_indices, col_indices)), 
                                    shape=(n_cells, n_cell_types))
        
        # 3.4 计算邻域组成矩阵
        # 使用稀疏矩阵乘法提高效率
        cell_by_neighborhood_matrix = connectivity_matrix.dot(one_hot_sparse)
        
        # 3.5 标准化矩阵
        row_sums = np.array(cell_by_neighborhood_matrix.sum(axis=1)).flatten()
        # 避免除以零
        row_sums[row_sums == 0] = 1
        # 使用稀疏矩阵除法
        normalized_matrix = cell_by_neighborhood_matrix.multiply(1 / row_sums[:, np.newaxis])
        
        # 3.6 转换为DataFrame
        cell_by_neighborhood_df = pd.DataFrame.sparse.from_spmatrix(
            normalized_matrix, columns=all_cell_types)
        
        # 3.7 添加元数据
        cell_by_neighborhood_df["cell_type"] = cts
        cell_by_neighborhood_df[["x", "y"]] = x_y_coordinates
        cell_by_neighborhood_df[sample] = s
        
        sample_data.append(cell_by_neighborhood_df)
        print(f"样本 {s} 处理完成 ({n_cells} 个细胞), 耗时: {time.time()-sample_start:.2f}秒")
    
    # 4. 合并所有样本数据
    all_neighbors_df = pd.concat(sample_data, ignore_index=True)
    print(f"合并所有样本数据, 总耗时: {time.time()-start_time:.2f}秒")
    
    # 5. 提取特征矩阵
    feature_cols = all_cell_types.tolist()
    values = all_neighbors_df[feature_cols].sparse.to_dense().values
    
    # 6. KMeans聚类
    km_start = time.time()
    km = MiniBatchKMeans(n_clusters=n_clusters, random_state=0, batch_size=1000)
    labelskm = km.fit_predict(values)
    labelskm = [f'CN{i}' for i in labelskm]
    adata.obs["CN"] = pd.Categorical(pd.Series(labelskm, index=adata.obs.index))
    k_centroids = km.cluster_centers_
    print(f"KMeans聚类完成, 耗时: {time.time()-km_start:.2f}秒")
    
    # 7. 计算fold change
    fc_start = time.time()
    tissue_avgs = values.mean(axis=0)
    fc = np.log2(((k_centroids+tissue_avgs)/(k_centroids+tissue_avgs).sum(axis = 1, keepdims = True))/tissue_avgs)
    fc_df = pd.DataFrame(
        fc, 
        columns=feature_cols, 
        index=[f"CN{i}" for i in range(n_clusters)]
    )
    fc_df.to_csv(os.path.join(out_dir, 'cn_heatmap_data.xls'), sep='\t')
    print(f"Fold change计算完成, 耗时: {time.time()-fc_start:.2f}秒")
    
    # 8. 绘制热图
    plot_start = time.time()
    plt.figure(figsize=(max(8, len(feature_cols)//2), max(6, n_clusters//3)))
    g = sns.clustermap(
        fc_df, 
        vmin=-2, 
        vmax=2, 
        cmap="vlag", 
        z_score=1,  # 对列进行标准化
        row_cluster=False, 
        col_cluster=True, 
        linewidths=0.5, 
        figsize=(6, 5)
    )
    g.ax_heatmap.tick_params(right=False, bottom=False)
    plt.savefig(os.path.join(out_dir, f"CNs_{n_clusters}.png"), dpi=300, bbox_inches='tight', facecolor='white')
    plt.savefig(os.path.join(out_dir, f"CNs_{n_clusters}.pdf"), dpi=300, bbox_inches='tight', facecolor='white')
    print(f"热图绘制完成, 耗时: {time.time()-plot_start:.2f}秒")
    
    # 9. 计算细胞邻域比例
    prop_start = time.time()
    cluster_counts = adata.obs.groupby([sample, 'CN']).size().unstack(fill_value=0)
    cluster_proportions = cluster_counts.div(cluster_counts.sum(axis=1), axis=0)
    cluster_proportions_long = cluster_proportions.reset_index().melt(
        id_vars=sample, var_name='Cluster', value_name='Proportion')
    
    # 添加分组信息（如果存在）
    if 'group' in adata.obs.columns:
        cluster_proportions_long = cluster_proportions_long.merge(
            adata.obs[[sample, 'group']].drop_duplicates(), on=sample)
    else:
        cluster_proportions_long['group'] = cluster_proportions_long['sample']
    cluster_proportions_long.to_csv(os.path.join(out_dir, 'cn_proportion.txt'), sep='\t', index=None)
    print(f"邻域比例计算完成, 耗时: {time.time()-prop_start:.2f}秒")
    
    total_time = time.time() - start_time
    print(f"分析完成! 总耗时: {total_time:.2f}秒")
    
    return adata

CN celltype enrichment heatmap

CN in spatial

2. CN细胞通讯

LIANA+ (LIANA+ provides an all-in-one framework for cell–cell communication inference)一个能够超越单一任务或数据类型进行细胞间通讯推断（CCC）的集成框架，重新实现并改进了八种配体-受体方法，用于从单细胞数据中推断相互作用，并提供了一个灵活的共识机制，可以整合这些方法的任意组合，适用于单细胞和空间分辨的多组学数据。还可以通过配体-受体推断弱反映共定位.

cpdb_res_sample = {}
for sample in adata.obs['sample'].drop_duplicates():
    tmp = adata[adata.obs['sample']==sample, :]
    cellphonedb(
        tmp,
        groupby='cell_type_clean',
        resource_name='mouseconsensus',
        # resource=resource,
        expr_prop=0.1,
        key_added='cpdb_res', 
        use_raw=True,
        verbose=True, 
    )
    cpdb_res_sample[sample] = tmp.uns['cpdb_res']
tmp = []
for k,v in cpdb_res_sample.items():
    v['sample'] = k
    v['source_target_celltype'] = v['source'] + "|" + v['target']
    v['L_R'] = v['ligand'] + "->" + v['receptor']
    tmp.append(v)
pd.concat(tmp).to_csv('cpdb_results.xls', sep='\t')
adata.uns['cpdb_res'] = pd.concat(tmp)

3. 邻近富集分析（Proximity Enrichment Analysis）

空间转录组数据分析中的一种重要方法，用于评估细胞类型对之间的空间关联性。该分析通过置换检验（permutation-based tests）来比较实际观察到的细胞类型对出现的频率与随机置换数据中的频率，从而确定细胞类型之间的空间关系是否具有统计学意义。分析结果以富集z分数（enrichment z-scores）的形式输出，这些分数反映了细胞类型对之间实际的邻近频率与随机预期频率的差异程度。详细介绍见：空间转录组数据下游分析(二)

def calculate_enrichment(
    adata: "AnnData", 
    groupby: str = "celltype", 
    n_permutations: int = 100, 
    niche_radius: float = 15.0,
    permute_radius: float = 50.0, 
    spatial_key: str = "spatial",
    seed: int = 123
) -> pd.DataFrame:
    
    categories_str_cat = list(adata.obs[groupby].cat.categories)
    categories_num_cat = range(len(categories_str_cat))
    map_dict = dict(zip(categories_num_cat, categories_str_cat))
    categories_str = adata.obs[groupby]
    categories_num = adata.obs[groupby].replace(categories_str_cat, categories_num_cat)
    labels = categories_num.to_numpy()
    print(f"Total number of {len(labels)} cells")
    
    max_index = len(categories_num_cat)
    pair_counts = np.zeros((max_index, max_index))
    
    # Calculate the true interactions
    con = radius_neighbors_graph(adata.obsm[spatial_key], niche_radius,  mode='connectivity', include_self=False)
    con_coo = coo_matrix(con)
    print(f"Calculating observed interactions...")
    for i, j, val in zip(con_coo.row, con_coo.col, con_coo.data):
        if i >= j:  
            continue
        type1 = labels[i]
        type2 = labels[j]
        if val:  
            if type1 != type2:
                pair_counts[type1, type2] += 1
                pair_counts[type2, type1] += 1
            else:
                pair_counts[type1, type2] += 1
                
    coords = adata.obsm[spatial_key]
    pair_counts_null = np.zeros((max_index, max_index, n_permutations))
    for perm in range(n_permutations):
        np.random.seed(seed=None if seed is None else seed + perm)
        
        permuted_coords = coords + np.random.uniform(-permute_radius, permute_radius, size=coords.shape)
        permuted_con = radius_neighbors_graph(permuted_coords, niche_radius, mode='connectivity', include_self=False)
        permuted_con_coo = coo_matrix(permuted_con)
        
        if (perm + 1) % 100 == 1:
            print(f"Permutation iterations: {perm}...")
        pair_counts_permuted = np.zeros((max_index, max_index))
        for i, j, val in zip(permuted_con_coo.row, permuted_con_coo.col, permuted_con_coo.data):
            if i >= j:  
                continue
            type1 = labels[i]
            type2 = labels[j]
            if val:  
                if type1 != type2:
                    pair_counts_permuted[type1, type2] += 1
                    pair_counts_permuted[type2, type1] += 1
                else:
                    pair_counts_permuted[type1, type2] += 1
        
        pair_counts_null[:, :, perm] = pair_counts_permuted
    pair_counts_permuted_means = np.mean(pair_counts_null, axis=2)
    pair_counts_permuted_stds = np.std(pair_counts_null, axis=2)
    
    z_scores = (pair_counts - pair_counts_permuted_means) / pair_counts_permuted_stds
    z_scores[z_scores < 0] = 0
    z_score_df = pd.DataFrame(z_scores, index=categories_str_cat, columns=categories_str_cat)
    print("Finished!")
    
    return z_score_df