作者，Evil Genius
基础分析真的希望引起大家的重视，打好基础再上项目，将来进入公司工作，如果岗位是生物信息，一开始的培训就是linux、R、python，培训一周的时间就考核，考核不过，转岗或者离开；即使是课题组，数据分析也离不开服务器，离不开R、python代码语言，番外培训2025番外--linux、R、python培训
https://mp.weixin.qq.com/s/8E1vYJMNhe5m0DXieBHfzA?payreadticket=HM3XzM3dAiSdVyi83oiJ3UWd-zPIOKO1EMMnYwy9ASu_scWj_WEje8Xwiwq_z8JvgcLDfd0&scene=1&click_id=1有条件一定要从头到尾梳理一下基础，基本上好好过一遍，基础就不成问题了。
今天我们更新脚本，高精度空转（Xenium、CosMx）空间距离分析，我们以量化fibroblast与tumor的空间距离关系为例。
分析一般会更进一步，依据fibroblast与tumor的距离关系划分为近距离、中距离、远距离三组，分析其中的基因表达变化,当然这种变化也代表了距离和基因的关系。
各种情况我们都要考虑到
分析思路很简单，fibroblast与tumor细胞二维坐标的距离计算，其中以fibroblast与最近的tumor细胞的距离为实际距离。
简单问大家一个问题，为什么写python代码的时候，有经验的人员为了实现一定的生物学功能通常会采用定义函数的方式么？
当然，实现起来需要很强的python基础
#!/usr/bin/env python3
####zhaoyunfei
####20250716
# -*- coding: utf-8 -*-
"""
Xenium/CosMx空间转录组数据分析脚本
分析两种细胞类型的空间距离关系及基因表达变化
"""import scanpy as sc
import squidpy as sq
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import warnings
warnings.filterwarnings('ignore')# 设置绘图风格
plt.rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = Falseclass SpatialDistanceAnalysis:def __init__(self, adata, celltype1='tumor', celltype2='fibroblast'):"""初始化分析类Parameters:-----------adata : AnnData包含空间转录组数据的AnnData对象celltype1 : str第一种细胞类型的标签（通常为参考细胞，如tumor）celltype2 : str第二种细胞类型的标签（通常为分析细胞，如fibroblast）"""self.adata = adataself.celltype1 = celltype1self.celltype2 = celltype2self.results = {}def check_data_requirements(self):"""检查数据是否满足分析要求"""print("=" * 60)print("检查数据要求...")print("=" * 60)# 检查必需的obs列if 'celltype' not in self.adata.obs.columns:raise ValueError("adata.obs中缺少'celltype'列")# 检查空间坐标spatial_keys = [key for key in self.adata.obsm_keys() if 'spatial' in key.lower() or 'X_spatial' in key]if not spatial_keys:raise ValueError("未找到空间坐标数据，请检查adata.obsm中的空间坐标")self.spatial_key = spatial_keys[0]print(f"使用空间坐标: {self.spatial_key}")# 检查细胞类型是否存在celltypes = self.adata.obs['celltype'].unique()if self.celltype1 not in celltypes:raise ValueError(f"未找到指定的细胞类型1: {self.celltype1}")if self.celltype2 not in celltypes:raise ValueError(f"未找到指定的细胞类型2: {self.celltype2}")print("✓ 数据要求检查通过")print(f"✓ 细胞类型1 ({self.celltype1}): {(self.adata.obs['celltype'] == self.celltype1).sum()} 个细胞")print(f"✓ 细胞类型2 ({self.celltype2}): {(self.adata.obs['celltype'] == self.celltype2).sum()} 个细胞")def calculate_spatial_distances(self):"""使用squidpy计算细胞类型2到细胞类型1的空间距离"""print("\n" + "=" * 60)print("使用squidpy计算空间距离...")print("=" * 60)# 获取细胞类型信息celltype1_mask = self.adata.obs['celltype'] == self.celltype1celltype2_mask = self.adata.obs['celltype'] == self.celltype2print(f"{self.celltype1}细胞数量: {celltype1_mask.sum()}")print(f"{self.celltype2}细胞数量: {celltype2_mask.sum()}")# 使用squidpy计算最近邻距离# 首先需要确保空间坐标在正确的格式中spatial_coords = self.adata.obsm[self.spatial_key]# 计算所有细胞到celltype1细胞的距离print("计算空间距离矩阵...")sq.gr.spatial_neighbors(self.adata, coord_type='generic', spatial_key=self.spatial_key)# 提取celltype1细胞的坐标celltype1_coords = spatial_coords[celltype1_mask]celltype2_coords = spatial_coords[celltype2_mask]# 使用squidpy的距离计算# 计算每个celltype2细胞到最近celltype1细胞的距离distances = []for i in range(len(celltype2_coords)):# 计算到所有celltype1细胞的欧几里得距离dists = np.sqrt(np.sum((celltype2_coords[i] - celltype1_coords) ** 2, axis=1))min_dist = np.min(dists)distances.append(min_dist)distances = np.array(distances)# 将距离添加到adata中self.adata.obs['distance_to_celltype1'] = np.nanself.adata.obs.loc[celltype2_mask, 'distance_to_celltype1'] = distances# 保存结果self.results['distances'] = distancesself.results['celltype1_coords'] = celltype1_coordsself.results['celltype2_coords'] = celltype2_coords# 显示距离统计print(f"\n{self.celltype2}到最近{self.celltype1}的距离统计:")print(f"最小距离: {distances.min():.2f}")print(f"最大距离: {distances.max():.2f}")print(f"平均距离: {distances.mean():.2f}")print(f"中位数距离: {np.median(distances):.2f}")print(f"标准差: {distances.std():.2f}")return distancesdef categorize_distances(self, method='quantile', custom_thresholds=None):"""对距离进行分层Parameters:-----------method : str分层方法: 'quantile' (分位数) 或 'custom' (自定义阈值)custom_thresholds : list自定义阈值列表，如 [10, 20] 表示 <10, 10-20, >20"""print("\n" + "=" * 60)print("距离分层...")print("=" * 60)celltype2_mask = self.adata.obs['celltype'] == self.celltype2celltype2_distances = self.adata.obs.loc[celltype2_mask, 'distance_to_celltype1'].valuesif method == 'quantile':# 基于分位数进行距离分层distance_quantiles = np.percentile(celltype2_distances, [33, 66])self.distance_thresholds = distance_quantilesprint(f"使用分位数分层:")print(f" 近距离: 0 - {distance_quantiles[0]:.2f}")print(f" 中距离: {distance_quantiles[0]:.2f} - {distance_quantiles[1]:.2f}")print(f" 远距离: {distance_quantiles[1]:.2f} - {celltype2_distances.max():.2f}")elif method == 'custom' and custom_thresholds is not None:# 使用自定义阈值self.distance_thresholds = custom_thresholdsprint(f"使用自定义阈值分层:")print(f" 近距离: 0 - {custom_thresholds[0]:.2f}")print(f" 中距离: {custom_thresholds[0]:.2f} - {custom_thresholds[1]:.2f}")print(f" 远距离: {custom_thresholds[1]:.2f} - {celltype2_distances.max():.2f}")else:raise ValueError("请提供有效的分层方法或自定义阈值")# 定义距离分层函数def categorize_distance(distance):if distance <= self.distance_thresholds[0]:return '近距离'elif distance <= self.distance_thresholds[1]:return '中距离'else:return '远距离'# 添加距离分层信息self.adata.obs['distance_category'] = '其他细胞'distance_categories = [categorize_distance(d) for d in celltype2_distances]self.adata.obs.loc[celltype2_mask, 'distance_category'] = distance_categories# 统计各层细胞数量distance_counts = self.adata.obs[self.adata.obs['celltype'] == self.celltype2]['distance_category'].value_counts()print(f"\n各距离层{self.celltype2}细胞数量:")for category in ['近距离', '中距离', '远距离']:count = distance_counts.get(category, 0)percentage = count / len(celltype2_distances) * 100print(f" {category}: {count}个细胞 ({percentage:.1f}%)")self.results['distance_categories'] = distance_categoriesself.results['distance_counts'] = distance_countsreturn distance_categoriesdef visualize_spatial_distribution(self, save_path=None):"""使用squidpy和scanpy可视化空间分布"""print("\n" + "=" * 60)print("可视化空间分布...")print("=" * 60)# 创建图形fig, axes = plt.subplots(2, 2, figsize=(16, 14))fig.suptitle(f'{self.celltype2}与{self.celltype1}的空间距离分析', fontsize=16, fontweight='bold')# 子图1: 使用scanpy绘制所有细胞类型ax1 = axes[0, 0]sc.pl.spatial(self.adata, color='celltype', ax=ax1, show=False, title='所有细胞类型空间分布', frameon=False)# 子图2: 突出显示两种目标细胞类型ax2 = axes[0, 1]# 创建临时列用于突出显示self.adata.obs['highlight'] = '其他细胞'self.adata.obs.loc[self.adata.obs['celltype'] == self.celltype1, 'highlight'] = self.celltype1self.adata.obs.loc[self.adata.obs['celltype'] == self.celltype2, 'highlight'] = self.celltype2sc.pl.spatial(self.adata, color='highlight', ax=ax2, show=False,title=f'{self.celltype1}(红) vs {self.celltype2}(蓝)', palette={self.celltype1: 'red', self.celltype2: 'blue', '其他细胞': 'lightgray'},frameon=False)# 子图3: 距离分层可视化ax3 = axes[1, 0]sc.pl.spatial(self.adata, color='distance_category', ax=ax3, show=False,title=f'{self.celltype2}距离分层',palette={'近距离': 'darkred', '中距离': 'orange', '远距离': 'lightblue', '其他细胞': 'lightgray'},frameon=False)# 子图4: 距离分布直方图ax4 = axes[1, 1]celltype2_mask = self.adata.obs['celltype'] == self.celltype2celltype2_distances = self.adata.obs.loc[celltype2_mask, 'distance_to_celltype1'].values# 绘制直方图n, bins, patches = ax4.hist(celltype2_distances, bins=30, alpha=0.7, color='skyblue', edgecolor='black')# 添加分层阈值线ax4.axvline(self.distance_thresholds[0], color='red', linestyle='--', label=f'近距离阈值: {self.distance_thresholds[0]:.2f}')ax4.axvline(self.distance_thresholds[1], color='orange', linestyle='--', label=f'中距离阈值: {self.distance_thresholds[1]:.2f}')ax4.set_xlabel(f'到最近{self.celltype1}细胞的距离')ax4.set_ylabel('细胞数量')ax4.set_title(f'{self.celltype2}细胞距离分布')ax4.legend()ax4.grid(True, alpha=0.3)plt.tight_layout()if save_path:plt.savefig(save_path, dpi=300, bbox_inches='tight')print(f"✓ 空间分布图已保存至: {save_path}")plt.show()def perform_differential_expression(self, min_cells_per_group=5):"""使用scanpy执行差异表达分析"""print("\n" + "=" * 60)print("使用scanpy执行差异表达分析...")print("=" * 60)# 提取celltype2细胞子集adata_celltype2 = self.adata[self.adata.obs['celltype'] == self.celltype2].copy()# 检查各组的细胞数量group_counts = adata_celltype2.obs['distance_category'].value_counts()valid_groups = group_counts[group_counts >= min_cells_per_group].index.tolist()print(f"各组细胞数量:")for group in ['近距离', '中距离', '远距离']:count = group_counts.get(group, 0)status = "✓" if count >= min_cells_per_group else "✗"print(f" {status} {group}: {count}个细胞")if len(valid_groups) < 2:print("警告: 有效组别数量不足，无法进行差异表达分析")return None# 进行差异表达分析print(f"\n进行差异表达分析，以'近距离'为参考组...")try:# 使用scanpy的rank_genes_groups进行差异表达分析sc.tl.rank_genes_groups(adata_celltype2, 'distance_category', method='wilcoxon', groups=valid_groups,reference='近距离')# 获取差异表达结果de_results = {}for group in valid_groups:if group != '近距离': # 跳过参考组de_df = sc.get.rank_genes_groups_df(adata_celltype2, group=group)de_results[f"{group}_vs_近距离"] = de_dfself.results['de_results'] = de_resultsself.results['adata_celltype2'] = adata_celltype2# 显示显著的差异表达基因print(f"\n显著差异表达基因 (p_val_adj < 0.05, |log2FC| > 0.5):")for comparison, de_df in de_results.items():sig_genes = de_df[(de_df['pvals_adj'] < 0.05) & (abs(de_df['logfoldchanges']) > 0.5)]print(f" {comparison}: {len(sig_genes)}个基因")if len(sig_genes) > 0:print(f" 上调基因: {len(sig_genes[sig_genes['logfoldchanges'] > 0])}个")print(f" 下调基因: {len(sig_genes[sig_genes['logfoldchanges'] < 0])}个")if len(sig_genes) <= 10: # 如果基因不多，显示具体基因for _, gene in sig_genes.head(5).iterrows():direction = "↑" if gene['logfoldchanges'] > 0 else "↓"print(f" {direction} {gene['names']} (log2FC: {gene['logfoldchanges']:.2f}, p_adj: {gene['pvals_adj']:.2e})")return de_resultsexcept Exception as e:print(f"差异表达分析失败: {e}")return Nonedef visualize_expression_analysis(self, top_genes=10, save_path=None):"""使用scanpy可视化表达分析结果"""if 'de_results' not in self.results:print("未找到差异表达分析结果，请先执行差异表达分析")returnprint("\n" + "=" * 60)print("可视化表达分析结果...")print("=" * 60)adata_celltype2 = self.results['adata_celltype2']de_results = self.results['de_results']# 创建图形fig, axes = plt.subplots(2, 2, figsize=(16, 12))fig.suptitle(f'{self.celltype2}细胞距离相关的基因表达变化', fontsize=16, fontweight='bold')# 子图1: 差异表达基因点图ax1 = axes[0, 0]# 收集所有比较组中的top基因top_genes_list = []for comparison, de_df in de_results.items():sig_genes = de_df[(de_df['pvals_adj'] < 0.05) & (abs(de_df['logfoldchanges']) > 0.5)]if len(sig_genes) > 0:top_genes_list.extend(sig_genes.head(top_genes)['names'].tolist())# 去重top_genes_list = list(set(top_genes_list))if len(top_genes_list) > 0:# 绘制点图sc.pl.dotplot(adata_celltype2, top_genes_list, groupby='distance_category',ax=ax1, show=False)ax1.set_title(f'Top {top_genes}差异表达基因')else:ax1.text(0.5, 0.5, '无显著差异表达基因', ha='center', va='center', transform=ax1.transAxes)ax1.set_title('差异表达基因')# 子图2: 基因表达violin图ax2 = axes[0, 1]if len(top_genes_list) >= 3:selected_genes = top_genes_list[:3]sc.pl.violin(adata_celltype2, selected_genes, groupby='distance_category',ax=ax2, rotation=45, show=False)ax2.set_title('代表性基因表达分布')else:ax2.text(0.5, 0.5, '基因数量不足', ha='center', va='center', transform=ax2.transAxes)ax2.set_title('基因表达分布')# 子图3: 差异表达基因热图ax3 = axes[1, 0]if len(top_genes_list) > 0:# 选择top基因进行热图展示heatmap_genes = top_genes_list[:min(20, len(top_genes_list))]# 对每个距离类别的细胞进行排序adata_sorted = adata_celltype2.copy()category_order = ['近距离', '中距离', '远距离']adata_sorted.obs['distance_category'] = pd.Categorical(adata_sorted.obs['distance_category'], categories=category_order, ordered=True)adata_sorted = adata_sorted[adata_sorted.obs['distance_category'].argsort()]sc.pl.heatmap(adata_sorted, heatmap_genes, groupby='distance_category',ax=ax3, show=False, cmap='viridis')ax3.set_title('差异表达基因热图')else:ax3.text(0.5, 0.5, '无显著差异表达基因', ha='center', va='center', transform=ax3.transAxes)ax3.set_title('差异表达基因热图')# 子图4: 各距离层表达特征总结ax4 = axes[1, 1]# 统计各比较组的差异基因数量comparisons = []up_genes = []down_genes = []for comparison, de_df in de_results.items():sig_genes = de_df[(de_df['pvals_adj'] < 0.05) & (abs(de_df['logfoldchanges']) > 0.5)]up_count = len(sig_genes[sig_genes['logfoldchanges'] > 0])down_count = len(sig_genes[sig_genes['logfoldchanges'] < 0])comparisons.append(comparison.replace('_vs_近距离', ''))up_genes.append(up_count)down_genes.append(down_count)if comparisons:x = np.arange(len(comparisons))width = 0.35ax4.bar(x - width/2, up_genes, width, label='上调基因', color='red', alpha=0.7)ax4.bar(x + width/2, down_genes, width, label='下调基因', color='blue', alpha=0.7)ax4.set_xlabel('比较组')ax4.set_ylabel('基因数量')ax4.set_title('各比较组差异表达基因数量')ax4.set_xticks(x)ax4.set_xticklabels(comparisons)ax4.legend()ax4.grid(True, alpha=0.3)else:ax4.text(0.5, 0.5, '无显著差异表达基因', ha='center', va='center', transform=ax4.transAxes)ax4.set_title('差异表达基因统计')plt.tight_layout()if save_path:plt.savefig(save_path, dpi=300, bbox_inches='tight')print(f"✓ 表达分析图已保存至: {save_path}")plt.show()def spatial_ligand_receptor_analysis(self):"""使用squidpy进行配体-受体相互作用分析"""print("\n" + "=" * 60)print("配体-受体相互作用分析...")print("=" * 60)try:# 安装squidpy的配体-受体数据库sq.gr.ligrec(self.adata, cluster_key='celltype', groups=[self.celltype1, self.celltype2])# 获取配体-受体相互作用结果lr_results = self.adata.uns['ligrec']['pvalues']print("配体-受体分析完成")# 可视化top相互作用fig, ax = plt.subplots(figsize=(10, 8))sq.pl.ligrec(self.adata, cluster_key='celltype', source_groups=[self.celltype1], target_groups=[self.celltype2],ax=ax)ax.set_title(f'{self.celltype1} -> {self.celltype2} 配体-受体相互作用')plt.tight_layout()plt.show()return lr_resultsexcept Exception as e:print(f"配体-受体分析失败: {e}")print("请确保已安装正确的squidpy版本和数据库")return Nonedef generate_report(self):"""生成分析报告"""print("\n" + "=" * 60)print("分析报告")print("=" * 60)celltype1_count = (self.adata.obs['celltype'] == self.celltype1).sum()celltype2_count = (self.adata.obs['celltype'] == self.celltype2).sum()print(f"分析总结:")print(f" - 细胞类型1 ({self.celltype1}): {celltype1_count} 个细胞")print(f" - 细胞类型2 ({self.celltype2}): {celltype2_count} 个细胞")if 'distances' in self.results:distances = self.results['distances']print(f" - {self.celltype2}到{self.celltype1}的平均距离: {distances.mean():.2f} ± {distances.std():.2f}")print(f" - 距离范围: {distances.min():.2f} - {distances.max():.2f}")if 'distance_counts' in self.results:distance_counts = self.results['distance_counts']print(f"\n距离分层统计:")for category in ['近距离', '中距离', '远距离']:count = distance_counts.get(category, 0)percentage = count / celltype2_count * 100print(f" - {category}: {count}个细胞 ({percentage:.1f}%)")if 'de_results' in self.results:de_results = self.results['de_results']print(f"\n差异表达分析结果:")for comparison, de_df in de_results.items():sig_genes = de_df[(de_df['pvals_adj'] < 0.05) & (abs(de_df['logfoldchanges']) > 0.5)]print(f" - {comparison}: {len(sig_genes)}个显著差异基因")print(f"\n✓ 分析完成!")# 使用示例
def main():"""主函数 - 使用示例"""# 示例数据加载 (请根据实际情况修改)# adata = sc.read('Xenium_data.h5ad')# 创建分析对象analyzer = SpatialDistanceAnalysis(adata, celltype1='tumor', celltype2='fibroblast')# 执行完整分析流程analyzer.check_data_requirements()analyzer.calculate_spatial_distances()analyzer.categorize_distances(method='quantile')analyzer.visualize_spatial_distribution(save_path='spatial_distribution.png')analyzer.perform_differential_expression(min_cells_per_group=5)analyzer.visualize_expression_analysis(save_path='expression_analysis.png')analyzer.generate_report()print("请根据您的数据修改main函数中的代码")if __name__ == "__main__":main()
生活很好，有你更好
