拿到cellranger输出的filtered矩阵那一刻大多数新手的第一反应是接下来呢读进来、画两张图、然后呢我见过太多人卡在这一步——矩阵读进来了看着那个几万行乘几千列的表格完全不知道哪些细胞该留、哪些该扔更不知道后面那堆聚类分群代码跑出来到底在说什么。这篇就接着上一篇的进度把scRNA-seq数据分析里最关键的“数据清洗”和“细胞群初步分型”这部分讲透。说白了这一阶段做的所有事情都是为了让下游结论站得住脚如果输入的是脏数据后面UMAP上任何一个“分群”都可能是假象。这篇文章适合刚跑完上游比对、准备自己完成“从表达矩阵到细胞类型注释”这条完整链路的人。我会按实际操作的顺序把每一步为什么要这么做、参数怎么定、踩过哪些坑一次说清。1. 从计数矩阵到分析就绪清洗前先看懂数据长什么样1.1 上游产物到底给了你什么不管用的是10x Genomics的Cell Ranger还是其他平台最终拿到手的通常是一个基因×细胞的计数矩阵。以Cell Ranger的filtered_feature_bc_matrix为例它已经帮你过滤掉了一部分明显空液滴但“过滤掉空液滴”和“数据干净”之间还差着十万八千里。这个矩阵里的每个数字代表一个UMIUnique Molecular Identifier计数——就是这个细胞里某个基因被捕捉并测序到了几次。读取时用Seurat的话代码不复杂library(Seurat) library(dplyr) library(Matrix) data_dir - path/to/filtered_feature_bc_matrix counts - Read10X(data.dir data_dir) # barcode是列名细胞基因是行名 obj - CreateSeuratObject(counts counts, project my_sample, min.cells 3, min.features 200)这里有两个容易忽略的参数min.cells 3表示一个基因至少在3个细胞里检测到才保留目的是去掉那种只在单个细胞里零星出现的背景噪声也可能是测序错误min.features 200则是最起码有200个基因表达的细胞才保留这一步会直接扔掉一些明显异常的空泡。这两个值不是死规矩但作为起步标准非常稳。我个人的习惯是往后看分布再决定要不要收紧而不是一开始就卡死。1.2 三个质控指标UMI数、基因数、线粒体比例读进来之后下一步就是看每个细胞的质量。Seurat会自动帮你算nCount_RNA这个细胞所有基因UMI的总和和nFeature_RNA这个细胞检测到了多少个基因但线粒体基因比例得自己算obj[[percent.mt]] - PercentageFeatureSet(obj, pattern ^MT-)注意这个正则表达式是有物种差异的人用^MT-小鼠用^mt-。大小写搞错了比例算出来全是0我身边至少两个人在这上面浪费过一下午。这三个指标背后都是有生物学意义的不是随便挑出来的统计量nCount_RNA反映测序深度和RNA捕获效率。如果某个细胞的UMI总数低得离谱它很可能是空液滴或细胞已经破裂如果高得反常就要怀疑是不是两个细胞被当成一个测了。nFeature_RNA反映细胞转录组的复杂程度。测序深度低会导致基因检出数低但反过来一个细胞如果检出的基因数异常高也可能是双细胞。percent.mt是最经典的质量标签。线粒体基因在活细胞中通常占一小部分一般是5%~15%当细胞膜受损、胞浆RNA泄漏后细胞质里的mRNA丢失剩下的主要是线粒体RNA所以这个比例会飙升。看到percent.mt超过20%的细胞基本可以直接视为垂死状态。1.3 阈值怎么定画图比翻教程有用新手最喜欢问“阈值到底填多少”答案永远是“看你的数据分布”。直接用VlnPlot和FeatureScatter看全景VlnPlot(obj, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3, pt.size 0.1) FeatureScatter(obj, feature1 nCount_RNA, feature2 nFeature_RNA)正常的样本里大部分细胞的nFeature_RNA会集中在一个主峰附近旁边拖一条低质量细胞的尾巴。我在多数人源组织样本里的经验是nFeature_RNA下限设在200~300上限设在2500~4000percent.mt上限设在10%~20%。但如果你做的是肝脏、心脏这类代谢活跃的组织线粒体比例天然偏高一刀切15%可能会把一群本来健康的细胞全扔掉。这时候宁可把阈值放宽一点等聚类之后再看哪个cluster的线粒体比例异常高、单独处理。我自己的经验法则是先画图找分布的“膝盖点”elbow point而不是拍脑袋填一个数。过滤时也只做一次别来回反复试阈值——你试一次阈值就会引入一次偏见后面统计检验会失真。2. 清洗不止过滤低质量细胞双细胞、环境RNA和批次效应2.1 双细胞藏在“高质量”细胞里的内鬼很多新手过滤完低质量细胞就觉得万事大吉了但最危险的一类污染恰恰躲在高UMI的细胞里——双细胞也就是一个液滴里包了两个细胞。从质控指标看它可能完全“健康”基因数高、UMI数高、线粒体比例也不超标。但实际上它是两个人或者两个不同类型细胞的转录组掺在一起在聚类时会在两个真实群体之间制造一个虚假的过渡群。双细胞的识别需要用专门的工具我常用的是DoubletFinder。它的核心思路是人为在真实数据里注入模拟的双细胞然后用机器学习判断真实细胞更像哪一个“模拟双细胞”。实际操作里有两个坑第一个坑是参数pK的选择。新手直接跑默认参数经常报错或者效果很差正确做法是先做一遍参数扫描library(DoubletFinder) # 先做标准归一化和降维后面会讲 obj_scaled - NormalizeData(obj) %% FindVariableFeatures() %% ScaleData() %% RunPCA() # param sweep sweep.res - paramSweep(obj_scaled, PCs 1:20, sct FALSE) sweep.stats - summarizeSweep(sweep.res, GT FALSE) bcmvn - find.pK(sweep.stats) # 选择BCmetric最高的pK值 optimal_pk - as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))第二个坑是预期双细胞比率nExp_poi。10x官方有个大致估算每1000个细胞大约有0.8%的双细胞率但这只是文库制备时的物理参数实际还会受上样量影响。对于3万个细胞左右的样本我通常会按1%~2%估算。设得低了会漏掉设得高了会把真实细胞误杀。如果你用的是更复杂的组织比如肿瘤双细胞率会更高可以适当上调到2%~4%。2.2 环境RNA背景污染比你想的更普遍环境RNAambient RNA是液滴里游离的RNA片段主要来自制备过程中细胞破裂释放的mRNA。它在数据里的表现是几乎所有细胞都“低表达”很多不该表达的基因。如果你发现某个分化标志基因在所有cluster里都有1~2的log-normalized表达量却没有任何一个cluster明显高表达它大概率就是环境RNA污染。处理环境RNA有专门的工具比如SoupX和DecontX。但我给新手的建议是先观察别急着矫正。环境RNA的影响在差异分析阶段通常可以容忍因为所有细胞都受到同等污染组间比较时背景会互相抵消。只有当你要分析某个稀有群体的特异表达、或者文库的背景污染肉眼可见地严重时再考虑用SoupX做减法。矫正不当反而会引入新偏差比如把低表达基因的信息抹掉。2.3 批次效应合并样本前必须想清楚的问题如果你的项目里有多于一个样本逃不掉的话题就是批次效应。所谓批次效应指的是技术差异不同天建库、不同测序批次、操作人员不同导致的系统性表达差异它和生物学差异混在一起时会让UMAP上相同类型的细胞因为“出自同一批次”而各自抱团。判断是否存在批次效应最直观的方法是跑完聚类后用样本身份给细胞着色。如果UMAP上的每个cluster都混合着所有样本那是理想状态如果出现“这个cluster几乎全是样本A那个cluster几乎全是样本B”就要警惕。处理批次效应有两个主流路线Seurat的CCA整合适合跨数据集整合和注释转移它能找到不同批次间共享的细胞状态。Seurat v5的写法是IntegrateLayers整合后生成一个新assay比如integrated下游PCA聚类都基于这个assay。Harmony速度极快适合大批量样本几十个以上的10x文库。我个人在大规模整合时更喜欢Harmony因为CCA在样本数很多时计算量会爆炸而Harmony的迭代矫正逻辑非常稳。不管用哪个都要记住整合是在“规范化高变基因”做完之后、降维之前做的。而且整合完别用原始RNA assay去跑PCA要用整合后的assay——这是一条新手最容易踩的坑跑完发现cluster完全跟着样本走回头一看代码原来是忘了切换assay。3. 标准化、高变基因与降维为聚类铺一条干净的路3.1 归一化到底在做什么为什么不直接比较原始计数归一化这个问题我见过太多人完全理解反了。原始UMI计数最大的问题是不同细胞测序深度不同有的细胞测到了3万UMI有的只有3000如果直接拿原始数字比较测序深度高的细胞里所有基因都“看起来更高表达”。所以标准做法是把每个细胞的总UMI数拉平到同一个水平obj - NormalizeData(obj, normalization.method LogNormalize, scale.factor 10000)这行代码做的事情是每个基因的表达量除以该细胞总UMI数再乘以10000然后做log1p变换也就是log(x1)。这个变换有两个目的一是把数据从偏态分布拉成接近正态分布方便后续假设检验二是把动态范围压缩几个数量级让低表达基因也能参与下游比较。这也是为什么我前面说环境RNA的背景表达“1~2”是log-normalized值——它对应的原始比例其实是0.01%的量级。另一个选择是SCTransformsct它在归一化的同时还能建模UMI计数深度对表达的影响选高变基因也比默认的vst方法更稳健。但它的计算量比LogNormalize大不少而且在Seurat里一旦用了sct后面的FindVariableFeatures就可以跳过因为sct已经选过了变基因并写在模型里了。我的建议是新手先老老实实用LogNormalize流程跑通整个pipeline搞清楚每个对象里data、scale.data、counts三个槽分别是什么。等熟练了再换成sct也不迟。特别提醒sct不适合用来做细胞周期矫正也不要在同一个数据集里混用两种流程对比结果。3.2 高变基因、ScaleData和PCA之间的因果关系归一化做完表达矩阵里依然有几万个基因。这里有个核心问题如果直接拿全基因去算细胞间距离结果会被那些在不同细胞间变化不大、只在个别细胞里爆发的基因主导比如线粒体基因真实的细胞类型信号反而被淹没了。所以需要先筛出高变基因HVGsobj - FindVariableFeatures(obj, selection.method vst, nfeatures 2000)为什么是2000这不是生物学定律纯粹是计算精度和效率的折中。2000个高变基因通常已经能覆盖绝大部分决定细胞身份的信息再用更多基因边际收益很低计算代价却直线上升。你可以扩展到3000或者5000但很少看到有人用全基因组做下游聚类。选完高变基因接下来是ScaleData。这一步做了两件事对每个基因做z-score标准化减去均值除以标准差并把数据存进scale.data槽。为什么要z-score因为PCA对方差极其敏感如果一个基因的表达量范围是0到5另一个是0到1000不归一化的话后者会主导主成分方向。这里还有一个隐藏问题z-score会让所有基因的方差基本相同低表达基因的噪声也随之被放大所以ScaleData通常只对高变基因做正好衔接上一步。后面就是标准的降维三部曲obj - RunPCA(obj, npcs 30, verbose FALSE) ElbowPlot(obj, ndims 30)RunPCA把高维2000维的表达谱压缩成几十个主成分。ElbowPlot画出来之后你会看到一条从高处快速下降然后趋于平坦的曲线拐点之后的主成分解释的方差占比已经很低。通常我会选拐点附近再加两三个主成分比如拐点在12~15之间就用15或者20个PC。不是PC越多越好——加上那些接近噪声的维度聚类反而会分裂出虚假群体。4. 聚类与UMAP可视化细胞群的初始分群4.1 KNN图聚类与resolution参数的直觉降维得到的几十个主成分并不是终点它们是用来算细胞相似度的。FindNeighbors基于PC坐标构建KNN图——每个细胞只跟它最相似的若干细胞连边。FindClusters则在这个图上运行社区发现算法Seurat默认是Louvain后来版本也支持Leiden把连接紧密的节点分成一个个社区也就是我们说的cluster。resolution这个参数是新手最容易纠结的它的作用简单说就是控制分群粒度数字越大分的cluster越多越细。我用过的经验值供你参考resolution 0.1会把整个数据集分成很少的几个大类适合第一眼看全局。resolution 0.5~0.8这是做初步分型最常见的区间一般能得到8~15个cluster和细胞类型基本对应。resolution 1.2~2.0会把已识别的细胞类型进一步切成亚群适合后续精细分析。注意FindClusters跑出来的cluster编号是随机的不代表任何生物学顺序。cluster 0不等于“最重要的细胞”它只是社区发现算法遍历时第一个碰到的社区。很多新手会因为“cluster 0”在UMAP上占据最大面积就想当然认为它是某种主要细胞类型这完全是误会。4.2 UMAP和tSNE怎么选聚类完成后的可视化主流是UMAP和tSNE。两者的区别用大白话说UMAP更在意“全局结构”速度快在大数据集上能保持不同细胞类型之间的相对距离关系tSNE则更擅长把局部邻域关系展示得极度清晰但压缩全局结构、组间距离几乎没有意义。实际使用上我通常用UMAP做总览和审阅分群用tSNE出最终论文图因为它看起来更“干净”细胞簇边界更明显。跑起来也很简单obj - RunUMAP(obj, dims 1:15) obj - RunTSNE(obj, dims 1:15) DimPlot(obj, reduction umap, group.by seurat_clusters, label TRUE)dims的参数要和前面选的PC数目保持一致否则你等于换了输入数据在画图。这是个细节但每次有人跟我抱怨UMAP图“怎么长得跟前一次完全不一样”十有八九是RunPCA之后换过npcs忘了同步RunUMAP里的dims。4.3 聚类之后先别急着注释先“质检”有太多人跑完聚类就急着找marker、给细胞命名。但在我自己的流程里聚类完成后的第一件事是检查分群质量。具体看三件事每个cluster的QC指标是否正常。如果某个cluster的percent.mt中位数明显高于其他cluster它很可能是一群垂死细胞应该标记出来而不是硬着头皮注释。每个cluster是否被某个样本主导。画出DimPlot(obj, group.by orig.ident)如果看到颜色按cluster整片划分说明批次效应没处理好回头处理批次而不是继续往下走。有没有“小得可疑”的cluster。UMAP上偶尔会出现那种只有几个细胞的孤立小团它们通常是低质量细胞或技术噪声可以人工剔除也可以暂时保留但别指望能注释出什么可靠结论。这个“聚类后质检”的习惯能帮你省下后面至少一半的返工时间。5. marker基因识别与细胞群注释给cluster贴上身份标签5.1 FindAllMarkers的正确打开方式三看原则注释细胞群的第一步是找出每个cluster的特征基因Seurat的标准接口是FindAllMarkersobj_markers - FindAllMarkers(obj, only.pos TRUE, min.pct 0.25, logfc.threshold 0.25, test.use wilcox)这里的关键参数逐一拆解only.pos TRUE只关注在本cluster中上调的基因因为负marker某基因在其他cluster高表达对注释的指示意义弱而且会让输出量大增。min.pct 0.25基因至少要在cluster内和cluster外的25%细胞中有检测。这是一个“表达广度”的过滤条件防止某个基因只靠几个极值细胞撑起一张显著面孔。实际项目里0.25偏宽松我在粗略初筛时用0.25筛选候选marker时会提高到0.4~0.5。logfc.threshold 0.25log2倍数变化的门槛。注意这个值是log2尺度0.25对应的其实只有约1.19倍的表达差异也就是个非常宽松的起步值。真正有价值的marker我一般会在后续人工筛选中要求avg_log2FC 1甚至更高。test.use wilcoxWilcoxon秩和检验。它是非参数检验不假设数据正态分布而单细胞数据普遍是零膨胀的用t检验很容易被少数细胞带偏。默认用Wilcoxon是合理的选择也是我在真实项目中的选择。输出结果里有一列avg_log2FC和一列p_val_adj。很多新手只看p值觉得p越小基因越“标志性”。但在单细胞数据里p值受细胞数影响极大——一个cluster有5000个细胞另一个有50个细胞两者比较时本来就容易“显著”。我给自己定的判断标准是“三看”一看p_val_adj是否小于0.05基本门槛二看pct.1与pct.2的差距是否够大表达覆盖度差异例如pct.10.9而pct.20.1才是真正的cluster特异性而pct.10.9、pct.20.85说明这个基因只是普遍表达高算不上特征三看avg_log2FC是否够大表达量差异至少大于0.5严格时用1。# 查看cluster 0的top markers top0 - obj_markers %% filter(cluster 0) %% arrange(desc(avg_log2FC)) %% head(10)得到候选marker后不要直接照着一篇文献就把细胞类型定死。打开一个已知的marker基因列表去FeaturePlot上视觉确认一遍高表达的细胞是不是正好落在这个cluster里表达的细胞范围是否和UMAP上的分群轮廓吻合这一步虽然“土”但永远是最可靠的验证。5.2 常见组织里值得优先查的marker面板不同组织的marker体系差异很大但有几个通行的“默认面板”覆盖了最常见的细胞类型。以人和小鼠的血液、肿瘤组织为例我做了一张速查表细胞类型经典marker基因人备注T细胞CD3D, CD2, CD3E通用T细胞标记CD8 T细胞CD8A, CD8B在CD3阳性基础上细看CD4 T细胞CD4, IL7R注意常常是低表达NK细胞NKG7, GNLY, KLRD1, KLRC1和T细胞易混淆看NKG7与CD3D的互斥表达B细胞MS4A1(CD20), CD79A, CD79BB细胞发育阶段不同marker有差异浆细胞MZB1, IGHG1, SDC1(CD138)经常被误认成B细胞单核细胞CD14, LYZ, FCN1经典标志分布通常较广巨噬细胞C1QA, C1QB, MARCO, CD68注意巨噬和单核有连续状态树突状细胞CLEC9A(CD141), CLEC10A(CD303), ITGAX(CD11c)亚群多注释时不要贪心上皮细胞EPCAM, KRT19, KRT8, KRT18实体瘤样本里的主要群体内皮细胞VWF, PECAM1(CD31), CLDN5血管相关成纤维细胞COL1A1, COL1A2, DCN在基质中非常常见中性粒细胞S100A8, S100A9, FCGR3B(CD16b), CSF3R中性粒细胞捕获难度高易被过滤掉这里要特别提醒一句没有哪个基因是100%专一的。比如CD14在巨噬细胞里也很高NKG7在部分T细胞亚群中也会表达。注释时永远要组合判断至少两三个marker共同支持再给结论。单一marker就断言细胞类型基本是注释事故的前兆。5.3 人工注释与自动注释的配合现在也有不少自动注释工具比如SingleR、Garnett以及基于参考图谱的映射工具。对新手来说我建议把这些工具当作“预注释助手”而不是“最终裁判”。以SingleR为例library(SingleR) library(celldex) ref - HumanPrimaryCellAtlasData() singler_pred - SingleR(test GetAssayData(obj, assay RNA, slot data), ref ref, labels ref$label.main) obj$SingleR_classification - singler_pred$labels跑完之后把自动注释结果对照DimPlot看一遍如果某个cluster里绝大多数细胞被标成同一种类型那这是个很好的起点如果同一个cluster被拆得七零八落、标出了五六种完全不同的类型就要怀疑它根本是个混合群体或双细胞残留。即使自动注释结果看起来很合理我也会做两件事做最终确认一是对每个cluster取出top markers人工和已知marker面板比对二是画一个DoHeatmap把每个cluster的top marker表达量热图展示出来看看是否“各群界限分明”。格式如下top_markers - obj_markers %% group_by(cluster) %% top_n(5, avg_log2FC) DoHeatmap(obj, features unique(top_markers$gene), group.by seurat_clusters)如果热图上有大片横向条纹——也就是某个marker不仅在目标cluster高表达其他cluster也一片红——那说明这个cluster的marker特异性不够注释前要重新审视。5.4 初步分型不等于最终结论当所有cluster都有了细胞类型标签之后通常会把细粒度的cluster合并成生物学上更粗的类别。例如cluster 1、5、9都是T细胞亚群可以合并成一个T_cell大群。合并之后再回看UMAP确认每个大群之间边界清晰。这一步其实是在校验“初步分型”的稳定性。我遇到过一种情况某个cluster在第一次聚类时被分成了两半合并后重新跑UMAP发现两半其实分散到不同区域了——这说明之前的切分依据并不是真正的生物学差异可能只是某个混杂因素比如细胞周期在作祟。另外初步分型阶段千万不要过度解读。测序深度不均、某些稀有细胞类型本来就只有几十个细胞它们的注释置信度天然就低。对这类“边缘群体”我通常在报告里标注为“疑似XX细胞”或“状态待定”而不是硬给它一个确定的名字。数据清洗和细胞分型本身就是一环套一环渐进的一次流程跑完能锁定主要细胞群已经是很大的成功了。最后分享一个小经验注释是一个“越做越觉得自己之前是错的”的活。我回头看自己在同一个数据集上三个月前做的注释总能发现几个当时误定的群体。这不是坏事反而说明你在进步。分析单细胞数据最忌讳的就是“得出一个漂亮的结论就停手”多跑几个分辨率、多试几套参考注释、多画几张FeaturePlot结论才立得住。建议在这一步把每个版本的cluster标记和marker列表保存成CSV给打开的每个RDS命名都带上日期。不然等你三个月后再回来看这个项目数据集还在但当时的判断逻辑早就忘了。
单细胞转录组数据分析：从数据清洗到细胞初步分型全流程实操指南
拿到cellranger输出的filtered矩阵那一刻大多数新手的第一反应是接下来呢读进来、画两张图、然后呢我见过太多人卡在这一步——矩阵读进来了看着那个几万行乘几千列的表格完全不知道哪些细胞该留、哪些该扔更不知道后面那堆聚类分群代码跑出来到底在说什么。这篇就接着上一篇的进度把scRNA-seq数据分析里最关键的“数据清洗”和“细胞群初步分型”这部分讲透。说白了这一阶段做的所有事情都是为了让下游结论站得住脚如果输入的是脏数据后面UMAP上任何一个“分群”都可能是假象。这篇文章适合刚跑完上游比对、准备自己完成“从表达矩阵到细胞类型注释”这条完整链路的人。我会按实际操作的顺序把每一步为什么要这么做、参数怎么定、踩过哪些坑一次说清。1. 从计数矩阵到分析就绪清洗前先看懂数据长什么样1.1 上游产物到底给了你什么不管用的是10x Genomics的Cell Ranger还是其他平台最终拿到手的通常是一个基因×细胞的计数矩阵。以Cell Ranger的filtered_feature_bc_matrix为例它已经帮你过滤掉了一部分明显空液滴但“过滤掉空液滴”和“数据干净”之间还差着十万八千里。这个矩阵里的每个数字代表一个UMIUnique Molecular Identifier计数——就是这个细胞里某个基因被捕捉并测序到了几次。读取时用Seurat的话代码不复杂library(Seurat) library(dplyr) library(Matrix) data_dir - path/to/filtered_feature_bc_matrix counts - Read10X(data.dir data_dir) # barcode是列名细胞基因是行名 obj - CreateSeuratObject(counts counts, project my_sample, min.cells 3, min.features 200)这里有两个容易忽略的参数min.cells 3表示一个基因至少在3个细胞里检测到才保留目的是去掉那种只在单个细胞里零星出现的背景噪声也可能是测序错误min.features 200则是最起码有200个基因表达的细胞才保留这一步会直接扔掉一些明显异常的空泡。这两个值不是死规矩但作为起步标准非常稳。我个人的习惯是往后看分布再决定要不要收紧而不是一开始就卡死。1.2 三个质控指标UMI数、基因数、线粒体比例读进来之后下一步就是看每个细胞的质量。Seurat会自动帮你算nCount_RNA这个细胞所有基因UMI的总和和nFeature_RNA这个细胞检测到了多少个基因但线粒体基因比例得自己算obj[[percent.mt]] - PercentageFeatureSet(obj, pattern ^MT-)注意这个正则表达式是有物种差异的人用^MT-小鼠用^mt-。大小写搞错了比例算出来全是0我身边至少两个人在这上面浪费过一下午。这三个指标背后都是有生物学意义的不是随便挑出来的统计量nCount_RNA反映测序深度和RNA捕获效率。如果某个细胞的UMI总数低得离谱它很可能是空液滴或细胞已经破裂如果高得反常就要怀疑是不是两个细胞被当成一个测了。nFeature_RNA反映细胞转录组的复杂程度。测序深度低会导致基因检出数低但反过来一个细胞如果检出的基因数异常高也可能是双细胞。percent.mt是最经典的质量标签。线粒体基因在活细胞中通常占一小部分一般是5%~15%当细胞膜受损、胞浆RNA泄漏后细胞质里的mRNA丢失剩下的主要是线粒体RNA所以这个比例会飙升。看到percent.mt超过20%的细胞基本可以直接视为垂死状态。1.3 阈值怎么定画图比翻教程有用新手最喜欢问“阈值到底填多少”答案永远是“看你的数据分布”。直接用VlnPlot和FeatureScatter看全景VlnPlot(obj, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3, pt.size 0.1) FeatureScatter(obj, feature1 nCount_RNA, feature2 nFeature_RNA)正常的样本里大部分细胞的nFeature_RNA会集中在一个主峰附近旁边拖一条低质量细胞的尾巴。我在多数人源组织样本里的经验是nFeature_RNA下限设在200~300上限设在2500~4000percent.mt上限设在10%~20%。但如果你做的是肝脏、心脏这类代谢活跃的组织线粒体比例天然偏高一刀切15%可能会把一群本来健康的细胞全扔掉。这时候宁可把阈值放宽一点等聚类之后再看哪个cluster的线粒体比例异常高、单独处理。我自己的经验法则是先画图找分布的“膝盖点”elbow point而不是拍脑袋填一个数。过滤时也只做一次别来回反复试阈值——你试一次阈值就会引入一次偏见后面统计检验会失真。2. 清洗不止过滤低质量细胞双细胞、环境RNA和批次效应2.1 双细胞藏在“高质量”细胞里的内鬼很多新手过滤完低质量细胞就觉得万事大吉了但最危险的一类污染恰恰躲在高UMI的细胞里——双细胞也就是一个液滴里包了两个细胞。从质控指标看它可能完全“健康”基因数高、UMI数高、线粒体比例也不超标。但实际上它是两个人或者两个不同类型细胞的转录组掺在一起在聚类时会在两个真实群体之间制造一个虚假的过渡群。双细胞的识别需要用专门的工具我常用的是DoubletFinder。它的核心思路是人为在真实数据里注入模拟的双细胞然后用机器学习判断真实细胞更像哪一个“模拟双细胞”。实际操作里有两个坑第一个坑是参数pK的选择。新手直接跑默认参数经常报错或者效果很差正确做法是先做一遍参数扫描library(DoubletFinder) # 先做标准归一化和降维后面会讲 obj_scaled - NormalizeData(obj) %% FindVariableFeatures() %% ScaleData() %% RunPCA() # param sweep sweep.res - paramSweep(obj_scaled, PCs 1:20, sct FALSE) sweep.stats - summarizeSweep(sweep.res, GT FALSE) bcmvn - find.pK(sweep.stats) # 选择BCmetric最高的pK值 optimal_pk - as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))第二个坑是预期双细胞比率nExp_poi。10x官方有个大致估算每1000个细胞大约有0.8%的双细胞率但这只是文库制备时的物理参数实际还会受上样量影响。对于3万个细胞左右的样本我通常会按1%~2%估算。设得低了会漏掉设得高了会把真实细胞误杀。如果你用的是更复杂的组织比如肿瘤双细胞率会更高可以适当上调到2%~4%。2.2 环境RNA背景污染比你想的更普遍环境RNAambient RNA是液滴里游离的RNA片段主要来自制备过程中细胞破裂释放的mRNA。它在数据里的表现是几乎所有细胞都“低表达”很多不该表达的基因。如果你发现某个分化标志基因在所有cluster里都有1~2的log-normalized表达量却没有任何一个cluster明显高表达它大概率就是环境RNA污染。处理环境RNA有专门的工具比如SoupX和DecontX。但我给新手的建议是先观察别急着矫正。环境RNA的影响在差异分析阶段通常可以容忍因为所有细胞都受到同等污染组间比较时背景会互相抵消。只有当你要分析某个稀有群体的特异表达、或者文库的背景污染肉眼可见地严重时再考虑用SoupX做减法。矫正不当反而会引入新偏差比如把低表达基因的信息抹掉。2.3 批次效应合并样本前必须想清楚的问题如果你的项目里有多于一个样本逃不掉的话题就是批次效应。所谓批次效应指的是技术差异不同天建库、不同测序批次、操作人员不同导致的系统性表达差异它和生物学差异混在一起时会让UMAP上相同类型的细胞因为“出自同一批次”而各自抱团。判断是否存在批次效应最直观的方法是跑完聚类后用样本身份给细胞着色。如果UMAP上的每个cluster都混合着所有样本那是理想状态如果出现“这个cluster几乎全是样本A那个cluster几乎全是样本B”就要警惕。处理批次效应有两个主流路线Seurat的CCA整合适合跨数据集整合和注释转移它能找到不同批次间共享的细胞状态。Seurat v5的写法是IntegrateLayers整合后生成一个新assay比如integrated下游PCA聚类都基于这个assay。Harmony速度极快适合大批量样本几十个以上的10x文库。我个人在大规模整合时更喜欢Harmony因为CCA在样本数很多时计算量会爆炸而Harmony的迭代矫正逻辑非常稳。不管用哪个都要记住整合是在“规范化高变基因”做完之后、降维之前做的。而且整合完别用原始RNA assay去跑PCA要用整合后的assay——这是一条新手最容易踩的坑跑完发现cluster完全跟着样本走回头一看代码原来是忘了切换assay。3. 标准化、高变基因与降维为聚类铺一条干净的路3.1 归一化到底在做什么为什么不直接比较原始计数归一化这个问题我见过太多人完全理解反了。原始UMI计数最大的问题是不同细胞测序深度不同有的细胞测到了3万UMI有的只有3000如果直接拿原始数字比较测序深度高的细胞里所有基因都“看起来更高表达”。所以标准做法是把每个细胞的总UMI数拉平到同一个水平obj - NormalizeData(obj, normalization.method LogNormalize, scale.factor 10000)这行代码做的事情是每个基因的表达量除以该细胞总UMI数再乘以10000然后做log1p变换也就是log(x1)。这个变换有两个目的一是把数据从偏态分布拉成接近正态分布方便后续假设检验二是把动态范围压缩几个数量级让低表达基因也能参与下游比较。这也是为什么我前面说环境RNA的背景表达“1~2”是log-normalized值——它对应的原始比例其实是0.01%的量级。另一个选择是SCTransformsct它在归一化的同时还能建模UMI计数深度对表达的影响选高变基因也比默认的vst方法更稳健。但它的计算量比LogNormalize大不少而且在Seurat里一旦用了sct后面的FindVariableFeatures就可以跳过因为sct已经选过了变基因并写在模型里了。我的建议是新手先老老实实用LogNormalize流程跑通整个pipeline搞清楚每个对象里data、scale.data、counts三个槽分别是什么。等熟练了再换成sct也不迟。特别提醒sct不适合用来做细胞周期矫正也不要在同一个数据集里混用两种流程对比结果。3.2 高变基因、ScaleData和PCA之间的因果关系归一化做完表达矩阵里依然有几万个基因。这里有个核心问题如果直接拿全基因去算细胞间距离结果会被那些在不同细胞间变化不大、只在个别细胞里爆发的基因主导比如线粒体基因真实的细胞类型信号反而被淹没了。所以需要先筛出高变基因HVGsobj - FindVariableFeatures(obj, selection.method vst, nfeatures 2000)为什么是2000这不是生物学定律纯粹是计算精度和效率的折中。2000个高变基因通常已经能覆盖绝大部分决定细胞身份的信息再用更多基因边际收益很低计算代价却直线上升。你可以扩展到3000或者5000但很少看到有人用全基因组做下游聚类。选完高变基因接下来是ScaleData。这一步做了两件事对每个基因做z-score标准化减去均值除以标准差并把数据存进scale.data槽。为什么要z-score因为PCA对方差极其敏感如果一个基因的表达量范围是0到5另一个是0到1000不归一化的话后者会主导主成分方向。这里还有一个隐藏问题z-score会让所有基因的方差基本相同低表达基因的噪声也随之被放大所以ScaleData通常只对高变基因做正好衔接上一步。后面就是标准的降维三部曲obj - RunPCA(obj, npcs 30, verbose FALSE) ElbowPlot(obj, ndims 30)RunPCA把高维2000维的表达谱压缩成几十个主成分。ElbowPlot画出来之后你会看到一条从高处快速下降然后趋于平坦的曲线拐点之后的主成分解释的方差占比已经很低。通常我会选拐点附近再加两三个主成分比如拐点在12~15之间就用15或者20个PC。不是PC越多越好——加上那些接近噪声的维度聚类反而会分裂出虚假群体。4. 聚类与UMAP可视化细胞群的初始分群4.1 KNN图聚类与resolution参数的直觉降维得到的几十个主成分并不是终点它们是用来算细胞相似度的。FindNeighbors基于PC坐标构建KNN图——每个细胞只跟它最相似的若干细胞连边。FindClusters则在这个图上运行社区发现算法Seurat默认是Louvain后来版本也支持Leiden把连接紧密的节点分成一个个社区也就是我们说的cluster。resolution这个参数是新手最容易纠结的它的作用简单说就是控制分群粒度数字越大分的cluster越多越细。我用过的经验值供你参考resolution 0.1会把整个数据集分成很少的几个大类适合第一眼看全局。resolution 0.5~0.8这是做初步分型最常见的区间一般能得到8~15个cluster和细胞类型基本对应。resolution 1.2~2.0会把已识别的细胞类型进一步切成亚群适合后续精细分析。注意FindClusters跑出来的cluster编号是随机的不代表任何生物学顺序。cluster 0不等于“最重要的细胞”它只是社区发现算法遍历时第一个碰到的社区。很多新手会因为“cluster 0”在UMAP上占据最大面积就想当然认为它是某种主要细胞类型这完全是误会。4.2 UMAP和tSNE怎么选聚类完成后的可视化主流是UMAP和tSNE。两者的区别用大白话说UMAP更在意“全局结构”速度快在大数据集上能保持不同细胞类型之间的相对距离关系tSNE则更擅长把局部邻域关系展示得极度清晰但压缩全局结构、组间距离几乎没有意义。实际使用上我通常用UMAP做总览和审阅分群用tSNE出最终论文图因为它看起来更“干净”细胞簇边界更明显。跑起来也很简单obj - RunUMAP(obj, dims 1:15) obj - RunTSNE(obj, dims 1:15) DimPlot(obj, reduction umap, group.by seurat_clusters, label TRUE)dims的参数要和前面选的PC数目保持一致否则你等于换了输入数据在画图。这是个细节但每次有人跟我抱怨UMAP图“怎么长得跟前一次完全不一样”十有八九是RunPCA之后换过npcs忘了同步RunUMAP里的dims。4.3 聚类之后先别急着注释先“质检”有太多人跑完聚类就急着找marker、给细胞命名。但在我自己的流程里聚类完成后的第一件事是检查分群质量。具体看三件事每个cluster的QC指标是否正常。如果某个cluster的percent.mt中位数明显高于其他cluster它很可能是一群垂死细胞应该标记出来而不是硬着头皮注释。每个cluster是否被某个样本主导。画出DimPlot(obj, group.by orig.ident)如果看到颜色按cluster整片划分说明批次效应没处理好回头处理批次而不是继续往下走。有没有“小得可疑”的cluster。UMAP上偶尔会出现那种只有几个细胞的孤立小团它们通常是低质量细胞或技术噪声可以人工剔除也可以暂时保留但别指望能注释出什么可靠结论。这个“聚类后质检”的习惯能帮你省下后面至少一半的返工时间。5. marker基因识别与细胞群注释给cluster贴上身份标签5.1 FindAllMarkers的正确打开方式三看原则注释细胞群的第一步是找出每个cluster的特征基因Seurat的标准接口是FindAllMarkersobj_markers - FindAllMarkers(obj, only.pos TRUE, min.pct 0.25, logfc.threshold 0.25, test.use wilcox)这里的关键参数逐一拆解only.pos TRUE只关注在本cluster中上调的基因因为负marker某基因在其他cluster高表达对注释的指示意义弱而且会让输出量大增。min.pct 0.25基因至少要在cluster内和cluster外的25%细胞中有检测。这是一个“表达广度”的过滤条件防止某个基因只靠几个极值细胞撑起一张显著面孔。实际项目里0.25偏宽松我在粗略初筛时用0.25筛选候选marker时会提高到0.4~0.5。logfc.threshold 0.25log2倍数变化的门槛。注意这个值是log2尺度0.25对应的其实只有约1.19倍的表达差异也就是个非常宽松的起步值。真正有价值的marker我一般会在后续人工筛选中要求avg_log2FC 1甚至更高。test.use wilcoxWilcoxon秩和检验。它是非参数检验不假设数据正态分布而单细胞数据普遍是零膨胀的用t检验很容易被少数细胞带偏。默认用Wilcoxon是合理的选择也是我在真实项目中的选择。输出结果里有一列avg_log2FC和一列p_val_adj。很多新手只看p值觉得p越小基因越“标志性”。但在单细胞数据里p值受细胞数影响极大——一个cluster有5000个细胞另一个有50个细胞两者比较时本来就容易“显著”。我给自己定的判断标准是“三看”一看p_val_adj是否小于0.05基本门槛二看pct.1与pct.2的差距是否够大表达覆盖度差异例如pct.10.9而pct.20.1才是真正的cluster特异性而pct.10.9、pct.20.85说明这个基因只是普遍表达高算不上特征三看avg_log2FC是否够大表达量差异至少大于0.5严格时用1。# 查看cluster 0的top markers top0 - obj_markers %% filter(cluster 0) %% arrange(desc(avg_log2FC)) %% head(10)得到候选marker后不要直接照着一篇文献就把细胞类型定死。打开一个已知的marker基因列表去FeaturePlot上视觉确认一遍高表达的细胞是不是正好落在这个cluster里表达的细胞范围是否和UMAP上的分群轮廓吻合这一步虽然“土”但永远是最可靠的验证。5.2 常见组织里值得优先查的marker面板不同组织的marker体系差异很大但有几个通行的“默认面板”覆盖了最常见的细胞类型。以人和小鼠的血液、肿瘤组织为例我做了一张速查表细胞类型经典marker基因人备注T细胞CD3D, CD2, CD3E通用T细胞标记CD8 T细胞CD8A, CD8B在CD3阳性基础上细看CD4 T细胞CD4, IL7R注意常常是低表达NK细胞NKG7, GNLY, KLRD1, KLRC1和T细胞易混淆看NKG7与CD3D的互斥表达B细胞MS4A1(CD20), CD79A, CD79BB细胞发育阶段不同marker有差异浆细胞MZB1, IGHG1, SDC1(CD138)经常被误认成B细胞单核细胞CD14, LYZ, FCN1经典标志分布通常较广巨噬细胞C1QA, C1QB, MARCO, CD68注意巨噬和单核有连续状态树突状细胞CLEC9A(CD141), CLEC10A(CD303), ITGAX(CD11c)亚群多注释时不要贪心上皮细胞EPCAM, KRT19, KRT8, KRT18实体瘤样本里的主要群体内皮细胞VWF, PECAM1(CD31), CLDN5血管相关成纤维细胞COL1A1, COL1A2, DCN在基质中非常常见中性粒细胞S100A8, S100A9, FCGR3B(CD16b), CSF3R中性粒细胞捕获难度高易被过滤掉这里要特别提醒一句没有哪个基因是100%专一的。比如CD14在巨噬细胞里也很高NKG7在部分T细胞亚群中也会表达。注释时永远要组合判断至少两三个marker共同支持再给结论。单一marker就断言细胞类型基本是注释事故的前兆。5.3 人工注释与自动注释的配合现在也有不少自动注释工具比如SingleR、Garnett以及基于参考图谱的映射工具。对新手来说我建议把这些工具当作“预注释助手”而不是“最终裁判”。以SingleR为例library(SingleR) library(celldex) ref - HumanPrimaryCellAtlasData() singler_pred - SingleR(test GetAssayData(obj, assay RNA, slot data), ref ref, labels ref$label.main) obj$SingleR_classification - singler_pred$labels跑完之后把自动注释结果对照DimPlot看一遍如果某个cluster里绝大多数细胞被标成同一种类型那这是个很好的起点如果同一个cluster被拆得七零八落、标出了五六种完全不同的类型就要怀疑它根本是个混合群体或双细胞残留。即使自动注释结果看起来很合理我也会做两件事做最终确认一是对每个cluster取出top markers人工和已知marker面板比对二是画一个DoHeatmap把每个cluster的top marker表达量热图展示出来看看是否“各群界限分明”。格式如下top_markers - obj_markers %% group_by(cluster) %% top_n(5, avg_log2FC) DoHeatmap(obj, features unique(top_markers$gene), group.by seurat_clusters)如果热图上有大片横向条纹——也就是某个marker不仅在目标cluster高表达其他cluster也一片红——那说明这个cluster的marker特异性不够注释前要重新审视。5.4 初步分型不等于最终结论当所有cluster都有了细胞类型标签之后通常会把细粒度的cluster合并成生物学上更粗的类别。例如cluster 1、5、9都是T细胞亚群可以合并成一个T_cell大群。合并之后再回看UMAP确认每个大群之间边界清晰。这一步其实是在校验“初步分型”的稳定性。我遇到过一种情况某个cluster在第一次聚类时被分成了两半合并后重新跑UMAP发现两半其实分散到不同区域了——这说明之前的切分依据并不是真正的生物学差异可能只是某个混杂因素比如细胞周期在作祟。另外初步分型阶段千万不要过度解读。测序深度不均、某些稀有细胞类型本来就只有几十个细胞它们的注释置信度天然就低。对这类“边缘群体”我通常在报告里标注为“疑似XX细胞”或“状态待定”而不是硬给它一个确定的名字。数据清洗和细胞分型本身就是一环套一环渐进的一次流程跑完能锁定主要细胞群已经是很大的成功了。最后分享一个小经验注释是一个“越做越觉得自己之前是错的”的活。我回头看自己在同一个数据集上三个月前做的注释总能发现几个当时误定的群体。这不是坏事反而说明你在进步。分析单细胞数据最忌讳的就是“得出一个漂亮的结论就停手”多跑几个分辨率、多试几套参考注释、多画几张FeaturePlot结论才立得住。建议在这一步把每个版本的cluster标记和marker列表保存成CSV给打开的每个RDS命名都带上日期。不然等你三个月后再回来看这个项目数据集还在但当时的判断逻辑早就忘了。
