1. 单细胞分析为什么选择Scanpy?
单细胞RNA测序技术近年来快速发展,已经成为生命科学研究的重要工具。相比传统的批量测序,单细胞测序能够揭示细胞异质性,发现新的细胞类型。在分析工具方面,Python生态中的Scanpy凭借其高效、灵活的特性,正在成为越来越多研究者的首选。
我刚开始接触单细胞分析时也用过R语言的Seurat,后来转用Scanpy后发现确实更方便。Scanpy基于Python的Anndata数据结构,整合了Numpy、Scipy、Pandas等科学计算库,处理大数据集时速度优势明显。更重要的是,它的API设计非常直观,从数据预处理到可视化的完整流程都能用简洁的代码实现。
Scanpy的核心优势在于:
- 内存效率高:使用稀疏矩阵存储表达数据,节省内存
- 计算速度快:底层基于C优化的数值计算库
- 可视化丰富:内置多种专业绘图函数
- 扩展性强:可以方便地与其他Python机器学习库集成
提示:如果你已经熟悉Seurat,学习Scanpy会很容易,两者的分析流程非常相似,只是实现语言不同。
2. 环境准备与数据导入
2.1 安装与配置
建议使用conda创建独立的Python环境:
conda create -n sc_analysis python=3.8 conda activate sc_analysis pip install scanpy umap-learn leidenalg我推荐使用Jupyter Notebook进行交互式分析,方便实时查看结果。配置绘图参数可以放在开头:
import scanpy as sc sc.settings.verbosity = 3 # 设置日志级别 sc.settings.set_figure_params(dpi=100, facecolor='white') # 设置图片分辨率2.2 数据加载实战
Scanpy支持多种数据格式,最常见的是10X Genomics的标准输出。假设数据存放在filtered_gene_bc_matrices/hg19/目录:
adata = sc.read_10x_mtx( './filtered_gene_bc_matrices/hg19/', var_names='gene_symbols', # 使用基因符号作为变量名 cache=True) # 缓存加速后续读取 adata.var_names_make_unique() # 确保基因名唯一如果是其他格式的表达矩阵,可以手动构建Anndata对象:
import pandas as pd df = pd.read_csv('expression_matrix.csv', index_col=0) adata = sc.AnnData(df)3. 数据预处理与质量控制
3.1 基础过滤
单细胞数据通常包含大量噪声,需要进行严格的质量控制:
# 过滤低质量细胞和基因 sc.pp.filter_cells(adata, min_genes=200) # 每个细胞至少表达200个基因 sc.pp.filter_genes(adata, min_cells=3) # 每个基因至少在3个细胞中表达3.2 线粒体基因过滤
线粒体基因占比是判断细胞质量的重要指标:
# 标记线粒体基因 adata.var['mt'] = adata.var_names.str.startswith('MT-') # 计算QC指标 sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, inplace=True) # 可视化QC指标 sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4, multi_panel=True)3.3 数据标准化
标准化是后续分析的关键步骤:
# 文库大小归一化 sc.pp.normalize_total(adata, target_sum=1e4) # 对数转换 sc.pp.log1p(adata)4. 特征选择与降维
4.1 高变基因筛选
识别高度变化的基因可以提高信噪比:
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5) sc.pl.highly_variable_genes(adata)4.2 PCA降维
主成分分析是最常用的线性降维方法:
# 使用高变基因 adata = adata[:, adata.var.highly_variable] # 去除技术噪声 sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt']) # 标准化 sc.pp.scale(adata, max_value=10) # PCA分析 sc.tl.pca(adata, svd_solver='arpack') # 可视化方差解释率 sc.pl.pca_variance_ratio(adata, log=True)5. 细胞聚类与可视化
5.1 邻域图构建
基于PCA结果构建细胞间的邻域关系:
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)5.2 UMAP/t-SNE可视化
非线性降维方法能更好展示细胞亚群:
sc.tl.umap(adata) sc.pl.umap(adata, color=['CST3', 'NKG7'])5.3 Leiden聚类
比传统Louvain算法更优的社区发现算法:
sc.tl.leiden(adata) sc.pl.umap(adata, color=['leiden'])6. 差异表达分析与注释
6.1 标记基因识别
使用Wilcoxon秩和检验找各群的标记基因:
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon') sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)6.2 细胞类型注释
根据标记基因手动注释细胞类型:
new_cluster_names = [ 'CD4 T', 'CD14 Mono', 'B', 'CD8 T', 'NK', 'FCGR3A Mono', 'DC', 'Platelet'] adata.rename_categories('leiden', new_cluster_names)7. 高级分析与技巧
7.1 轨迹推断
使用PAGA分析细胞发育轨迹:
sc.tl.paga(adata) sc.pl.paga(adata, plot=False) # 基础布局 sc.tl.umap(adata, init_pos='paga')7.2 细胞周期分析
评估细胞周期阶段的影响:
# 加载周期基因 cell_cycle_genes = [x.strip() for x in open('cell_cycle_genes.txt')] # 计算周期分数 sc.tl.score_genes_cell_cycle(adata, cell_cycle_genes)在实际项目中,我发现Scanpy的灵活性和Python生态的结合让它特别适合复杂分析流程的开发。比如可以很方便地整合scVelo进行RNA速率分析,或者用CellPhoneDB进行细胞间互作研究。