单细胞测序质控分析（QC）实战指南：从数据加载到低质量细胞识别

1. 单细胞测序质控分析的重要性第一次接触单细胞测序数据时我被那些密密麻麻的数字表格搞得头晕眼花。直到导师指着某个样本说这个细胞已经死了我才恍然大悟——原来原始数据里藏着这么多坑。单细胞测序就像给每个细胞拍X光片但有些患者在检测前就已经负伤我们需要先找出这些不合格的片子。常见的低质量细胞主要有三类症状首先是营养不良型这些细胞中检测到的基因总数library size少得可怜其次是表达缺陷型能检测到的基因数量number of expressed genes严重不足最狡猾的是线粒体泄露型它们的线粒体基因占比异常升高。就像体检时血红蛋白、白细胞、血小板三个指标同时异常这个细胞大概率出了问题。我曾分析过一组大脑组织数据原始3005个细胞中竟有799个线粒体基因占比超过10%。如果不做质控这些细胞会像噪音发生器一样干扰后续分析。更可怕的是有些低质量细胞会伪装成稀有细胞类型导致我们误判细胞亚群比例。去年有个合作项目就因此走了弯路花了两个月才发现在差异表达分析中显著的基因其实是技术假象。2. 数据加载与初步探索用R语言处理单细胞数据就像组装乐高积木需要先准备好各个模块。我习惯用scRNAseq包里的测试数据集上手比如经典的Zeisel脑细胞数据library(scRNAseq) example_sce - ZeiselBrainData()运行后控制台会打印这个数据集的体检报告包含20006个基因和3005个细胞主要信息存储在counts这个仓库里。这里有个新手容易忽略的细节——rowData和colData这两个档案柜前者记录基因特征后者存储细胞元信息。查看数据维度时我总会多做个检查dim(counts(example_sce)) # 确认基因×细胞矩阵 head(colData(example_sce)) # 查看细胞注释信息最近帮同事调试代码时发现有些数据集会把线粒体基因标记为MT-开头而非mt-这种大小写差异会导致后续计算错误。建议先用以下命令确认grep(^mt-, rownames(example_sce), ignore.caseTRUE, valueTRUE)3. 核心质控指标详解3.1 文库大小与基因数量计算质控指标就像给细胞做体检scater包的perCellQCMetrics是现成的体检仪器library(scater) is.mt - grep(^mt-, rownames(example_sce)) df - perCellQCMetrics(example_sce, subsetslist(Mitois.mt))得到的df数据框包含12项指标其中最重要的是sum细胞总UMI数相当于细胞RNA总量detected表达量超过检测阈值的基因数subsets_Mito_percent线粒体基因占比有个经验公式健康哺乳动物细胞的线粒体基因占比通常在5-15%之间。但去年分析斑马鱼数据时我们发现正常值范围是10-25%这说明阈值需要根据物种调整。3.2 线粒体基因的警示作用线粒体基因比例升高就像细胞的发烧症状可能意味着细胞膜破损导致胞质RNA流失细胞处于凋亡早期阶段样本处理时机械损伤但要注意特殊情况心肌细胞天生线粒体丰富造血干细胞则比例较低。我常用双指标联合作战plot(df$sum, df$subsets_Mito_percent, xlabLibrary size, ylabMT%) abline(h10, colred) # 线粒体阈值 abline(vquantile(df$sum, 0.1), colblue) # 文库大小10%分位数4. 低质量细胞识别实战4.1 固定阈值法固定阈值就像统一录取分数线qc.lib - df$sum 100000 qc.nexprs - df$detected 5000 qc.mito - df$subsets_Mito_percent 10 discard - qc.lib | qc.nexprs | qc.mito但这种方法有个致命缺陷——容易误伤稀有细胞类型。有次分析肿瘤微环境数据固定阈值竟然过滤掉了80%的循环肿瘤细胞就是因为这类细胞天然RNA含量低。4.2 自适应阈值法scater包的isOutlier函数更智能它会根据数据分布动态划界qc.lib2 - isOutlier(df$sum, logTRUE, typelower) qc.mito2 - isOutlier(df$subsets_Mito_percent, typehigher) discard2 - qc.lib2 | qc.mito2这个方法的本质是寻找落单的离群值默认以中位数±3MAD为界。不过要注意两个前提大部分细胞是高质量的技术变异大于生物变异4.3 批次敏感型阈值当数据存在批次效应时需要分批次计算阈值batch - example_sce$tissue # 假设tissue列记录批次 qc.batch - isOutlier(df$sum, batchbatch, logTRUE)曾处理过跨实验室整合数据发现某个批次的线粒体基因阈值竟是其他批次的2倍。后来追查发现是该实验室使用了不同的消化酶。5. 质控结果可视化好的可视化能揭示隐藏问题。我必做的三张图图1指标分布直方图hist(df$sum, breaks100, mainLibrary size distribution) abline(vattr(qc.lib2, thresholds), colred)图2指标间散点图plotColData(example_sce, xsum, ysubsets_Mito_percent, colour_bydiscard)图3过滤细胞统计upsetR::upset(data.frame(qc.lib2, qc.mito2))最近发现scater的plotHighestExprs函数也很有用它能显示表达量最高的基因帮助发现异常高表达的污染基因。6. 常见陷阱与解决方案陷阱1过度过滤解决方案先保留疑似低质量细胞在下游分析时观察其聚类位置陷阱2忽略批次差异解决方案用modelGeneVarByBatch检查批次效应陷阱3参数僵化解决方案准备阳性对照已知细胞类型测试过滤效果上周有个典型案例用户设置线粒体基因阈值5%结果过滤后完全找不到星形胶质细胞。调整到8%后细胞类型比例与文献报道一致。这提醒我们阈值需要结合生物学知识调整。7. 进阶技巧对于复杂数据集我会采用组合策略先使用自适应阈值初筛人工检查离群细胞在t-SNE图中的位置对临界细胞进行差异表达分析# 保留临界细胞示例 uncertain - df$subsets_Mito_percent 8 df$subsets_Mito_percent 12 example_sce$uncertain - uncertain存储质控结果时我习惯用addPerCellQC直接写入对象example_sce - addPerCellQC(example_sce, subsetslist(Mitois.mt))这样在后续分析中可以随时调用这些指标。比如做差异表达时排除低质量细胞的影响de_results - findMarkers(example_sce[,!discard2], ...)记住质控不是一锤子买卖。在聚类、轨迹分析等步骤后都应该回头检查低质量细胞的分布情况。有次做发育轨迹分析发现某个分支全是高线粒体比例的细胞后来证实这些是应激细胞而非真正的发育状态。