单细胞分析入门十个有八个把注意力都放在后面的聚类、找marker基因、做拟时序这些“高光步骤”上结果一跑起来就发现数据一团糟后面所有的图都解释不通。我最早做单细胞数据质控的时候就吃过这个亏拿到一个公共数据集看都没细看就直接往下走最后的UMAP图上细胞群体边界模糊得没法看重新回头补课才发现问题全出在源头。数据质控这件事听着基础实际上直接决定你后面所有分析靠不靠谱。这篇内容我不打算给你讲那些大而全的理论而是围绕一个完整的实操链路——从GEO下载数据、整理表达矩阵、跑质控指标、定阈值过滤到质控后的验证——把每一步为什么要这么做、怎么判断结果合不合理讲明白。适合刚接触单细胞分析、或者跑过流程但心里没底的朋友照着做能少走很多弯路。1. 为什么我把质控当作单细胞分析的“第一面墙”1.1 你的数据里混着什么质控要解决的三个真实问题单细胞测序和普通转录组测序有一个本质区别普通转录组测的是“一堆细胞平均下来”的表达量而单细胞测序拿到的是“一个一个细胞”的表达量。这个“一个一个”的实现方式决定了数据里天然会存在噪声。以10x Genomics平台的droplet-based方案为例细胞悬液和带有barcode的凝胶珠被一起包裹进微滴里每个微滴理论上只包一个细胞和一个凝胶珠。但实际运行的时候会有一批微滴包进去两个细胞doublet也有一批微滴根本没包到细胞只包到了游离的RNA也就是ambient RNA环境RNA。还有一些细胞在被捕获的时候就已经破裂死亡细胞膜不再完整细胞质里的RNA跑出去了一大部分测序仪最后捕获到的只是线粒体基因的残留片段。这三个问题对应到数据上就是三个不同维度的异常空微滴和低质量细胞会让一部分barcode对应的UMI数极低它们不配被当作一个“细胞”进入后续分析。doublet会造成一个barcode下混杂了两个不同细胞类型的表达谱聚类的时候它们可能单独形成一个奇怪的群体也可能把两个群体的边界模糊掉。破裂细胞和死亡细胞会因为mRNA降解呈现出线粒体基因比例异常偏高的特征。我见过很多人在质控这一步图省事用Seurat的默认参数跑一遍就过去了结果后面找marker基因的时候发现一堆未知基因注释细胞类型的时候怎么都对不上。实际上这些异常在质控阶段都有非常清晰的信号只是你没去看。1.2 质控不是“清垃圾”而是在定义“什么是细胞”我一直觉得质控这个翻译容易让人产生误解好像是把“坏”的东西扔掉。“坏”这个概念本身就很主观——一个线粒体基因比例20%的细胞在肿瘤组织里可能是一个真实存在的、处于应激状态的细胞你把它当成“坏死细胞”丢掉等于把真实的生物学信息也扔掉了。所以更准确的理解是质控是在定义“你的数据里什么样的barcode配叫一个可以进入分析的细胞”。你设定nGene 500意思就是“低于500个基因检测到的barcode我没有信心认为它是一个完整的细胞所以不纳入”。你设定percent.mt 20%意思就是“线粒体基因占比超过20%的barcode我倾向于认为是垂死细胞或者空微滴污染所以排除”。这个理解方式有两个实际好处。第一你不容易被网上的“标准阈值”带偏因为阈值本身是基于你的数据分布决定的而不是某个教程拍脑袋定的数字。第二当你跟审稿人或者合作者讨论的时候你能清楚地说出“我的细胞定义标准是什么”而不是含糊地说“按标准流程过滤”。我跟很多审稿人打过交道在方法学部分把质控阈值怎么定的写清楚是减少被质疑“数据质量存疑”的最有效方式。2. 从GEO下载到表达矩阵走到质控之前的路2.1 取数前先学会“读货”看懂GEO条目里的三个关键字段很多人一上来就想跑代码但GEO数据集很多时候并不会给你一个现成的“表达矩阵.csv”。你得花点时间把它的数据结构搞清楚不然下载完发现根本不是你要的东西白白浪费一个下午。看一个GEO条目比如GSE编号对应的页面我一般重点看三个字段。第一个是Overall design它告诉你这个数据集整体是怎么设计的有几组样本每组几个重复用的是哪个平台。第二个是Supplementary file这里才是真正的数据仓库——10x Genomics的数据通常会上传filtered_feature_bc_matrix、raw_feature_bc_matrix或者Cell Ranger处理好的h5文件如果是Smart-seq2这类全长转录组方案你拿到的可能是每个样本一个表达矩阵文件。第三个是SRA/ENA的原始数据入口如果你需要自己重新跑比对定量就得从这里拿FASTQ。尤其要注意的是GEO上有些数据集只有RAW data原始计数有些只提供了processed data已经过滤过的矩阵。我之前就踩过坑下载了一个已经做了过滤的矩阵自己又做了一轮过滤结果阳性细胞比例明显偏低好在发现得早。用别人的数据前先看清楚README和data processing描述看它标没标“already filtered”。2.2 三种常见数据形态以及各自的接法我去GEO上找单细胞数据最长遇到的是下面这三种形态10x标准输出的文件夹包含barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三个文件这是Cell Ranger的filtered_feature_bc_matrix输出。读到Seurat里很简单Read10X(data.dir 路径)一行搞定。h5格式文件10x的h5文件里同时包含了所有的barcode和feature信息。用Read10X_h5()读取注意函数会根据文件内容自动判断是gene还是gene expression读进来以后建议检查一下rownames有没有重复。纯文本格式的矩阵有的数据集会提供一个表达矩阵.csv或者表达矩阵.txt行是基因列是细胞barcode。这种情况下我习惯先read.table()读进来然后用CreateSeuratObject(counts mat, project 项目名)创建对象。这里有一个经常被忽略的细节GEO上很多数据集的表达矩阵提供给用户的可能是CPM、TPM或者log-normalized的形式而不一定是原始UMI count。如果你想自己走一遍完整的质控流程最好用原始counts数据重新构建。因为像PercentageFeatureSet算线粒体基因比例、CalculateBarcodeInflections判断空微滴阈值这些操作都是基于UMI count的统计性质来设计的你用标准化以后的数据跑结果会失真。注意如果数据集页面上已经给了“pipeline信息”和“filtered”字样说明数据提供方已经做过一轮质控。这种数据你仍然可以做二次质控但要把阈值适当放宽否则容易把本来是好的细胞也过滤掉。2.3 一个常用工具链的快速落地方案从下载到建好Seurat对象我常用的完整流程是这样的R环境下先下载数据文件用GEOquery包的getGEO()可以拿到GSE的软格式和平台信息但真正的大矩阵文件一般还是用download.file()手动下载比较可控。拿到压缩包以后解压如果是10x的三件套直接在R里跑library(Seurat) library(dplyr) # 读取10x标准输出 counts - Read10X(data.dir ./GSE123456_RAW/filtered_feature_bc_matrix/) # 创建Seurat对象min.cells 3表示基因至少在3个细胞里有表达 sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)min.cells 3这一步不是质控它只是在减少矩阵的稀疏程度——如果一个基因只在1个细胞里有表达它几乎没有统计效力留着还会拖慢后面的计算。min.features 200相当于最粗的一道门槛把连200个基因都没有的barcode先拿掉后面还会用更精细的指标来处理。如果是h5文件counts - Read10X_h5(./GSE123456_RAW/filtered_feature_bc_matrix.h5) sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)建好对象以后不要急着画图先加两个QC指标到meta.data里sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-)这里的^MT-是人类基因名的写法小鼠的话要用^mt-。为什么要单独算线粒体基因比例我下一节细讲。这里先记住一个原则线粒体比例是判断细胞状态的一个核心窗口不是在“找茬”而是在“察言观色”。如果是R环境中同时读了多个样本准备合并之后再统一质控建议先各自建好Seurat对象并算好percent.mt最后用merge()合并再一起做过滤。合并以后batch effect的问题另说但QC指标在合并前算好更准确。3. 三个核心指标nGene、nUMI、percent.mt到底在告诉你什么3.1 nUMI测到了多少分子nUMI指的是一个细胞内所有基因的UMIUnique Molecular Identifier总数也就是技术上去重以后这个细胞里总共被测到的转录本分子数量。UMI是10x平台实现“绝对定量”的关键——每个mRNA分子在反转录的时候都会被贴上一条随机的UMI序列测序以后比对上同一个基因、UMI相同的reads就合并成一条。所以nUMI的高低直接反映了这个细胞总体的转录本丰度。nUMI特别低通常意味着两种可能一是这个微滴里根本没包到细胞只有少量环境RNA二是包到了细胞但细胞已经破裂mRNA大部分降解了。不管是哪种这个barcode都不能算作一个“有效细胞”。反过来nUMI特别高也要注意后面我会讲到doublet的问题。为什么我不建议只靠nUMI一个指标做过滤因为不同类型的细胞转录本丰度天然差异很大。比如静止的T细胞和活跃的浆细胞前者nUMI可能只有后者的几分之一。你要是设一个很高的nUMI阈值等于把一群真实的小细胞全部干掉了。3.2 nGene检测到了多少基因nGene是检测到的基因数量也就是这个细胞里有多少种不同的基因被表达。它和nUMI高度相关但不完全是一回事。打个比方一个细胞测到了1000个UMI全部都来自Actb这一个基因那它的nGene是1另一个细胞同样测到1000个UMI但平均分布在1000个基因上那它的nGene是1000。显然前一个细胞大概率有问题——正常细胞不会只表达一个基因还这么极端。nGene偏低往往代表细胞质量不好或测序深度不够。但同样要注意不同细胞类型nGene的基线水平不一样。对一些低RNA含量的细胞类型比如静息状态的中性粒细胞nGene本来就不高硬设一个很高的阈值会误杀。我习惯把nGene的目标定在“大于500到1000”这个范围然后根据实际分布的拐点做调整。nGene还有一个用途是作为“测序深度”的估算参考。同一批细胞如果大部分nGene集中在4000到6000说明测序深度足够如果你看到nGene整体都在1000上下先别急着过滤更可能是建库或测序环节出了问题需要回到上游去看看。3.3 percent.mt线粒体基因比例为什么这么关键线粒体基因比例可能是单细胞质控里最有信息量、也最容易被误用的一个指标。原理是这样的正常活细胞的mRNA大部分来自核基因组编码的基因线粒体基因组的转录本只占一小部分。当细胞开始凋亡或者发生坏死时细胞膜通透性增加胞质内的大分子mRNA会降解而线粒体因为结构相对稳定其转录本反而保留得更多。所以凋亡/坏死细胞里的线粒体基因比例会异常升高。更微妙的一个情况是某些代谢旺盛的细胞比如心肌细胞、肝细胞天生线粒体转录本比例就偏高。所以线粒体比例阈值不是“一个标准走天下”一般在5%到20%之间浮动。我在处理肿瘤样本的时候常常会把阈值定到10%或15%因为肿瘤微环境里缺氧区域多细胞应激状态普遍死守“20%”这个数字可能会丢掉很多真实信号。提示算percent.mt时一定要确认基因命名规范。人类的MT-、小鼠的mt-大小写不同正则表达式写错就直接算不出来了。如果你处理的是非模式物种建议先去查线粒体基因组的注释文件把对应的基因名列出来再算。4. 阈值不能照抄要看数据自己说话4.1 为什么统一阈值常常翻车网上一搜单细胞质控最多的就是“keep cells with nFeature_RNA 500 percent.mt 20%”这种代码模板。不能说错但不同数据集之间差异极大直接用统一阈值很危险。我处理过一批小鼠脑组织的公共数据nGene中位数大约在1800percent.mt中位数低到2%这时候用500/20的阈值基本等于没过滤。但有一批某肿瘤单细胞数据nGene中位数只有600percent.mt中位数高达15%——因为样本是冻存组织细胞状态本来就差一些。这种情况下还用500/20的阈值哪怕过滤完后面的聚类里也会混入大量的低质量细胞。判断阈值应该看什么看分布。每个样本的nFeature_RNA、nCount_RNA、percent.mt分布都有自己的形状你要找的是分布里的“异常尾巴”。4.2 用分布图“读”数据而不是用眼睛猜我一般先把三个指标的分布画出来再决定阈值library(patchwork) p1 - VlnPlot(sce, features nFeature_RNA, pt.size 0) geom_hline(yintercept 500, linetype dashed, color red) p2 - VlnPlot(sce, features nCount_RNA, pt.size 0) geom_hline(yintercept 1000, linetype dashed, color red) p3 - VlnPlot(sce, features percent.mt, pt.size 0) geom_hline(yintercept 20, linetype dashed, color red) p1 p2 p3pt.size 0是关键反正单细胞数据几千上万个点全画出来就是一团黑干扰判断。小提琴图加一个阈值参考线一眼就能看出你的阈值线和分布主体隔了多远。除了小提琴图我还惯用FeatureScatter画出nFeature_RNA和nCount_RNA的关系FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA)正常情况下这两个指标呈明显的正相关——测序更深、测到的UMI更多检测到的基因数也会更多。如果有一批细胞nCount很高但nFeature很低那几乎可以肯定是doublet或者异常高表达某个基因的细胞。这个图上的离群点就是你需要review的重点对象。再有一个隐藏用法是看UMI和基因数量关系的“斜率”。不同测序深度下斜率会有变化但如果你发现某个样本的斜率明显偏离其他样本说明这个样本的测序质量或者细胞活性跟别的不一样后续分析需要区别对待。4.3 第一次过滤保守永远是更好的策略定阈值有一个原则第一轮过滤宁可保守不要激进。因为单细胞数据分析里有太多“后面还能补救”的环节但你一旦把一群真实细胞过滤掉后面就没有回来的路了。我一般用“双标准手动微调”的方式先用分布图确定一个大致的过滤区间比如nFeature_RNA 200 nFeature_RNA 6000 percent.mt 25跑完第一轮聚类之后再结合社群特异的QC指标做一轮“二次过滤”。这时候如果某个cluster普遍表现为nGene很低、percent.mt很高大概率是碎片或垂死细胞形成的cluster可以考虑单独去除。第一次过滤的代码长这样sce - subset(sce, subset nFeature_RNA 500 nFeature_RNA 6000 percent.mt 20)这里的上限6000不是绝对标准。我在一些高质量样本里见过nFeature_RNA到8000以上的真实巨噬细胞为什么因为巨噬细胞吞噬能力强胞内包含了大量外源RNA。这种时候上限要放宽到8000甚至10000免得把这类功能特殊的细胞全干掉。判断标准很简单——看你的分布图那个长尾巴到底是连续延伸的还是和主体明显断开。主体之外孤立的一小撮才是你真正需要警惕的“异常人群”。过滤完了不要立刻欢呼先看一眼过滤前后细胞数目的变化心里有个数。如果从一个样本3万个barcode过滤完只剩5000个你要先确认是不是数据本身有大量空微滴而不是你的阈值过严。过滤比例太大或太小都说明前面对数据的理解可能出了问题。5. 质控完别急着跑差异分析先做两件事5.1 数据质量的自检可视化过滤完成以后我建议先做一轮完整的可视化自检确认数据“在视觉上”是健康的。最基础的两个图一个是过滤前后三个指标的小提琴图对比另一个是过滤后的UMAP/tSNE初版聚类图。UMAP在过滤后立刻跑可能为时过早因为还没做归一化和高变基因筛选但你可以用Seurat快速走一遍标准pipeline看初版的cluster分布是否符合基本预期。做法很简单sce - NormalizeData(sce) sce - FindVariableFeatures(sce, selection.method vst, nfeatures 2000) sce - ScaleData(sce) sce - RunPCA(sce, npcs 30) sce - RunUMAP(sce, dims 1:20) sce - FindNeighbors(sce, dims 1:20) sce - FindClusters(sce, resolution 0.5) DimPlot(sce, label TRUE)这时候注意看两个信号第一有没有某个cluster几乎单独聚在图形边缘而它的QC指标明显和其他群体不同如果有回到VlnPlot去看这个cluster的nFeature和percent.mt大概率能找到问题。第二整体cluster分布是否过于碎片化如果resolution 0.5就分出30多个cluster而你的样本类型预期只有几种主要细胞类型先不要急着调resolution先把数据质量再看一遍。还有一个容易被忽略的自检项文库复杂度。用scuttle包或者手写代码估算一下每个细胞的cDNA文库复杂度如果整体偏低说明反转录效率不高这种数据后面做拟时序分析会很痛苦因为表达量的动态范围被压缩了。5.2 doublet检测一个经常被跳过的坑doublet问题在10x数据里非常普遍。10x官方给的doublet rate大约是每1000个细胞回收0.8%但实际建库过程中影响因素很多真实值往往更高。尤其是在上样量偏大的时候doublet率可以轻松超过5%。一个包含doublet的数据最典型的特征是某个细胞亚群表达两种完全不同的细胞类型marker。如果这种细胞数量够多它们会在UMAP上形成一个介于两个群体之间的“桥”导致你的细胞类型注释出现偏差——不是归到这一类就是归到那一类怎么解释都很别扭。质控阶段有没有办法处理doublet有但不是靠nFeature阈值。目前最常用的两个工具是DoubletFinderR包和ScrubletPython包原理类似用人工合成的doublet把两个细胞的表达谱按不同比例混合作为基准对每个barcode计算一个doublet score再用阈值判定哪些barcode最像doublet。以DoubletFinder为例在归一化和PCA之后调用library(DoubletFinder) # pK的选择是DoubletFinder最关键的一步 # 先跑一遍参数搜索选出pK值 sweep.res - paramSweep(sce, PCs 1:20, sct FALSE) sweep.stats - summarizeSweep(sweep.res, GT FALSE) bcmvn - find.pK(sweep.stats) # 根据bcmvn结果选最优pK然后预测doublet nExp_poi - round(0.05 * ncol(sce)) # 假设5% doublet率 sce - doubletFinder(sce, PCs 1:20, pN 0.25, pK 0.09, nExp nExp_poi, sct FALSE) # doublet列名一般是DF.classifications_0.25_0.09_xxx sce$doublet - scemeta.data[, grep(DF.classifications, colnames(scemeta.data))] sce - subset(sce, subset doublet Singlet)这里有一个经验教训DoubletFinder对pK值很敏感不跑参数搜索就直接用默认值效果差很多。而pK的搜索结果完全取决于你的数据分布不同样本要从头跑一遍。还有一个更朴素但很有效的doublet线索一个barcode的nUMI远超样本中位数比如5倍以上同时nFeature也超高但单独看每个指标又都在过滤范围内——这种“看着没问题但异常偏高”的barcode优先怀疑doublet。处理完doublet之后再看UMAP很多原本模糊的群体边界会变得清晰。关于doublet的阈值我的建议是“宁多杀勿放过”。因为doublet混在数据里对下游的影响是结构性的它会污染差异分析、干扰细胞类型鉴定。即便误杀了一小部分真实细胞损失也比留下doublet造成的伪迹小。这一步做完了你的数据才算真正达到了可以做下游分析的门槛。 单细胞数据质控实战：从GEO下载到阈值设定的完整指南
单细胞分析入门十个有八个把注意力都放在后面的聚类、找marker基因、做拟时序这些“高光步骤”上结果一跑起来就发现数据一团糟后面所有的图都解释不通。我最早做单细胞数据质控的时候就吃过这个亏拿到一个公共数据集看都没细看就直接往下走最后的UMAP图上细胞群体边界模糊得没法看重新回头补课才发现问题全出在源头。数据质控这件事听着基础实际上直接决定你后面所有分析靠不靠谱。这篇内容我不打算给你讲那些大而全的理论而是围绕一个完整的实操链路——从GEO下载数据、整理表达矩阵、跑质控指标、定阈值过滤到质控后的验证——把每一步为什么要这么做、怎么判断结果合不合理讲明白。适合刚接触单细胞分析、或者跑过流程但心里没底的朋友照着做能少走很多弯路。1. 为什么我把质控当作单细胞分析的“第一面墙”1.1 你的数据里混着什么质控要解决的三个真实问题单细胞测序和普通转录组测序有一个本质区别普通转录组测的是“一堆细胞平均下来”的表达量而单细胞测序拿到的是“一个一个细胞”的表达量。这个“一个一个”的实现方式决定了数据里天然会存在噪声。以10x Genomics平台的droplet-based方案为例细胞悬液和带有barcode的凝胶珠被一起包裹进微滴里每个微滴理论上只包一个细胞和一个凝胶珠。但实际运行的时候会有一批微滴包进去两个细胞doublet也有一批微滴根本没包到细胞只包到了游离的RNA也就是ambient RNA环境RNA。还有一些细胞在被捕获的时候就已经破裂死亡细胞膜不再完整细胞质里的RNA跑出去了一大部分测序仪最后捕获到的只是线粒体基因的残留片段。这三个问题对应到数据上就是三个不同维度的异常空微滴和低质量细胞会让一部分barcode对应的UMI数极低它们不配被当作一个“细胞”进入后续分析。doublet会造成一个barcode下混杂了两个不同细胞类型的表达谱聚类的时候它们可能单独形成一个奇怪的群体也可能把两个群体的边界模糊掉。破裂细胞和死亡细胞会因为mRNA降解呈现出线粒体基因比例异常偏高的特征。我见过很多人在质控这一步图省事用Seurat的默认参数跑一遍就过去了结果后面找marker基因的时候发现一堆未知基因注释细胞类型的时候怎么都对不上。实际上这些异常在质控阶段都有非常清晰的信号只是你没去看。1.2 质控不是“清垃圾”而是在定义“什么是细胞”我一直觉得质控这个翻译容易让人产生误解好像是把“坏”的东西扔掉。“坏”这个概念本身就很主观——一个线粒体基因比例20%的细胞在肿瘤组织里可能是一个真实存在的、处于应激状态的细胞你把它当成“坏死细胞”丢掉等于把真实的生物学信息也扔掉了。所以更准确的理解是质控是在定义“你的数据里什么样的barcode配叫一个可以进入分析的细胞”。你设定nGene 500意思就是“低于500个基因检测到的barcode我没有信心认为它是一个完整的细胞所以不纳入”。你设定percent.mt 20%意思就是“线粒体基因占比超过20%的barcode我倾向于认为是垂死细胞或者空微滴污染所以排除”。这个理解方式有两个实际好处。第一你不容易被网上的“标准阈值”带偏因为阈值本身是基于你的数据分布决定的而不是某个教程拍脑袋定的数字。第二当你跟审稿人或者合作者讨论的时候你能清楚地说出“我的细胞定义标准是什么”而不是含糊地说“按标准流程过滤”。我跟很多审稿人打过交道在方法学部分把质控阈值怎么定的写清楚是减少被质疑“数据质量存疑”的最有效方式。2. 从GEO下载到表达矩阵走到质控之前的路2.1 取数前先学会“读货”看懂GEO条目里的三个关键字段很多人一上来就想跑代码但GEO数据集很多时候并不会给你一个现成的“表达矩阵.csv”。你得花点时间把它的数据结构搞清楚不然下载完发现根本不是你要的东西白白浪费一个下午。看一个GEO条目比如GSE编号对应的页面我一般重点看三个字段。第一个是Overall design它告诉你这个数据集整体是怎么设计的有几组样本每组几个重复用的是哪个平台。第二个是Supplementary file这里才是真正的数据仓库——10x Genomics的数据通常会上传filtered_feature_bc_matrix、raw_feature_bc_matrix或者Cell Ranger处理好的h5文件如果是Smart-seq2这类全长转录组方案你拿到的可能是每个样本一个表达矩阵文件。第三个是SRA/ENA的原始数据入口如果你需要自己重新跑比对定量就得从这里拿FASTQ。尤其要注意的是GEO上有些数据集只有RAW data原始计数有些只提供了processed data已经过滤过的矩阵。我之前就踩过坑下载了一个已经做了过滤的矩阵自己又做了一轮过滤结果阳性细胞比例明显偏低好在发现得早。用别人的数据前先看清楚README和data processing描述看它标没标“already filtered”。2.2 三种常见数据形态以及各自的接法我去GEO上找单细胞数据最长遇到的是下面这三种形态10x标准输出的文件夹包含barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三个文件这是Cell Ranger的filtered_feature_bc_matrix输出。读到Seurat里很简单Read10X(data.dir 路径)一行搞定。h5格式文件10x的h5文件里同时包含了所有的barcode和feature信息。用Read10X_h5()读取注意函数会根据文件内容自动判断是gene还是gene expression读进来以后建议检查一下rownames有没有重复。纯文本格式的矩阵有的数据集会提供一个表达矩阵.csv或者表达矩阵.txt行是基因列是细胞barcode。这种情况下我习惯先read.table()读进来然后用CreateSeuratObject(counts mat, project 项目名)创建对象。这里有一个经常被忽略的细节GEO上很多数据集的表达矩阵提供给用户的可能是CPM、TPM或者log-normalized的形式而不一定是原始UMI count。如果你想自己走一遍完整的质控流程最好用原始counts数据重新构建。因为像PercentageFeatureSet算线粒体基因比例、CalculateBarcodeInflections判断空微滴阈值这些操作都是基于UMI count的统计性质来设计的你用标准化以后的数据跑结果会失真。注意如果数据集页面上已经给了“pipeline信息”和“filtered”字样说明数据提供方已经做过一轮质控。这种数据你仍然可以做二次质控但要把阈值适当放宽否则容易把本来是好的细胞也过滤掉。2.3 一个常用工具链的快速落地方案从下载到建好Seurat对象我常用的完整流程是这样的R环境下先下载数据文件用GEOquery包的getGEO()可以拿到GSE的软格式和平台信息但真正的大矩阵文件一般还是用download.file()手动下载比较可控。拿到压缩包以后解压如果是10x的三件套直接在R里跑library(Seurat) library(dplyr) # 读取10x标准输出 counts - Read10X(data.dir ./GSE123456_RAW/filtered_feature_bc_matrix/) # 创建Seurat对象min.cells 3表示基因至少在3个细胞里有表达 sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)min.cells 3这一步不是质控它只是在减少矩阵的稀疏程度——如果一个基因只在1个细胞里有表达它几乎没有统计效力留着还会拖慢后面的计算。min.features 200相当于最粗的一道门槛把连200个基因都没有的barcode先拿掉后面还会用更精细的指标来处理。如果是h5文件counts - Read10X_h5(./GSE123456_RAW/filtered_feature_bc_matrix.h5) sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)建好对象以后不要急着画图先加两个QC指标到meta.data里sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-)这里的^MT-是人类基因名的写法小鼠的话要用^mt-。为什么要单独算线粒体基因比例我下一节细讲。这里先记住一个原则线粒体比例是判断细胞状态的一个核心窗口不是在“找茬”而是在“察言观色”。如果是R环境中同时读了多个样本准备合并之后再统一质控建议先各自建好Seurat对象并算好percent.mt最后用merge()合并再一起做过滤。合并以后batch effect的问题另说但QC指标在合并前算好更准确。3. 三个核心指标nGene、nUMI、percent.mt到底在告诉你什么3.1 nUMI测到了多少分子nUMI指的是一个细胞内所有基因的UMIUnique Molecular Identifier总数也就是技术上去重以后这个细胞里总共被测到的转录本分子数量。UMI是10x平台实现“绝对定量”的关键——每个mRNA分子在反转录的时候都会被贴上一条随机的UMI序列测序以后比对上同一个基因、UMI相同的reads就合并成一条。所以nUMI的高低直接反映了这个细胞总体的转录本丰度。nUMI特别低通常意味着两种可能一是这个微滴里根本没包到细胞只有少量环境RNA二是包到了细胞但细胞已经破裂mRNA大部分降解了。不管是哪种这个barcode都不能算作一个“有效细胞”。反过来nUMI特别高也要注意后面我会讲到doublet的问题。为什么我不建议只靠nUMI一个指标做过滤因为不同类型的细胞转录本丰度天然差异很大。比如静止的T细胞和活跃的浆细胞前者nUMI可能只有后者的几分之一。你要是设一个很高的nUMI阈值等于把一群真实的小细胞全部干掉了。3.2 nGene检测到了多少基因nGene是检测到的基因数量也就是这个细胞里有多少种不同的基因被表达。它和nUMI高度相关但不完全是一回事。打个比方一个细胞测到了1000个UMI全部都来自Actb这一个基因那它的nGene是1另一个细胞同样测到1000个UMI但平均分布在1000个基因上那它的nGene是1000。显然前一个细胞大概率有问题——正常细胞不会只表达一个基因还这么极端。nGene偏低往往代表细胞质量不好或测序深度不够。但同样要注意不同细胞类型nGene的基线水平不一样。对一些低RNA含量的细胞类型比如静息状态的中性粒细胞nGene本来就不高硬设一个很高的阈值会误杀。我习惯把nGene的目标定在“大于500到1000”这个范围然后根据实际分布的拐点做调整。nGene还有一个用途是作为“测序深度”的估算参考。同一批细胞如果大部分nGene集中在4000到6000说明测序深度足够如果你看到nGene整体都在1000上下先别急着过滤更可能是建库或测序环节出了问题需要回到上游去看看。3.3 percent.mt线粒体基因比例为什么这么关键线粒体基因比例可能是单细胞质控里最有信息量、也最容易被误用的一个指标。原理是这样的正常活细胞的mRNA大部分来自核基因组编码的基因线粒体基因组的转录本只占一小部分。当细胞开始凋亡或者发生坏死时细胞膜通透性增加胞质内的大分子mRNA会降解而线粒体因为结构相对稳定其转录本反而保留得更多。所以凋亡/坏死细胞里的线粒体基因比例会异常升高。更微妙的一个情况是某些代谢旺盛的细胞比如心肌细胞、肝细胞天生线粒体转录本比例就偏高。所以线粒体比例阈值不是“一个标准走天下”一般在5%到20%之间浮动。我在处理肿瘤样本的时候常常会把阈值定到10%或15%因为肿瘤微环境里缺氧区域多细胞应激状态普遍死守“20%”这个数字可能会丢掉很多真实信号。提示算percent.mt时一定要确认基因命名规范。人类的MT-、小鼠的mt-大小写不同正则表达式写错就直接算不出来了。如果你处理的是非模式物种建议先去查线粒体基因组的注释文件把对应的基因名列出来再算。4. 阈值不能照抄要看数据自己说话4.1 为什么统一阈值常常翻车网上一搜单细胞质控最多的就是“keep cells with nFeature_RNA 500 percent.mt 20%”这种代码模板。不能说错但不同数据集之间差异极大直接用统一阈值很危险。我处理过一批小鼠脑组织的公共数据nGene中位数大约在1800percent.mt中位数低到2%这时候用500/20的阈值基本等于没过滤。但有一批某肿瘤单细胞数据nGene中位数只有600percent.mt中位数高达15%——因为样本是冻存组织细胞状态本来就差一些。这种情况下还用500/20的阈值哪怕过滤完后面的聚类里也会混入大量的低质量细胞。判断阈值应该看什么看分布。每个样本的nFeature_RNA、nCount_RNA、percent.mt分布都有自己的形状你要找的是分布里的“异常尾巴”。4.2 用分布图“读”数据而不是用眼睛猜我一般先把三个指标的分布画出来再决定阈值library(patchwork) p1 - VlnPlot(sce, features nFeature_RNA, pt.size 0) geom_hline(yintercept 500, linetype dashed, color red) p2 - VlnPlot(sce, features nCount_RNA, pt.size 0) geom_hline(yintercept 1000, linetype dashed, color red) p3 - VlnPlot(sce, features percent.mt, pt.size 0) geom_hline(yintercept 20, linetype dashed, color red) p1 p2 p3pt.size 0是关键反正单细胞数据几千上万个点全画出来就是一团黑干扰判断。小提琴图加一个阈值参考线一眼就能看出你的阈值线和分布主体隔了多远。除了小提琴图我还惯用FeatureScatter画出nFeature_RNA和nCount_RNA的关系FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA)正常情况下这两个指标呈明显的正相关——测序更深、测到的UMI更多检测到的基因数也会更多。如果有一批细胞nCount很高但nFeature很低那几乎可以肯定是doublet或者异常高表达某个基因的细胞。这个图上的离群点就是你需要review的重点对象。再有一个隐藏用法是看UMI和基因数量关系的“斜率”。不同测序深度下斜率会有变化但如果你发现某个样本的斜率明显偏离其他样本说明这个样本的测序质量或者细胞活性跟别的不一样后续分析需要区别对待。4.3 第一次过滤保守永远是更好的策略定阈值有一个原则第一轮过滤宁可保守不要激进。因为单细胞数据分析里有太多“后面还能补救”的环节但你一旦把一群真实细胞过滤掉后面就没有回来的路了。我一般用“双标准手动微调”的方式先用分布图确定一个大致的过滤区间比如nFeature_RNA 200 nFeature_RNA 6000 percent.mt 25跑完第一轮聚类之后再结合社群特异的QC指标做一轮“二次过滤”。这时候如果某个cluster普遍表现为nGene很低、percent.mt很高大概率是碎片或垂死细胞形成的cluster可以考虑单独去除。第一次过滤的代码长这样sce - subset(sce, subset nFeature_RNA 500 nFeature_RNA 6000 percent.mt 20)这里的上限6000不是绝对标准。我在一些高质量样本里见过nFeature_RNA到8000以上的真实巨噬细胞为什么因为巨噬细胞吞噬能力强胞内包含了大量外源RNA。这种时候上限要放宽到8000甚至10000免得把这类功能特殊的细胞全干掉。判断标准很简单——看你的分布图那个长尾巴到底是连续延伸的还是和主体明显断开。主体之外孤立的一小撮才是你真正需要警惕的“异常人群”。过滤完了不要立刻欢呼先看一眼过滤前后细胞数目的变化心里有个数。如果从一个样本3万个barcode过滤完只剩5000个你要先确认是不是数据本身有大量空微滴而不是你的阈值过严。过滤比例太大或太小都说明前面对数据的理解可能出了问题。5. 质控完别急着跑差异分析先做两件事5.1 数据质量的自检可视化过滤完成以后我建议先做一轮完整的可视化自检确认数据“在视觉上”是健康的。最基础的两个图一个是过滤前后三个指标的小提琴图对比另一个是过滤后的UMAP/tSNE初版聚类图。UMAP在过滤后立刻跑可能为时过早因为还没做归一化和高变基因筛选但你可以用Seurat快速走一遍标准pipeline看初版的cluster分布是否符合基本预期。做法很简单sce - NormalizeData(sce) sce - FindVariableFeatures(sce, selection.method vst, nfeatures 2000) sce - ScaleData(sce) sce - RunPCA(sce, npcs 30) sce - RunUMAP(sce, dims 1:20) sce - FindNeighbors(sce, dims 1:20) sce - FindClusters(sce, resolution 0.5) DimPlot(sce, label TRUE)这时候注意看两个信号第一有没有某个cluster几乎单独聚在图形边缘而它的QC指标明显和其他群体不同如果有回到VlnPlot去看这个cluster的nFeature和percent.mt大概率能找到问题。第二整体cluster分布是否过于碎片化如果resolution 0.5就分出30多个cluster而你的样本类型预期只有几种主要细胞类型先不要急着调resolution先把数据质量再看一遍。还有一个容易被忽略的自检项文库复杂度。用scuttle包或者手写代码估算一下每个细胞的cDNA文库复杂度如果整体偏低说明反转录效率不高这种数据后面做拟时序分析会很痛苦因为表达量的动态范围被压缩了。5.2 doublet检测一个经常被跳过的坑doublet问题在10x数据里非常普遍。10x官方给的doublet rate大约是每1000个细胞回收0.8%但实际建库过程中影响因素很多真实值往往更高。尤其是在上样量偏大的时候doublet率可以轻松超过5%。一个包含doublet的数据最典型的特征是某个细胞亚群表达两种完全不同的细胞类型marker。如果这种细胞数量够多它们会在UMAP上形成一个介于两个群体之间的“桥”导致你的细胞类型注释出现偏差——不是归到这一类就是归到那一类怎么解释都很别扭。质控阶段有没有办法处理doublet有但不是靠nFeature阈值。目前最常用的两个工具是DoubletFinderR包和ScrubletPython包原理类似用人工合成的doublet把两个细胞的表达谱按不同比例混合作为基准对每个barcode计算一个doublet score再用阈值判定哪些barcode最像doublet。以DoubletFinder为例在归一化和PCA之后调用library(DoubletFinder) # pK的选择是DoubletFinder最关键的一步 # 先跑一遍参数搜索选出pK值 sweep.res - paramSweep(sce, PCs 1:20, sct FALSE) sweep.stats - summarizeSweep(sweep.res, GT FALSE) bcmvn - find.pK(sweep.stats) # 根据bcmvn结果选最优pK然后预测doublet nExp_poi - round(0.05 * ncol(sce)) # 假设5% doublet率 sce - doubletFinder(sce, PCs 1:20, pN 0.25, pK 0.09, nExp nExp_poi, sct FALSE) # doublet列名一般是DF.classifications_0.25_0.09_xxx sce$doublet - scemeta.data[, grep(DF.classifications, colnames(scemeta.data))] sce - subset(sce, subset doublet Singlet)这里有一个经验教训DoubletFinder对pK值很敏感不跑参数搜索就直接用默认值效果差很多。而pK的搜索结果完全取决于你的数据分布不同样本要从头跑一遍。还有一个更朴素但很有效的doublet线索一个barcode的nUMI远超样本中位数比如5倍以上同时nFeature也超高但单独看每个指标又都在过滤范围内——这种“看着没问题但异常偏高”的barcode优先怀疑doublet。处理完doublet之后再看UMAP很多原本模糊的群体边界会变得清晰。关于doublet的阈值我的建议是“宁多杀勿放过”。因为doublet混在数据里对下游的影响是结构性的它会污染差异分析、干扰细胞类型鉴定。即便误杀了一小部分真实细胞损失也比留下doublet造成的伪迹小。这一步做完了你的数据才算真正达到了可以做下游分析的门槛。