从数据到洞见：用Scanpy搞定单细胞测序分析的完整实战流程（附代码）-港品优选

从数据到洞见：用Scanpy搞定单细胞测序分析的完整实战流程（附代码）

单细胞测序技术正在彻底改变我们对复杂生物系统的理解能力。想象一下，你手中握有一份来自10x Genomics平台的原始数据，里面蕴含着数千个细胞的基因表达信息。如何将这些看似杂乱的数据转化为有意义的生物学发现？这就是Scanpy——基于Python的单细胞分析工具链要解决的问题。

对于刚接触这一领域的研究者来说，最大的挑战往往不是写代码，而是理解整个分析流程的逻辑链条：为什么要做质控？如何判断过滤阈值？聚类结果怎么解释？本文将用一个真实的PBMC（外周血单个核细胞）数据集，带你走完从原始数据到可发表结果的完整流程，重点解决那些教程里很少提及的"为什么"和"怎么判断"的问题。

1. 环境准备与数据加载

1.1 工具链配置

单细胞分析需要一套完整的Python生态工具。推荐使用conda创建独立环境以避免依赖冲突：

conda create -n sc_analysis python=3.8 conda activate sc_analysis pip install scanpy leidenalg

关键组件说明：

Scanpy：核心分析框架
leidenalg：社区检测算法实现
umap-learn：降维可视化（自动依赖安装）

1.2 数据加载技巧

从10x Genomics的Cell Ranger输出目录加载数据时，注意var_names参数的选择：

import scanpy as sc adata = sc.read_10x_mtx( 'filtered_gene_bc_matrices/hg19/', var_names='gene_symbols', # 使用基因符号而非ENSEMBL ID cache=True # 启用缓存加速后续加载 )

注意：当处理多个样本时，建议使用sc.read_10x_h5直接读取h5文件，并通过concatenate合并

加载后立即检查数据的基本结构：

print(adata)

输出示例：

AnnData object with n_obs × n_vars = 2638 × 1838 obs: 'n_genes', 'percent_mito' var: 'gene_ids', 'feature_types', 'genome', 'mt'

2. 数据质控与过滤策略

2.1 质控指标可视化

质控是单细胞分析的基石，需要同时考虑三个核心指标：

每个细胞的基因检出数（n_genes_by_counts）
每个细胞的总UMI数（total_counts）
线粒体基因占比（pct_counts_mt）

使用组合图表展示分布：

sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4, multi_panel=True)

2.2 动态阈值确定

过滤阈值不应固定，而应根据数据分布动态确定：

# 计算动态阈值 mt_threshold = adata.obs.pct_counts_mt.quantile(0.95) gene_threshold = adata.obs.n_genes_by_counts.quantile(0.99) # 应用过滤 adata = adata[adata.obs.pct_counts_mt < mt_threshold, :] adata = adata[adata.obs.n_genes_by_counts < gene_threshold, :]

提示：保留过滤前的数据副本有助于后续回溯分析

2.3 双标图辅助决策

通过散点图观察指标间关系，识别异常细胞群体：

sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', color='n_genes_by_counts')

典型异常模式包括：

高线粒体+低基因数：濒死细胞
高UMI+高基因数：双细胞或多细胞

3. 数据标准化与特征选择

3.1 标准化策略对比

方法	公式	适用场景
CPM	$X_{ij} = \frac{X_{ij}}{\sum_j X_{ij}} \times 10^4$	初步标准化
Log1p	$X_{ij} = \log(X_{ij} + 1)$	方差稳定化
SCTransform	负二项回归	大样本数据

推荐基础流程：

sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) adata.raw = adata # 保存原始数据

3.2 高变基因筛选

通过分散度分析选择信息量最大的基因：

sc.pp.highly_variable_genes( adata, min_mean=0.0125, max_mean=3, min_disp=0.5 ) sc.pl.highly_variable_genes(adata)

筛选后保留高变基因：

adata = adata[:, adata.var.highly_variable]

4. 降维与聚类分析

4.1 PCA深度解读

运行PCA后，务必检查方差解释率：

sc.tl.pca(adata, svd_solver='arpack') sc.pl.pca_variance_ratio(adata, log=True)

选择PC数量的经验法则：

肘部法则：累计解释率曲线的拐点
随机矩阵理论：PC特征值超过随机分布预期

4.2 聚类优化技巧

UMAP参数对结果影响显著，推荐采用多分辨率分析：

# 邻域图构建 sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30) # 多分辨率聚类 for res in [0.3, 0.5, 0.7]: sc.tl.leiden(adata, resolution=res, key_added=f"leiden_{res}") # 可视化比较 sc.pl.umap(adata, color=['leiden_0.3', 'leiden_0.5', 'leiden_0.7'])

5. 标记基因分析与注释

5.1 差异表达分析

Wilcoxon检验是最稳健的选择：

sc.tl.rank_genes_groups( adata, 'leiden', method='wilcoxon', pts=True # 计算表达比例 )

结果解读要点：

logFC：表达量变化倍数（对数）
pvals：统计显著性
pts：在目标群中的表达比例

5.2 细胞类型注释

建立标记基因字典辅助注释：

marker_dict = { 'T cells': ['CD3D', 'CD3E'], 'B cells': ['MS4A1', 'CD79A'], 'Monocytes': ['CD14', 'LYZ'], 'NK cells': ['GNLY', 'NKG7'] }

使用点图验证注释：

sc.pl.dotplot(adata, marker_dict, groupby='leiden_0.5', dendrogram=True)

6. 高级分析与可视化

6.1 轨迹推断

使用PAGA分析细胞状态转换：

sc.tl.paga(adata) sc.pl.paga(adata, color=['leiden_0.5', 'CST3'], threshold=0.03)

6.2 基因共表达网络

识别基因模块：

sc.tl.dendrogram(adata, groupby='leiden_0.5') sc.tl.correlation_matrix(adata) sc.pl.correlation_matrix(adata, 'leiden_0.5')

7. 结果导出与报告生成

7.1 交互式可视化

使用Cellxgene工具导出可交互结果：

adata.write('results/pbmc_processed.cxg', compression='gzip')

7.2 自动化报告

生成HTML分析报告：

sc.settings.report_images = True sc.report(adata, 'results/report.html', title='PBMC单细胞分析报告')

在实际项目中，我发现最耗时的往往不是分析本身，而是参数调整和结果解释。例如，当UMAP图中出现"拉面状"结构时，可能需要调整n_neighbors参数；而聚类分辨率的选择应该基于已知的生物学预期——外周血PBMC通常会有6-10个主要细胞类型。

企业官网建设流程全解析