Scanpy 1.10 单细胞质控实战3指标联合过滤与Scrublet双细胞识别附Python代码1. 单细胞RNA测序质控的核心逻辑单细胞RNA测序(scRNA-seq)数据的质量控制(QC)是分析流程中至关重要的第一步。与bulk RNA-seq不同单细胞数据具有更高的技术噪音和生物学异质性这使得质控环节成为确保下游分析可靠性的关键屏障。在单细胞数据中常见的低质量细胞来源包括破损细胞细胞膜破裂导致RNA大量流失凋亡细胞正在经历程序性死亡的细胞双细胞/多细胞(doublets/multiplets)一个液滴中包含多个细胞空液滴未捕获细胞的液滴这些低质量细胞会引入系统性偏差影响后续的细胞分群、差异表达分析等结果。因此建立科学的质控标准对单细胞分析至关重要。2. 核心质控指标与阈值设定2.1 三大核心质控指标在Scanpy中我们主要依赖三个核心指标进行质控指标类型计算方式生物学意义异常表现每个细胞的基因数n_genes_by_counts细胞转录组复杂度过低破损细胞过高双细胞每个细胞的UMI总数total_counts测序深度过低低质量细胞过高双细胞线粒体基因占比pct_counts_mt细胞活性过高破损或凋亡细胞线粒体基因的高表达通常意味着细胞应激或凋亡因为当细胞膜破损时胞质RNA会流失而线粒体RNA由于位置固定得以保留。2.2 基于MAD的自动阈值设定传统质控通常使用固定阈值但这种方法无法适应不同数据集的特异性。Scanpy 1.10推荐使用中值绝对偏差(MAD)方法自动计算阈值# 计算MAD自动阈值 def mad_threshold(series, n_mads3): median series.median() mad (series - median).abs().median() * 1.4826 return median n_mads * mad # 应用自动阈值 min_genes mad_threshold(adata.obs[n_genes_by_counts], n_mads-3) max_genes mad_threshold(adata.obs[n_genes_by_counts], n_mads3) min_counts mad_threshold(adata.obs[total_counts], n_mads-3) max_counts mad_threshold(adata.obs[total_counts], n_mads3) max_mt mad_threshold(adata.obs[pct_counts_mt], n_mads3)这种方法能根据数据分布自动调整阈值特别适合处理不同来源或不同实验条件的数据集。3. 质控流程实战演示3.1 数据加载与初步质控首先加载数据并计算基础指标import scanpy as sc import numpy as np # 读取数据 adata sc.read_10x_mtx(filtered_feature_bc_matrix/, var_namesgene_symbols) # 计算基础QC指标 adata.var[mt] adata.var_names.str.startswith(MT-) # 标记线粒体基因 sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue) # 可视化QC指标分布 sc.pl.violin(adata, [n_genes_by_counts, total_counts, pct_counts_mt], jitter0.4, multi_panelTrue)3.2 三指标联合过滤基于计算得到的阈值进行过滤# 设置过滤条件 min_genes 500 # 最少检测到的基因数 max_genes 6000 # 最多检测到的基因数 max_mt 15 # 最大线粒体基因占比 # 应用过滤 sc.pp.filter_cells(adata, min_genesmin_genes) sc.pp.filter_cells(adata, max_genesmax_genes) adata adata[adata.obs[pct_counts_mt] max_mt, :] # 记录过滤后数据 print(f过滤后细胞数: {adata.n_obs})3.3 Scrublet双细胞识别双细胞(doublets)是单细胞测序中的常见问题Scrublet通过模拟双细胞表达谱来识别真实数据中的双细胞# 安装Scrublet !pip install scrublet # 运行Scrublet import scrublet as scr scrub scr.Scrublet(adata.X) doublet_scores, predicted_doublets scrub.scrub_doublets() # 将结果添加到adata对象 adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets # 可视化双细胞得分 sc.pl.scatter(adata, xtotal_counts, ydoublet_score, colorpredicted_doublet) # 过滤双细胞 adata adata[~adata.obs[predicted_doublet], :].copy()4. 质控后数据标准化与可视化完成质控后我们需要对数据进行标准化处理# 基本预处理流程 sc.pp.normalize_total(adata, target_sum1e4) # CPM标准化 sc.pp.log1p(adata) # 对数转换 sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5) # 筛选高变基因 adata adata[:, adata.var.highly_variable] # 保留高变基因 # PCA降维 sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack) # UMAP可视化 sc.pp.neighbors(adata, n_pcs30, n_neighbors20) sc.tl.umap(adata) sc.pl.umap(adata, color[n_genes_by_counts, total_counts, pct_counts_mt])5. 质控结果解读与注意事项5.1 质控决策要点阈值设置的灵活性不同组织类型、不同实验protocol可能需要调整阈值指标间的相关性高UMI数常伴随高基因数但线粒体占比应独立评估批次效应影响跨批次数据应分别质控后再整合5.2 常见问题解决方案问题现象可能原因解决方案过滤后细胞数过少阈值设置过严放宽MAD系数或手动调整阈值双细胞比例异常高细胞悬液浓度过高检查实验记录考虑重新分析线粒体基因普遍偏高样本处理不当检查样本采集到建库的时间间隔6. 完整代码实现以下是整合所有步骤的完整代码# 单细胞质控全流程 import scanpy as sc import numpy as np import scrublet as scr # 1. 数据加载 adata sc.read_10x_mtx(data/, var_namesgene_symbols) # 2. 计算QC指标 adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue) # 3. 自动阈值计算 def mad_threshold(series, n_mads3): median series.median() mad (series - median).abs().median() * 1.4826 return median n_mads * mad min_genes mad_threshold(adata.obs[n_genes_by_counts], n_mads-3) max_genes mad_threshold(adata.obs[n_genes_by_counts], n_mads3) min_counts mad_threshold(adata.obs[total_counts], n_mads-3) max_counts mad_threshold(adata.obs[total_counts], n_mads3) max_mt 15 # 线粒体阈值通常固定 # 4. 细胞过滤 sc.pp.filter_cells(adata, min_genesmin_genes) sc.pp.filter_cells(adata, max_genesmax_genes) adata adata[adata.obs[pct_counts_mt] max_mt, :] # 5. 双细胞检测 scrub scr.Scrublet(adata.X) doublet_scores, predicted_doublets scrub.scrub_doublets() adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets adata adata[~adata.obs[predicted_doublet], :].copy() # 6. 数据标准化与可视化 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5) adata adata[:, adata.var.highly_variable] sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack) sc.pp.neighbors(adata, n_pcs30, n_neighbors20) sc.tl.umap(adata) # 7. 结果可视化 sc.pl.violin(adata, [n_genes_by_counts, total_counts, pct_counts_mt], jitter0.4) sc.pl.scatter(adata, xtotal_counts, ypct_counts_mt, colorn_genes_by_counts) sc.pl.umap(adata, color[n_genes_by_counts, total_counts, pct_counts_mt])在实际项目中我们还需要考虑样本特异性、实验批次等因素可能需要针对不同样本单独进行质控后再整合分析。Scanpy提供的这套质控流程已经过大量实际项目验证能够为后续分析提供可靠的数据基础。
