尧图网页设计 
做单细胞数据分析也有几年了从最早拿到10X数据手足无措到现在能一口气从FASTQ原始数据跑到UMAP出图中间踩过的坑不少但也沉淀了一套相对稳定、可以复用的完整流程。这篇文章我就把整套思路和实操细节写出来从最原始的FASTQ文件开始一步步走到最终的可视化结果希望能让刚入坑的朋友少走点弯路。我先说一下这篇文章适合谁看。如果你刚拿到10X Genomics的测序数据对Cell Ranger和Seurat只是听说过但不知道怎么串起来或者你已经跑通了流程但总觉得结果质量不太对劲、UMAP图片不理想那这篇文章应该很适合你。我默认你有基本的Linux命令行操作能力R语言有一定基础但即便是小白按着步骤一步步来也能跑出比较靠谱的结果。1. 整体流程拆解从FASTQ到UMAP中间到底发生了什么1.1 单细胞数据的特殊之处在讲具体操作之前我得先花点篇幅讲讲单细胞数据为什么不能像普通转录组那样直接处理。普通转录组RNA-seq拿到的是一堆细胞混在一起的平均表达值而10X Genomics技术解决的核心问题是给每一个细胞打上独特的身份证再给每一条转录本打上分子条码这样测序数据回来后我们就能从海量reads里还原出每个细胞里各有多少条mRNA。这里面有两个概念你必须彻底理解因为后面所有质控都围绕它们展开Cell Barcode细胞条码一段特异的序列用来标记mRNA来自哪一个细胞。UMIUnique Molecular Identifier唯一分子标识符一段随机序列用来区分同一条mRNA的重复扩增拷贝。简单打个比方细胞条码是门牌号告诉你这条数据是哪一户的UMI是购物小票就算同一件商品同一条mRNA被复制了好多份你也知道它本质上只买过一次。Cell Ranger这个工具干的事情就是把测序得到的海量reads按照这两重信息拆解成一张细胞 x 基因的数字矩阵。1.2 为什么选Cell Ranger而不是自己写脚本比对可能有人会问我能不能拿常规的STAR、hisat2去比对然后自己用featureCounts计数理论上可以但实践上非常不建议。原因有三个第一10X数据的文库结构特殊read 1的前16个碱基是细胞条码接下来10个碱基是UMI后面才是转录本序列。如果不用专门工具你得自己花很大力气去处理这三段信息的拆分而且容易出错。第二Cell Ranger内置了一个重要的纠错逻辑它会对细胞条码和UMI做校正。测序过程中会出现碱基错读如果一个细胞条码只差一个碱基Cell Ranger会认为它们来自同一个细胞这能显著降低数据噪音。自己写的脚本基本没有这个能力。第三Cell Ranger比对用的STAR比对器做了专门优化比对速度极快。我之前试过用常规STAR去比对一个8G左右的FASTQ文件跑了快十个小时而Cell Ranger只用不到三个小时差距非常大。从我的经验来看你就老老实实用官方工具跑通整个pipeline把精力省下来去做下游分析这才是正路。1.3 整体流程的五个核心节点从原始FASTQ到最终的UMAP可视化整个流程可以拆成五个关键节点节点输入核心工具主要产出数据质控与格式检查FASTQ原始文件fastqc / multiqcQC报告定量分析FASTQ 参考转录组Cell Ranger表达矩阵、BAM文件、质控指标数据过滤与预处理表达矩阵Seurat过滤后的Seurat对象降维聚类预处理后的矩阵Seurat harmony可选PCA、聚类结果可视化与解释聚类结果Seurat ggplot2UMAP图、marker基因图这个流程的主干是不变的每个环节都有它自己容易踩的坑。下面我逐个环节展开讲重点说那些我在实际项目中反复遇到、容易出问题的地方。2. 起点决定终点FASTQ数据准备与质控细节2.1 拿到数据的第一个动作检查完整性很多新手拿到数据就急着跑Cell Ranger结果跑到一半发现FASTQ文件不完整或者文件名不对白白浪费时间。我习惯的做法是拿到数据后先做三件事。第一确认FASTQ文件数量是否正确。10X的数据通常每个样本有四个文件R1、R2分别对应两条测序readI1、I2是index read。如果你用了双index的建库方式I1和I2都必须存在。少一个文件整个样本就作废了。第二检查文件大小是否合理。一般的10X单细胞数据一个样本的FASTQ文件加起来少则有几个GB多则几十个GB。如果你拿到手的文件只有几百MB非常可疑建议用zcat查看一下reads数是否达标。第三检查文件名是否符合Cell Ranger的命名规则。如果你用官方的cellranger mkfastq做过拆分文件名一般是SampleName_S1_L001_R1_001.fastq.gz这样的格式。如果是从其他渠道拿来的数据文件名可能比较乱建议先规范化命名免得后面启动cellranger count的时候识别不到样本。2.2 测序质量到底要不要先过滤FASTQ拿到手之后一个常见疑问是要不要先用fastp、trimmomatic之类的工具做一遍质量过滤再跑Cell Ranger我的答案是不需要甚至可以明确说不要这么做。原因在于Cell Ranger自身的pipeline中已经内置了read质量过滤的步骤而且它的过滤逻辑是专门针对单细胞文库结构设计的。你在外面再做一遍过滤反而可能把带有细胞条码信息的read误杀导致有效细胞数下降。但我不反对你做质控报告。用fastqc和multiqc看看整体测序质量的分布从大局上有帮助如果Q30比例过低比如低于80%你心里要有数后面检查比对率和有效barcode比例的时候重点留意。测序质量不好的数据不是简单的过滤能救回来的这一点要有心理准备。注意Cell Ranger自带的速度很快且专门针对单细胞测序的文库结构做了优化不要在外部再做额外的预处理。2.3 参考转录组的选择和准备参考基因组的选择直接影响比对结果和基因定量准确性这个环节值得仔细对待。人源样本我用的是10X官网提供的refdata-cellranger-GRCh38-2020-A。这个版本的参考转录组包含GENCODE v32的基因注释是官方构建好的里面已经包含pre-mRNA序列对于检测内含子区域很有帮助。为什么单细胞要用包含pre-mRNA的参考因为单细胞测序捕获的是细胞裂解出来的RNA大量mRNA还没来得及剪接完毕如果参考转录组里没有内含子序列这些未剪接转录本会被白白浪费掉直接导致检测到的基因数和UMI数偏低。具体在Cell Ranger里怎么构建参考转录组官方文档写得很清楚我在这里就不重复了。只说一个关键点优先用官方构建好的参考转录组而不是自己从GENCODE/Gencode下载GTF和FASTA再跑cellranger mkref。因为官方构建过程做了很多隐含的优化比如处理了基因ID和转录本版本的对应关系这些都是自己构建时容易出错的地方。人源之外的其他物种官方没有现成的参考那就只能自己跑cellranger mkref了。构建参考时要注意GTF文件与FASTA文件的版本一致性比如GRCm39就配GRCm39的GTF混搭不同版本很容易出现比对率极低、注释率极低的问题而且这类问题非常隐蔽不为新手常见。3. Cell Ranger实战一次跑通的6个关键设置3.1 cellranger count命令的完整参数解析假设你的工作目录是/data/scRNAFASTQ文件在/data/scRNA/FASTQ目录下样本名为Tumor01我们要跑的命令大概长这样cd /data/scRNA cellranger count \ --idTumor01_hg38 \ --transcriptome/data/ref/refdata-cellranger-GRCh38-2020-A \ --fastqs/data/scRNA/FASTQ \ --sampleTumor01 \ --expect-cells5000 \ --localcores16 \ --localmem64这里我逐个参数讲一下我的设置习惯。--expect-cells这个参数比较关键。它告诉Cell Ranger这个样本大概有多少个细胞会影响Cell Ranger判断有效barcode的阈值。我一般会先看看建库时候的目标捕获数通常目标在8000左右我就填--expect-cells8000。如果填得太保守比如实际有10000个细胞但只填了3000那高表达的细胞没什么问题但低表达的细胞群体可能会被当成背景噪音丢掉。反之如果填得比实际高很多背景的barcode会被大量保留下来下游分析会非常痛苦。拿不准的时候宁可填得稍微高一点也不要低估因为下游过滤时低质量细胞还能通过QC步骤去掉但如果细胞压根没被识别到那就彻底丢了。--localcores和--localmem是控制计算资源分配的。这里有个容易出问题的点--localmem的单位是GB不是MB。我有段时间没注意填了--localmem32000结果Cell Ranger把它当成32TB内存来用直接跑到OOM被系统kill了。正确做法是看机器空闲内存比如机器有128GB内存就给Cell Ranger分配64到96GB。核心数同理不要贪心留出几个核给系统的其他进程免得互相抢资源。3.2 比对环节发生了什么为什么它这么重要cellranger count内部的比对环节是基于STAR比对器的。这一步干的事是把read 2上测到的转录本序列比对到参考基因组上而read 1上的细胞条码和UMI基本上直接被保留和使用不做太多处理。这里有个非常重要的细节Cell Ranger默认允许一个转录本比对到基因组的多处位置。这个设计是深思熟虑的——因为参考基因组里有很多来自同一个基因家族的重复序列如果一个read只比对到一个位置而不管其他位置很多基因的定量就会偏低。Cell Ranger在比对完成后会根据UMI数去判断这条read到底更可能来自哪个位置这个逻辑比简单粗暴的唯一比对要合理得多。比对率是考察数据质量的核心指标之一。正常情况下人源样本的比对率应该在85%到95%之间。如果你的比对率低于80%先检查参考转录组是否存在版本问题其次考虑样本是否有其他物种污染。3.3 看完web_summary.html才算跑完跑完cellranger count并不是万事大吉。每个样本跑完后都会在outs目录下生成一个web_summary.html文件这个文件包含了全套质控指标是我必看的内容。我一般会重点看几个数Estimated Number of Cells评估出的细胞数和预期差距过大就要警惕。Mean Reads per Cell每个细胞的平均reads数人源样本要做到20000以上比较理想但如果你测序深度不高只要后续分析能出结果也不是绝对的。Median Genes per Cell每个细胞检测到的中位基因数人源样本一般要在2000到5000之间如果低于1000说明数据质量堪忧。Q30 Bases in Barcode细胞条码区域的Q30比例这个数值低说明建库环节可能有问题。Valid Barcodes有效barcode比例低于80%说明可能有index错配或者建库污染。需要说明的是这些指标没有一个绝对的标准需要结合你实际的实验设计、物种、建库方案来综合判断。我之前有一个样本中位基因数只有1200左右一开始以为质量很差后来发现是做了神经元分选后的样本——神经元这种细胞本身就是低复杂度RNA表达这个数值是正常的。3.4 Cell Ranger的结果目录里到底有什么跑完Cell Ranger后outs目录下会生成多个文件。对下游分析我主要用到这几个文件路径内容用途outs/filtered_feature_bc_matrix/过滤后的细胞-基因表达矩阵Seurat分析的主要输入outs/raw_feature_bc_matrix/未过滤的全部barcode矩阵用于自定义过滤不常用outs/analysis/聚类、差异表达等初步分析快速预览结果outs/bam/比对后的BAM文件重新定量或可视化时用下游分析的核心输入基本就是filtered_feature_bc_matrix这个目录它包含三个文件matrix.mtx.gz、features.tsv.gz和barcodes.tsv.gz。这三个文件是一个标准的稀疏矩阵存储格式Seurat可以直接读取。你要理解这三样东西的匹配关系barcodes.tsv是细胞ID列表features.tsv是基因ID列表matrix.mtx.gz是它们的行和列对应的计数矩阵。4. 进入Seurat表达矩阵读取与质量控制在细节里4.1 用Read10X快速读取矩阵数据进入R环境后读取数据这一步很简单但有些细节值得注意。我习惯这样写library(Seurat) library(dplyr) library(ggplot2) data_dir - outs/filtered_feature_bc_matrix data - Read10X(data.dir data_dir) # 这里特别注意如果有多个样本需要合并先各自创建对象再merge expr - CreateSeuratObject( counts data, project Tumor01, min.cells 3, min.features 200 )min.cells 3的意思是一个基因至少在3个细胞中有表达才被保留这个参数能过滤掉那些只在极少数细胞里出现的低质量基因min.features 200的意思是一个细胞至少检测到200个基因才被保留这个阈值过滤掉了那些可能没有细胞核或RNA含量极低的空液滴。特别强调一下读取路径的问题Read10X函数的data.dir参数要求的是包含三个文件的目录路径不是matrix.mtx.gz的文件路径。很多人第一次用的时候把路径写到了具体的文件层级结果报错找不到数据这个坑不算大但很烦。如果你有多个样本要合并分析我建议先分别创建Seurat对象再用merge()函数合并。在合并后加group信息后面分析批次效应时好区分。4.2 线粒体基因比例、核糖体基因比例和双细胞识别Seurat对象创建好之后第一个核心步骤是计算质控指标。这里最常用的三个指标是nCount_RNA每个细胞的UMI总数、nFeature_RNA每个细胞检测到的基因数和percent.mt线粒体基因占所有UMI的比例。# 计算线粒体基因比例人源样本用MT-前缀小鼠用mt-前缀 expr[[percent.mt]] - PercentageFeatureSet(expr, pattern ^MT-) # 核糖体基因比例可以过滤掉某些极端情况 expr[[percent.ribo]] - PercentageFeatureSet(expr, pattern ^RP[SL]) # 查看QC指标的分布 VlnPlot(expr, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)这三个指标背后的生物学意义你要理解nFeature_RNA和nCount_RNA过低细胞可能已经破裂RNA流出去了或者这个液滴是空的。nCount_RNA很高但nFeature_RNA很低可能是某个低复杂度细胞类型比如红细胞但也可能是文库污染。percent.mt过高细胞处于应激或即将凋亡的状态。因为胞质mRNA在破裂时会被降解而线粒体mRNA由于包裹在线粒体内部反而相对稳定。所以一个濒临死亡的细胞会表现出高线粒体基因比例。但是做好功课某些组织类型比如脑组织里的少突胶质细胞天然线粒体比例偏高一刀切过滤容易误杀。这时候多看几组数据结合解剖位置判断。关于双细胞doublets简单说一下两个细胞被包在同一个液滴里被测序这种情况下nCount和nFeature通常都会异常高。你可以先用简单的阈值卡掉nFeature超过一定值的细胞比如人源样本6000但如果需要更严格的过滤建议用DoubletFinder这样的专门工具。后面我会在常见问题里再展开。4.3 过滤阈值怎么定看分布而不是拍脑袋过滤阈值是很多新手特别纠结的地方。我见过有人直接拿教程里的nFeature_RNA 500 nFeature_RNA 5000 percent.mt 20套用到自己的数据上结果效果很差。正确的做法是先画图看分布再根据分布形态确定阈值。# 画QC指标分布 plot1 - FeatureScatter(expr, feature1 nCount_RNA, feature2 percent.mt) plot2 - FeatureScatter(expr, feature1 nCount_RNA, feature2 nFeature_RNA) plot1 plot2比如某个样本的nFeature_RNA分布主峰在2000到4000之间然后有一个小尾巴延伸到8000以上那这个8000以上就很可能是双细胞阈值应该设在8000附近或略低一些。再比如percent.mt在主峰附近如果集中在10%以下那阈值设在15%或20%就比较合理。但如果你的样本有某种细胞类型天然线粒体高你看到分布的主峰直接在15%到20%之间那阈值设在25%甚至30%也合理。你得记住一个核心理念过滤器只是去除明显的异常值而不是做精细生物学划分。真正决定细胞类型界限的是后续的聚类分析不要在过滤环节下手太狠。5. 归一化与高变基因选择让数据说话的前提5.1 LogNormalize和SCTransform怎么选过滤完之后下一步是归一化。我用Seurat的时候面对两个选项传统的LogNormalize和推荐的SCTransform。LogNormalize的逻辑很直观每个细胞的UMI总数标准化到相同的深度默认10000然后做log变换。这个方法的优点是简单、快、结果也稳定缺点是它没有完全考虑转录本的技术噪音——高表达基因在细胞间的波动天然比低表达基因大LogNormalize没法很好地处理这个问题。SCTransform做了一件更聪明的事它对每个基因建立了一个广义线性模型把UMI数和建库深度之间的关系拟合出来然后把深度效应直接去掉。这样归一化后的数据能更好地反映真实的生物学差异而不是技术操作差异。我的经验是如果你的数据量不大比如小于2万个细胞直接用SCTransform替代整个NormalizeData和ScaleData步骤通常结果更干净。如果数据量很大超过5万个细胞SCTransform的计算时间会明显拉长这时候也可以先跑LogNormalize走一遍流程看看聚类结果再说。还有一个折中方案用LogNormalize做初步的聚类拿到结果后对感兴趣的细胞亚群用SCTransform再精修一遍这样速度和精度都能兼顾。# 方法一传统LogNormalize流程 expr - NormalizeData(expr, normalization.method LogNormalize, scale.factor 10000) expr - FindVariableFeatures(expr, selection.method vst, nfeatures 2000) expr - ScaleData(expr) # 方法二SCTransform流程 expr - SCTransform(expr, vars.to.regress percent.mt, verbose FALSE)SCTransform里vars.to.regress参数可以顺带回归掉技术变量比如线粒体基因比例。这里要注意一个原则不要随便回归掉你认为不重要的生物学变量。很多人看到percent.mt高第一反应是把它回归掉但线粒体压力所反映的细胞状态本身可能是有生物学意义的把它回归掉相当于抹掉了这部分信息。有争议的情况下先不做回归看看聚类结果再说。5.2 高变基因2000个还是更多的选择逻辑FindVariableFeatures这一步选出的高变基因HVGs决定了你后续PCA分析能看到什麼样的坐标系。默认是找2000个高变基因这对于大部分样本已经够用。但也有例外如果你研究的样本组织复杂度非常高比如脑、肝脏这种有多种细胞类型且每种都有特异表达基因的组织可以考虑提高到3000或4000个让更多稀有细胞类型的特异基因入选。这个逻辑不难理解高变基因就是PCA降维时用来构建主成分的特征集如果你的稀有细胞类型的关键特异基因不在这个特征集里那这个细胞群在聚类时就可能被其他细胞群的特征淹没难以分离。选基因之后ScaleData这一步会把基因表达值标准化为z-score均值为0方差为1并且默认对每个基因做一次线性回归去掉细胞测序深度的影响。这一步在SCTransform流程里已经被包含进去了所以如果你用了SCTransform就不要再重复调用ScaleData否则会出问题。# 传统流程需要单独ScaleData expr - ScaleData(expr, vars.to.regress nCount_RNA)需要留意的是如果数据中某些基因的表达量全为0比如过滤后有些基因存在于矩阵中但所有细胞都不表达ScaleData会提示错误或给出警告这个问题出现频率不低。解决办法是在ScaleData之前先去掉全0基因Seurat现在的版本通常会帮你处理但如果你用很老的版本还是要手动检查。6. 降维聚类的完整实操从PCA到UMAP的关键决策6.1 PCA分析主成分数量不该拍脑袋决定PCA降维是整个流程中最技术的一步因为它的结果直接影响后续聚类质量。核心问题是保留多少个主成分PC用于后续分析这个决定不能拍脑袋要靠证据。首先明确PCA在干什么我们有了高变基因矩阵但维度还是很高几千个基因直接做聚类计算量巨大且容易受噪音干扰。PCA把几千个基因维度压缩成几十个综合变量主成分每个主成分都解释了数据方差的一部分。你可以想象成把一张几千维的地图压扁到几十维同时尽量保留了各个数据点之间的距离关系。主成分数量的选择我一般看三个证据# 方法一ElbowPlot直接看拐点 ElbowPlot(expr, ndims 40) # 方法二JackStraw做统计检验数据量大时比较耗时 expr - JackStraw(expr, num.replicate 100) expr - ScoreJackStraw(expr, dims 1:30) JackStrawPlot(expr, dims 1:30) # 方法三查看累计方差占比 std_devs - Stdev(object expr, reduction pca) var_explained - (std_devs^2) / sum(std_devs^2) cumulative_var - cumsum(var_explained)ElbowPlot是最直观的主成分数量在图中形成一个从陡峭到平缓的拐点通常取拐点附近的值。如果拐点不明显就参考累计方差占比是否达到70%到80%。一般经验是10到30个PC之间对大部分数据集都是合理的范围。这里有一个非常常见的反面案例有人不管三七二十一PC数量固定填20或者30。对于低复杂度样本30个PC里可能后面10个全在解释噪音这会导致聚类结果混乱分出来的几个cluster其实是同一个群体被强行拆碎对于高复杂度样本PC数量不够也会导致稀有细胞群无法分离。所以我强烈建议每个样本至少花个三分钟看一眼ElbowPlot再定PC数不要图省事。6.2 聚类分辨率UMAP图好不好看的幕后推手确定主成分后就是聚类FindClusters。这里有一个核心参数resolution分辨率。它决定了聚类的粒度——resolution越大产生的cluster数量越多。# 初步尝试多个分辨率 for (res in c(0.2, 0.5, 0.8, 1.2)) { expr - FindClusters(expr, resolution res, verbose FALSE) } # 用Idents查看不同分辨率下的cluster数量 table(exprmeta.data$RNA_snn_res.0.2) table(exprmeta.data$RNA_snn_res.0.5)怎么选分辨率我个人的原则是先粗后细结合细胞类型marker检验。先用低分辨率跑得到一个比较保守的分群比如10到15群看看major cell type能不能区分开然后提高分辨率把感兴趣的亚群再细分。最终分辨率的选择不是越大越好或者数字漂亮就好而是应该服务一个核心问题分出来的每一群细胞能否被已知的marker gene所解释。举个例子我在分析肿瘤组织样本时先用resolution 0.2跑出了大概10个cluster其中有一个cluster高表达CD3D、CD3E这些T细胞marker这就是一个T细胞的大群。但我注意到这个群里其实同时包含CD4和CD8阳性的细胞为了进一步区分我提高resolution到0.6T细胞群就会进一步分裂成CD4 T和CD8 T两个子群。用marker gene去判断cluster的生物学意义比单纯看cluster数量要有意义得多。6.3 运行UMAP并理解关键参数的影响聚类完成之后UMAP可视化登场。这是我个人认为整个流程中最具visual冲击力的环节也给不熟悉的读者解释一下本质。UMAP的目标是把高维空间中的数据点降维到二维平面上展示同时尽可能保持数据点之间的局部邻域关系。我习惯用RunUMAP做这一步expr - RunUMAP(expr, dims 1:20, umap.method uwot) DimPlot(expr, reduction umap, label TRUE)但UMAP有两个参数比很多人意识到的更容易影响出图效果n.neighbors决定UMAP在构图时考虑多大范围的邻居默认15。这个值越小局部结构越清晰但也会放大噪音值越大全局结构越明显但会模糊局部边界。如果样本量很大比如超过5万细胞建议适当调大到30到50否则小范围的随机波动可能导致图形碎掉。如果样本量很小几千个细胞默认值15就够用。min.dist决定点在低维空间中的最小距离默认0.3。这个值越小聚类群间的空隙就越大视觉上群与群之间分得越开值越大细胞群会挤在一起边界模糊。如果你发现UMAP图里各群挤成一团难以区分试着把min.dist调小到0.1或0.2。# 我的常用出图参数调整方案 expr - RunUMAP(expr, dims 1:20, n.neighbors 30, min.dist 0.2)注意一个常见误解UMAP图上的距离并不等于真实的细胞状态距离。UMAP为了保证局部邻域的呈现效果会故意拉近或压缩某些区域所以你看到的群之间的距离只能说明它们在高维空间中有差异不能说明差异大小。解释数据时一定要回到marker基因表达来支撑结论不能搞看图说话。7. 可视化进阶让UMAP图说清楚你想表达的话7.1 基础可视化的四板斧代码层面的功能都给你了viz步骤挺丰富的。最基本的四个可视化分别是DimPlot展示聚类结果标注每个cluster的编号适合俯瞰整体细胞分群结构。FeaturePlot展示一个或几个基因的表达量在UMAP上的分布是验证cluster生物学身份的金标准。VlnPlot展示某个或某几个marker基因在不同cluster中的表达分布比FeaturePlot更适合做定量比较。DoHeatmap展示指定marker基因在所有cluster中的表达热图适合做整体对比和文章插图。# 典型marker验证代码 FeaturePlot(expr, features c(CD3D, CD8A, CD4, CD68, NKG7, MS4A1)) VlnPlot(expr, features c(CD3D, CD8A, CD4), ncol 3)一张好的UMAP图配色要克制、标注要清晰、中文字体要处理得当。如果要做发表级的图微调DimPlot里的cols参数或者直接用ggplot2的扩展接口去做二次加工都是不错的选择。7.2 marker基因判断cluster身份的实际经验判断一个cluster是什么细胞类型需要至少2到3个正面marker和一个负面marker来佐证。只靠一个基因下结论容易翻车。我给自己定了一个检查清单T细胞CD3D、CD3E阳性同时CD68巨噬细胞marker阴性B细胞MS4A1CD20、CD79A阳性NK细胞NKG7、GNLY、KLRD1阳性CD3D为低表达或阴性单核/巨噬细胞CD68、LYZ阳性中性粒细胞FCGR3B、CSF3R阳性内皮细胞PECAM1CD31、VWF阳性成纤维细胞COL1A1、DCN阳性上皮细胞EPCAM、KRT19阳性实际分析中cluster身份的判定很少一蹴而就。更常见的是先挑出几个明确的核心marker锁定大方向再通过SetIdent把某个群单独拿出来做差异表达找到该群特异高表达的全部基因综合判断亚群身份。这个思路比marker撞库式地逐个尝试要高效得多。7.3 多组学整合时可视化的注意事项如果你做的是多组学数据比如同时有10X的转录组和表面蛋白Seurat的WNNWeighted Nearest Neighbor分析可以整合多模态信息。这个模式下运行UAMP的代码会变成expr - RunUMAP(expr, nn.method annoy, reduction weighted.nn)用WNN得到的UMAP图通常会比单独用RNA数据做出的图在细胞类型边界上更清晰因为额外的蛋白模态数据提供了更强的分辨信息。这个功能在我做CITE-seq数据时非常有用强烈推荐。8. 实战中的坑从数据异常到批次效应的一次性排查清单8.1 细胞数异常偏少或偏多如果你预期的细胞数是8000但Cell Ranger评估出来的只有2000问题很可能出在测序深度。你可以看一眼Mean Reads per Cell如果只有5000左右说明每个细胞的有效测序深度太低很多低表达的转录本没有被捕获到过滤后大量细胞被当成背景剔除了。这种情况下可以降级使用raw_feature_bc_matrix自己把过滤阈值放宽后再看。如果评估出的细胞数远高于预期比如预期5000但评估出15000先不要高兴这可能是因为有大量空液滴空载体背景被纳入了分析。这时候要重点看Median Genes per Cell如果这个值很低比如几百高度怀疑是背景信号过大。可以尝试把--expect-cells调低重新跑Cell Ranger或者在下游分析时用更严格的过滤阈值。8.2 线粒体基因比例整体偏高某个样本的percent.mt中位数超过20%甚至更高意味着大部分细胞都处于比较差的状态。这通常不是数据分析能解决的而是实验端的问题——可能是细胞分离过程太长细胞在制备过程中大量应激死亡。这种情况我会如实报告数据质量分析结果谨慎解读。但如果只是少数cluster线粒体比例高那可能是某种细胞类型的特点或某个细胞亚群处于应激状态这种情况可以考虑把相关cluster单独拿出来做差异分析看看它是否是有意义的生物学状态。8.3 批次效应怎么办更多时候你不是分析一个样本而是多个样本合并分析。如果合并后细胞不是按细胞类型聚集而是一个样本一团那就是明显的批次效应。处理方法我一般从两个层面考虑一是如果批次差异明确来自技术因素比如不同时间测序、不同芯片可以用harmony做批次校正。它的原理是迭代地把不同批次的数据在PCA空间中对齐校正后再去聚类效果通常立竿见影。library(harmony) expr - RunHarmony(expr, group.by.vars sample) expr - RunUMAP(expr, reduction harmony, dims 1:20)二是先不要急着校正先看看批次差异是不是由真实的生物学差异驱动比如肿瘤样本和癌旁样本免疫细胞组成本来就不同。这种情况直接校正反而会抹掉真实的病生理差异。我的建议是先合并、分群、用marker判断细胞类型看看类型分布如果同一类型的细胞还是明显分开成两群再做批次校正也不迟。8.4 UMAP图分不开怎么办如果你的UMAP图所有细胞挤成一团没有清晰的边界最常见的原因是主成分数量选得不对或者没有做充分的预处理。优先检查三件事是否做过ScaleData高变基因是否有意义PCA是否选对了维度按我踩坑经验80%的分不开都出在这三步。还有一个小概率但值得一查的情况如果你直接把counts矩阵喂给RunUMAP而没有经过PCA或没有用正确的reduction参数UMAP也会跑出一团糊。确认一下你是用reduction pca或harmony作为UMAP的输入而不是默认的原始数据空间。8.5 运行速度慢、内存不足的应对大样本比如10万个细胞以上跑Cell Ranger或者Seurat的某些步骤内存很容易不够。这里分享几个实战技巧Cell Ranger阶段如果基因组比对步骤跑得很慢检查一下--localcores是否真的生效了。另外确认机器有足够的临时磁盘空间用于STAR比对STAR比对会产生大量的临时文件磁盘满了会直接报错。Seurat阶段对大样本超过5万细胞我会采取两个策略第一把nfeatures降低到2000以下减少高变基因数量这会线性降低后续计算开销第二考虑使用Seurat的future包做并行计算同时把BPPARAM设置为SnowParam或MulticoreParam。ScaleData是最耗内存的步骤之一可以配合blocksize参数控制单次处理量避免一次性把所有数据加载到内存里。library(future) plan(multicore, workers 8) options(future.globals.maxSize 50 * 1024^3)9. 实操心得一些没人告诉你的小细节最后的最后分享几个我做了很多次之后才体会到的细节。这些内容不算核心技术但往往决定了你的分析过程顺畅不顺畅。Cell Ranger的--id参数既是输出目录名也会出现在后续所有输出文件的命名里。最好起一个能一眼认出样本信息的名字比如Tumor01_5prime_v2别偷懒用test或者output这种名字。不然等你要同时管理十几个样本的时候光看目录就头大。Seurat的SaveRDS和LoadRDS是保存分析中间结果的神器。每次完成一个阶段的分析就把Seurat对象保存一次。尤其是跑完SCTransform、harmony这种比较耗时的步骤后一定要保存。为什么因为这些步骤的随机种子如果不固定每次跑出的结果可能略有不同保存好结果可以保证你论文里的图和之后补跑的图一致。saveRDS(expr, file Tumor01_sct_harmony.rds) # 下次加载直接 expr - readRDS(Tumor01_sct_harmony.rds)关于随机种子UMAP和聚类都带有随机性。如果希望结果可复现请在运行之前设置set.seed(42)。一个很尴尬的场景是你之前做的聚类结果和marker基因图都对得上但重跑一次后cluster编号全变了比如原来的cluster 3现在是cluster 5如果没保存结果论文里图和代码对不上会折腾很久。最后一个建议整个分析过程中一定要养成记录的习惯。哪些参数用了什么值、为什么这么选、结果怎么看都记下来。我早期做的时候嫌麻烦不记录后期写方法部分时经常要回忆当初到底用的哪个版本的Cell Ranger、哪个版本的Seurat、filter参数是什么非常痛苦。项目多了之后我发现这些记录不仅是对自己负责更是对一篇论文可复现性的基本要求。现在我把每个样本的参数记录在一个Markdown文件里跟分析代码放在同一个目录下方便随时回溯。这个方法看似笨拙但长期坚持下来收益非常大。 单细胞转录组数据分析全流程：从FASTQ到UMAP的完整指南
做单细胞数据分析也有几年了从最早拿到10X数据手足无措到现在能一口气从FASTQ原始数据跑到UMAP出图中间踩过的坑不少但也沉淀了一套相对稳定、可以复用的完整流程。这篇文章我就把整套思路和实操细节写出来从最原始的FASTQ文件开始一步步走到最终的可视化结果希望能让刚入坑的朋友少走点弯路。我先说一下这篇文章适合谁看。如果你刚拿到10X Genomics的测序数据对Cell Ranger和Seurat只是听说过但不知道怎么串起来或者你已经跑通了流程但总觉得结果质量不太对劲、UMAP图片不理想那这篇文章应该很适合你。我默认你有基本的Linux命令行操作能力R语言有一定基础但即便是小白按着步骤一步步来也能跑出比较靠谱的结果。1. 整体流程拆解从FASTQ到UMAP中间到底发生了什么1.1 单细胞数据的特殊之处在讲具体操作之前我得先花点篇幅讲讲单细胞数据为什么不能像普通转录组那样直接处理。普通转录组RNA-seq拿到的是一堆细胞混在一起的平均表达值而10X Genomics技术解决的核心问题是给每一个细胞打上独特的身份证再给每一条转录本打上分子条码这样测序数据回来后我们就能从海量reads里还原出每个细胞里各有多少条mRNA。这里面有两个概念你必须彻底理解因为后面所有质控都围绕它们展开Cell Barcode细胞条码一段特异的序列用来标记mRNA来自哪一个细胞。UMIUnique Molecular Identifier唯一分子标识符一段随机序列用来区分同一条mRNA的重复扩增拷贝。简单打个比方细胞条码是门牌号告诉你这条数据是哪一户的UMI是购物小票就算同一件商品同一条mRNA被复制了好多份你也知道它本质上只买过一次。Cell Ranger这个工具干的事情就是把测序得到的海量reads按照这两重信息拆解成一张细胞 x 基因的数字矩阵。1.2 为什么选Cell Ranger而不是自己写脚本比对可能有人会问我能不能拿常规的STAR、hisat2去比对然后自己用featureCounts计数理论上可以但实践上非常不建议。原因有三个第一10X数据的文库结构特殊read 1的前16个碱基是细胞条码接下来10个碱基是UMI后面才是转录本序列。如果不用专门工具你得自己花很大力气去处理这三段信息的拆分而且容易出错。第二Cell Ranger内置了一个重要的纠错逻辑它会对细胞条码和UMI做校正。测序过程中会出现碱基错读如果一个细胞条码只差一个碱基Cell Ranger会认为它们来自同一个细胞这能显著降低数据噪音。自己写的脚本基本没有这个能力。第三Cell Ranger比对用的STAR比对器做了专门优化比对速度极快。我之前试过用常规STAR去比对一个8G左右的FASTQ文件跑了快十个小时而Cell Ranger只用不到三个小时差距非常大。从我的经验来看你就老老实实用官方工具跑通整个pipeline把精力省下来去做下游分析这才是正路。1.3 整体流程的五个核心节点从原始FASTQ到最终的UMAP可视化整个流程可以拆成五个关键节点节点输入核心工具主要产出数据质控与格式检查FASTQ原始文件fastqc / multiqcQC报告定量分析FASTQ 参考转录组Cell Ranger表达矩阵、BAM文件、质控指标数据过滤与预处理表达矩阵Seurat过滤后的Seurat对象降维聚类预处理后的矩阵Seurat harmony可选PCA、聚类结果可视化与解释聚类结果Seurat ggplot2UMAP图、marker基因图这个流程的主干是不变的每个环节都有它自己容易踩的坑。下面我逐个环节展开讲重点说那些我在实际项目中反复遇到、容易出问题的地方。2. 起点决定终点FASTQ数据准备与质控细节2.1 拿到数据的第一个动作检查完整性很多新手拿到数据就急着跑Cell Ranger结果跑到一半发现FASTQ文件不完整或者文件名不对白白浪费时间。我习惯的做法是拿到数据后先做三件事。第一确认FASTQ文件数量是否正确。10X的数据通常每个样本有四个文件R1、R2分别对应两条测序readI1、I2是index read。如果你用了双index的建库方式I1和I2都必须存在。少一个文件整个样本就作废了。第二检查文件大小是否合理。一般的10X单细胞数据一个样本的FASTQ文件加起来少则有几个GB多则几十个GB。如果你拿到手的文件只有几百MB非常可疑建议用zcat查看一下reads数是否达标。第三检查文件名是否符合Cell Ranger的命名规则。如果你用官方的cellranger mkfastq做过拆分文件名一般是SampleName_S1_L001_R1_001.fastq.gz这样的格式。如果是从其他渠道拿来的数据文件名可能比较乱建议先规范化命名免得后面启动cellranger count的时候识别不到样本。2.2 测序质量到底要不要先过滤FASTQ拿到手之后一个常见疑问是要不要先用fastp、trimmomatic之类的工具做一遍质量过滤再跑Cell Ranger我的答案是不需要甚至可以明确说不要这么做。原因在于Cell Ranger自身的pipeline中已经内置了read质量过滤的步骤而且它的过滤逻辑是专门针对单细胞文库结构设计的。你在外面再做一遍过滤反而可能把带有细胞条码信息的read误杀导致有效细胞数下降。但我不反对你做质控报告。用fastqc和multiqc看看整体测序质量的分布从大局上有帮助如果Q30比例过低比如低于80%你心里要有数后面检查比对率和有效barcode比例的时候重点留意。测序质量不好的数据不是简单的过滤能救回来的这一点要有心理准备。注意Cell Ranger自带的速度很快且专门针对单细胞测序的文库结构做了优化不要在外部再做额外的预处理。2.3 参考转录组的选择和准备参考基因组的选择直接影响比对结果和基因定量准确性这个环节值得仔细对待。人源样本我用的是10X官网提供的refdata-cellranger-GRCh38-2020-A。这个版本的参考转录组包含GENCODE v32的基因注释是官方构建好的里面已经包含pre-mRNA序列对于检测内含子区域很有帮助。为什么单细胞要用包含pre-mRNA的参考因为单细胞测序捕获的是细胞裂解出来的RNA大量mRNA还没来得及剪接完毕如果参考转录组里没有内含子序列这些未剪接转录本会被白白浪费掉直接导致检测到的基因数和UMI数偏低。具体在Cell Ranger里怎么构建参考转录组官方文档写得很清楚我在这里就不重复了。只说一个关键点优先用官方构建好的参考转录组而不是自己从GENCODE/Gencode下载GTF和FASTA再跑cellranger mkref。因为官方构建过程做了很多隐含的优化比如处理了基因ID和转录本版本的对应关系这些都是自己构建时容易出错的地方。人源之外的其他物种官方没有现成的参考那就只能自己跑cellranger mkref了。构建参考时要注意GTF文件与FASTA文件的版本一致性比如GRCm39就配GRCm39的GTF混搭不同版本很容易出现比对率极低、注释率极低的问题而且这类问题非常隐蔽不为新手常见。3. Cell Ranger实战一次跑通的6个关键设置3.1 cellranger count命令的完整参数解析假设你的工作目录是/data/scRNAFASTQ文件在/data/scRNA/FASTQ目录下样本名为Tumor01我们要跑的命令大概长这样cd /data/scRNA cellranger count \ --idTumor01_hg38 \ --transcriptome/data/ref/refdata-cellranger-GRCh38-2020-A \ --fastqs/data/scRNA/FASTQ \ --sampleTumor01 \ --expect-cells5000 \ --localcores16 \ --localmem64这里我逐个参数讲一下我的设置习惯。--expect-cells这个参数比较关键。它告诉Cell Ranger这个样本大概有多少个细胞会影响Cell Ranger判断有效barcode的阈值。我一般会先看看建库时候的目标捕获数通常目标在8000左右我就填--expect-cells8000。如果填得太保守比如实际有10000个细胞但只填了3000那高表达的细胞没什么问题但低表达的细胞群体可能会被当成背景噪音丢掉。反之如果填得比实际高很多背景的barcode会被大量保留下来下游分析会非常痛苦。拿不准的时候宁可填得稍微高一点也不要低估因为下游过滤时低质量细胞还能通过QC步骤去掉但如果细胞压根没被识别到那就彻底丢了。--localcores和--localmem是控制计算资源分配的。这里有个容易出问题的点--localmem的单位是GB不是MB。我有段时间没注意填了--localmem32000结果Cell Ranger把它当成32TB内存来用直接跑到OOM被系统kill了。正确做法是看机器空闲内存比如机器有128GB内存就给Cell Ranger分配64到96GB。核心数同理不要贪心留出几个核给系统的其他进程免得互相抢资源。3.2 比对环节发生了什么为什么它这么重要cellranger count内部的比对环节是基于STAR比对器的。这一步干的事是把read 2上测到的转录本序列比对到参考基因组上而read 1上的细胞条码和UMI基本上直接被保留和使用不做太多处理。这里有个非常重要的细节Cell Ranger默认允许一个转录本比对到基因组的多处位置。这个设计是深思熟虑的——因为参考基因组里有很多来自同一个基因家族的重复序列如果一个read只比对到一个位置而不管其他位置很多基因的定量就会偏低。Cell Ranger在比对完成后会根据UMI数去判断这条read到底更可能来自哪个位置这个逻辑比简单粗暴的唯一比对要合理得多。比对率是考察数据质量的核心指标之一。正常情况下人源样本的比对率应该在85%到95%之间。如果你的比对率低于80%先检查参考转录组是否存在版本问题其次考虑样本是否有其他物种污染。3.3 看完web_summary.html才算跑完跑完cellranger count并不是万事大吉。每个样本跑完后都会在outs目录下生成一个web_summary.html文件这个文件包含了全套质控指标是我必看的内容。我一般会重点看几个数Estimated Number of Cells评估出的细胞数和预期差距过大就要警惕。Mean Reads per Cell每个细胞的平均reads数人源样本要做到20000以上比较理想但如果你测序深度不高只要后续分析能出结果也不是绝对的。Median Genes per Cell每个细胞检测到的中位基因数人源样本一般要在2000到5000之间如果低于1000说明数据质量堪忧。Q30 Bases in Barcode细胞条码区域的Q30比例这个数值低说明建库环节可能有问题。Valid Barcodes有效barcode比例低于80%说明可能有index错配或者建库污染。需要说明的是这些指标没有一个绝对的标准需要结合你实际的实验设计、物种、建库方案来综合判断。我之前有一个样本中位基因数只有1200左右一开始以为质量很差后来发现是做了神经元分选后的样本——神经元这种细胞本身就是低复杂度RNA表达这个数值是正常的。3.4 Cell Ranger的结果目录里到底有什么跑完Cell Ranger后outs目录下会生成多个文件。对下游分析我主要用到这几个文件路径内容用途outs/filtered_feature_bc_matrix/过滤后的细胞-基因表达矩阵Seurat分析的主要输入outs/raw_feature_bc_matrix/未过滤的全部barcode矩阵用于自定义过滤不常用outs/analysis/聚类、差异表达等初步分析快速预览结果outs/bam/比对后的BAM文件重新定量或可视化时用下游分析的核心输入基本就是filtered_feature_bc_matrix这个目录它包含三个文件matrix.mtx.gz、features.tsv.gz和barcodes.tsv.gz。这三个文件是一个标准的稀疏矩阵存储格式Seurat可以直接读取。你要理解这三样东西的匹配关系barcodes.tsv是细胞ID列表features.tsv是基因ID列表matrix.mtx.gz是它们的行和列对应的计数矩阵。4. 进入Seurat表达矩阵读取与质量控制在细节里4.1 用Read10X快速读取矩阵数据进入R环境后读取数据这一步很简单但有些细节值得注意。我习惯这样写library(Seurat) library(dplyr) library(ggplot2) data_dir - outs/filtered_feature_bc_matrix data - Read10X(data.dir data_dir) # 这里特别注意如果有多个样本需要合并先各自创建对象再merge expr - CreateSeuratObject( counts data, project Tumor01, min.cells 3, min.features 200 )min.cells 3的意思是一个基因至少在3个细胞中有表达才被保留这个参数能过滤掉那些只在极少数细胞里出现的低质量基因min.features 200的意思是一个细胞至少检测到200个基因才被保留这个阈值过滤掉了那些可能没有细胞核或RNA含量极低的空液滴。特别强调一下读取路径的问题Read10X函数的data.dir参数要求的是包含三个文件的目录路径不是matrix.mtx.gz的文件路径。很多人第一次用的时候把路径写到了具体的文件层级结果报错找不到数据这个坑不算大但很烦。如果你有多个样本要合并分析我建议先分别创建Seurat对象再用merge()函数合并。在合并后加group信息后面分析批次效应时好区分。4.2 线粒体基因比例、核糖体基因比例和双细胞识别Seurat对象创建好之后第一个核心步骤是计算质控指标。这里最常用的三个指标是nCount_RNA每个细胞的UMI总数、nFeature_RNA每个细胞检测到的基因数和percent.mt线粒体基因占所有UMI的比例。# 计算线粒体基因比例人源样本用MT-前缀小鼠用mt-前缀 expr[[percent.mt]] - PercentageFeatureSet(expr, pattern ^MT-) # 核糖体基因比例可以过滤掉某些极端情况 expr[[percent.ribo]] - PercentageFeatureSet(expr, pattern ^RP[SL]) # 查看QC指标的分布 VlnPlot(expr, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)这三个指标背后的生物学意义你要理解nFeature_RNA和nCount_RNA过低细胞可能已经破裂RNA流出去了或者这个液滴是空的。nCount_RNA很高但nFeature_RNA很低可能是某个低复杂度细胞类型比如红细胞但也可能是文库污染。percent.mt过高细胞处于应激或即将凋亡的状态。因为胞质mRNA在破裂时会被降解而线粒体mRNA由于包裹在线粒体内部反而相对稳定。所以一个濒临死亡的细胞会表现出高线粒体基因比例。但是做好功课某些组织类型比如脑组织里的少突胶质细胞天然线粒体比例偏高一刀切过滤容易误杀。这时候多看几组数据结合解剖位置判断。关于双细胞doublets简单说一下两个细胞被包在同一个液滴里被测序这种情况下nCount和nFeature通常都会异常高。你可以先用简单的阈值卡掉nFeature超过一定值的细胞比如人源样本6000但如果需要更严格的过滤建议用DoubletFinder这样的专门工具。后面我会在常见问题里再展开。4.3 过滤阈值怎么定看分布而不是拍脑袋过滤阈值是很多新手特别纠结的地方。我见过有人直接拿教程里的nFeature_RNA 500 nFeature_RNA 5000 percent.mt 20套用到自己的数据上结果效果很差。正确的做法是先画图看分布再根据分布形态确定阈值。# 画QC指标分布 plot1 - FeatureScatter(expr, feature1 nCount_RNA, feature2 percent.mt) plot2 - FeatureScatter(expr, feature1 nCount_RNA, feature2 nFeature_RNA) plot1 plot2比如某个样本的nFeature_RNA分布主峰在2000到4000之间然后有一个小尾巴延伸到8000以上那这个8000以上就很可能是双细胞阈值应该设在8000附近或略低一些。再比如percent.mt在主峰附近如果集中在10%以下那阈值设在15%或20%就比较合理。但如果你的样本有某种细胞类型天然线粒体高你看到分布的主峰直接在15%到20%之间那阈值设在25%甚至30%也合理。你得记住一个核心理念过滤器只是去除明显的异常值而不是做精细生物学划分。真正决定细胞类型界限的是后续的聚类分析不要在过滤环节下手太狠。5. 归一化与高变基因选择让数据说话的前提5.1 LogNormalize和SCTransform怎么选过滤完之后下一步是归一化。我用Seurat的时候面对两个选项传统的LogNormalize和推荐的SCTransform。LogNormalize的逻辑很直观每个细胞的UMI总数标准化到相同的深度默认10000然后做log变换。这个方法的优点是简单、快、结果也稳定缺点是它没有完全考虑转录本的技术噪音——高表达基因在细胞间的波动天然比低表达基因大LogNormalize没法很好地处理这个问题。SCTransform做了一件更聪明的事它对每个基因建立了一个广义线性模型把UMI数和建库深度之间的关系拟合出来然后把深度效应直接去掉。这样归一化后的数据能更好地反映真实的生物学差异而不是技术操作差异。我的经验是如果你的数据量不大比如小于2万个细胞直接用SCTransform替代整个NormalizeData和ScaleData步骤通常结果更干净。如果数据量很大超过5万个细胞SCTransform的计算时间会明显拉长这时候也可以先跑LogNormalize走一遍流程看看聚类结果再说。还有一个折中方案用LogNormalize做初步的聚类拿到结果后对感兴趣的细胞亚群用SCTransform再精修一遍这样速度和精度都能兼顾。# 方法一传统LogNormalize流程 expr - NormalizeData(expr, normalization.method LogNormalize, scale.factor 10000) expr - FindVariableFeatures(expr, selection.method vst, nfeatures 2000) expr - ScaleData(expr) # 方法二SCTransform流程 expr - SCTransform(expr, vars.to.regress percent.mt, verbose FALSE)SCTransform里vars.to.regress参数可以顺带回归掉技术变量比如线粒体基因比例。这里要注意一个原则不要随便回归掉你认为不重要的生物学变量。很多人看到percent.mt高第一反应是把它回归掉但线粒体压力所反映的细胞状态本身可能是有生物学意义的把它回归掉相当于抹掉了这部分信息。有争议的情况下先不做回归看看聚类结果再说。5.2 高变基因2000个还是更多的选择逻辑FindVariableFeatures这一步选出的高变基因HVGs决定了你后续PCA分析能看到什麼样的坐标系。默认是找2000个高变基因这对于大部分样本已经够用。但也有例外如果你研究的样本组织复杂度非常高比如脑、肝脏这种有多种细胞类型且每种都有特异表达基因的组织可以考虑提高到3000或4000个让更多稀有细胞类型的特异基因入选。这个逻辑不难理解高变基因就是PCA降维时用来构建主成分的特征集如果你的稀有细胞类型的关键特异基因不在这个特征集里那这个细胞群在聚类时就可能被其他细胞群的特征淹没难以分离。选基因之后ScaleData这一步会把基因表达值标准化为z-score均值为0方差为1并且默认对每个基因做一次线性回归去掉细胞测序深度的影响。这一步在SCTransform流程里已经被包含进去了所以如果你用了SCTransform就不要再重复调用ScaleData否则会出问题。# 传统流程需要单独ScaleData expr - ScaleData(expr, vars.to.regress nCount_RNA)需要留意的是如果数据中某些基因的表达量全为0比如过滤后有些基因存在于矩阵中但所有细胞都不表达ScaleData会提示错误或给出警告这个问题出现频率不低。解决办法是在ScaleData之前先去掉全0基因Seurat现在的版本通常会帮你处理但如果你用很老的版本还是要手动检查。6. 降维聚类的完整实操从PCA到UMAP的关键决策6.1 PCA分析主成分数量不该拍脑袋决定PCA降维是整个流程中最技术的一步因为它的结果直接影响后续聚类质量。核心问题是保留多少个主成分PC用于后续分析这个决定不能拍脑袋要靠证据。首先明确PCA在干什么我们有了高变基因矩阵但维度还是很高几千个基因直接做聚类计算量巨大且容易受噪音干扰。PCA把几千个基因维度压缩成几十个综合变量主成分每个主成分都解释了数据方差的一部分。你可以想象成把一张几千维的地图压扁到几十维同时尽量保留了各个数据点之间的距离关系。主成分数量的选择我一般看三个证据# 方法一ElbowPlot直接看拐点 ElbowPlot(expr, ndims 40) # 方法二JackStraw做统计检验数据量大时比较耗时 expr - JackStraw(expr, num.replicate 100) expr - ScoreJackStraw(expr, dims 1:30) JackStrawPlot(expr, dims 1:30) # 方法三查看累计方差占比 std_devs - Stdev(object expr, reduction pca) var_explained - (std_devs^2) / sum(std_devs^2) cumulative_var - cumsum(var_explained)ElbowPlot是最直观的主成分数量在图中形成一个从陡峭到平缓的拐点通常取拐点附近的值。如果拐点不明显就参考累计方差占比是否达到70%到80%。一般经验是10到30个PC之间对大部分数据集都是合理的范围。这里有一个非常常见的反面案例有人不管三七二十一PC数量固定填20或者30。对于低复杂度样本30个PC里可能后面10个全在解释噪音这会导致聚类结果混乱分出来的几个cluster其实是同一个群体被强行拆碎对于高复杂度样本PC数量不够也会导致稀有细胞群无法分离。所以我强烈建议每个样本至少花个三分钟看一眼ElbowPlot再定PC数不要图省事。6.2 聚类分辨率UMAP图好不好看的幕后推手确定主成分后就是聚类FindClusters。这里有一个核心参数resolution分辨率。它决定了聚类的粒度——resolution越大产生的cluster数量越多。# 初步尝试多个分辨率 for (res in c(0.2, 0.5, 0.8, 1.2)) { expr - FindClusters(expr, resolution res, verbose FALSE) } # 用Idents查看不同分辨率下的cluster数量 table(exprmeta.data$RNA_snn_res.0.2) table(exprmeta.data$RNA_snn_res.0.5)怎么选分辨率我个人的原则是先粗后细结合细胞类型marker检验。先用低分辨率跑得到一个比较保守的分群比如10到15群看看major cell type能不能区分开然后提高分辨率把感兴趣的亚群再细分。最终分辨率的选择不是越大越好或者数字漂亮就好而是应该服务一个核心问题分出来的每一群细胞能否被已知的marker gene所解释。举个例子我在分析肿瘤组织样本时先用resolution 0.2跑出了大概10个cluster其中有一个cluster高表达CD3D、CD3E这些T细胞marker这就是一个T细胞的大群。但我注意到这个群里其实同时包含CD4和CD8阳性的细胞为了进一步区分我提高resolution到0.6T细胞群就会进一步分裂成CD4 T和CD8 T两个子群。用marker gene去判断cluster的生物学意义比单纯看cluster数量要有意义得多。6.3 运行UMAP并理解关键参数的影响聚类完成之后UMAP可视化登场。这是我个人认为整个流程中最具visual冲击力的环节也给不熟悉的读者解释一下本质。UMAP的目标是把高维空间中的数据点降维到二维平面上展示同时尽可能保持数据点之间的局部邻域关系。我习惯用RunUMAP做这一步expr - RunUMAP(expr, dims 1:20, umap.method uwot) DimPlot(expr, reduction umap, label TRUE)但UMAP有两个参数比很多人意识到的更容易影响出图效果n.neighbors决定UMAP在构图时考虑多大范围的邻居默认15。这个值越小局部结构越清晰但也会放大噪音值越大全局结构越明显但会模糊局部边界。如果样本量很大比如超过5万细胞建议适当调大到30到50否则小范围的随机波动可能导致图形碎掉。如果样本量很小几千个细胞默认值15就够用。min.dist决定点在低维空间中的最小距离默认0.3。这个值越小聚类群间的空隙就越大视觉上群与群之间分得越开值越大细胞群会挤在一起边界模糊。如果你发现UMAP图里各群挤成一团难以区分试着把min.dist调小到0.1或0.2。# 我的常用出图参数调整方案 expr - RunUMAP(expr, dims 1:20, n.neighbors 30, min.dist 0.2)注意一个常见误解UMAP图上的距离并不等于真实的细胞状态距离。UMAP为了保证局部邻域的呈现效果会故意拉近或压缩某些区域所以你看到的群之间的距离只能说明它们在高维空间中有差异不能说明差异大小。解释数据时一定要回到marker基因表达来支撑结论不能搞看图说话。7. 可视化进阶让UMAP图说清楚你想表达的话7.1 基础可视化的四板斧代码层面的功能都给你了viz步骤挺丰富的。最基本的四个可视化分别是DimPlot展示聚类结果标注每个cluster的编号适合俯瞰整体细胞分群结构。FeaturePlot展示一个或几个基因的表达量在UMAP上的分布是验证cluster生物学身份的金标准。VlnPlot展示某个或某几个marker基因在不同cluster中的表达分布比FeaturePlot更适合做定量比较。DoHeatmap展示指定marker基因在所有cluster中的表达热图适合做整体对比和文章插图。# 典型marker验证代码 FeaturePlot(expr, features c(CD3D, CD8A, CD4, CD68, NKG7, MS4A1)) VlnPlot(expr, features c(CD3D, CD8A, CD4), ncol 3)一张好的UMAP图配色要克制、标注要清晰、中文字体要处理得当。如果要做发表级的图微调DimPlot里的cols参数或者直接用ggplot2的扩展接口去做二次加工都是不错的选择。7.2 marker基因判断cluster身份的实际经验判断一个cluster是什么细胞类型需要至少2到3个正面marker和一个负面marker来佐证。只靠一个基因下结论容易翻车。我给自己定了一个检查清单T细胞CD3D、CD3E阳性同时CD68巨噬细胞marker阴性B细胞MS4A1CD20、CD79A阳性NK细胞NKG7、GNLY、KLRD1阳性CD3D为低表达或阴性单核/巨噬细胞CD68、LYZ阳性中性粒细胞FCGR3B、CSF3R阳性内皮细胞PECAM1CD31、VWF阳性成纤维细胞COL1A1、DCN阳性上皮细胞EPCAM、KRT19阳性实际分析中cluster身份的判定很少一蹴而就。更常见的是先挑出几个明确的核心marker锁定大方向再通过SetIdent把某个群单独拿出来做差异表达找到该群特异高表达的全部基因综合判断亚群身份。这个思路比marker撞库式地逐个尝试要高效得多。7.3 多组学整合时可视化的注意事项如果你做的是多组学数据比如同时有10X的转录组和表面蛋白Seurat的WNNWeighted Nearest Neighbor分析可以整合多模态信息。这个模式下运行UAMP的代码会变成expr - RunUMAP(expr, nn.method annoy, reduction weighted.nn)用WNN得到的UMAP图通常会比单独用RNA数据做出的图在细胞类型边界上更清晰因为额外的蛋白模态数据提供了更强的分辨信息。这个功能在我做CITE-seq数据时非常有用强烈推荐。8. 实战中的坑从数据异常到批次效应的一次性排查清单8.1 细胞数异常偏少或偏多如果你预期的细胞数是8000但Cell Ranger评估出来的只有2000问题很可能出在测序深度。你可以看一眼Mean Reads per Cell如果只有5000左右说明每个细胞的有效测序深度太低很多低表达的转录本没有被捕获到过滤后大量细胞被当成背景剔除了。这种情况下可以降级使用raw_feature_bc_matrix自己把过滤阈值放宽后再看。如果评估出的细胞数远高于预期比如预期5000但评估出15000先不要高兴这可能是因为有大量空液滴空载体背景被纳入了分析。这时候要重点看Median Genes per Cell如果这个值很低比如几百高度怀疑是背景信号过大。可以尝试把--expect-cells调低重新跑Cell Ranger或者在下游分析时用更严格的过滤阈值。8.2 线粒体基因比例整体偏高某个样本的percent.mt中位数超过20%甚至更高意味着大部分细胞都处于比较差的状态。这通常不是数据分析能解决的而是实验端的问题——可能是细胞分离过程太长细胞在制备过程中大量应激死亡。这种情况我会如实报告数据质量分析结果谨慎解读。但如果只是少数cluster线粒体比例高那可能是某种细胞类型的特点或某个细胞亚群处于应激状态这种情况可以考虑把相关cluster单独拿出来做差异分析看看它是否是有意义的生物学状态。8.3 批次效应怎么办更多时候你不是分析一个样本而是多个样本合并分析。如果合并后细胞不是按细胞类型聚集而是一个样本一团那就是明显的批次效应。处理方法我一般从两个层面考虑一是如果批次差异明确来自技术因素比如不同时间测序、不同芯片可以用harmony做批次校正。它的原理是迭代地把不同批次的数据在PCA空间中对齐校正后再去聚类效果通常立竿见影。library(harmony) expr - RunHarmony(expr, group.by.vars sample) expr - RunUMAP(expr, reduction harmony, dims 1:20)二是先不要急着校正先看看批次差异是不是由真实的生物学差异驱动比如肿瘤样本和癌旁样本免疫细胞组成本来就不同。这种情况直接校正反而会抹掉真实的病生理差异。我的建议是先合并、分群、用marker判断细胞类型看看类型分布如果同一类型的细胞还是明显分开成两群再做批次校正也不迟。8.4 UMAP图分不开怎么办如果你的UMAP图所有细胞挤成一团没有清晰的边界最常见的原因是主成分数量选得不对或者没有做充分的预处理。优先检查三件事是否做过ScaleData高变基因是否有意义PCA是否选对了维度按我踩坑经验80%的分不开都出在这三步。还有一个小概率但值得一查的情况如果你直接把counts矩阵喂给RunUMAP而没有经过PCA或没有用正确的reduction参数UMAP也会跑出一团糊。确认一下你是用reduction pca或harmony作为UMAP的输入而不是默认的原始数据空间。8.5 运行速度慢、内存不足的应对大样本比如10万个细胞以上跑Cell Ranger或者Seurat的某些步骤内存很容易不够。这里分享几个实战技巧Cell Ranger阶段如果基因组比对步骤跑得很慢检查一下--localcores是否真的生效了。另外确认机器有足够的临时磁盘空间用于STAR比对STAR比对会产生大量的临时文件磁盘满了会直接报错。Seurat阶段对大样本超过5万细胞我会采取两个策略第一把nfeatures降低到2000以下减少高变基因数量这会线性降低后续计算开销第二考虑使用Seurat的future包做并行计算同时把BPPARAM设置为SnowParam或MulticoreParam。ScaleData是最耗内存的步骤之一可以配合blocksize参数控制单次处理量避免一次性把所有数据加载到内存里。library(future) plan(multicore, workers 8) options(future.globals.maxSize 50 * 1024^3)9. 实操心得一些没人告诉你的小细节最后的最后分享几个我做了很多次之后才体会到的细节。这些内容不算核心技术但往往决定了你的分析过程顺畅不顺畅。Cell Ranger的--id参数既是输出目录名也会出现在后续所有输出文件的命名里。最好起一个能一眼认出样本信息的名字比如Tumor01_5prime_v2别偷懒用test或者output这种名字。不然等你要同时管理十几个样本的时候光看目录就头大。Seurat的SaveRDS和LoadRDS是保存分析中间结果的神器。每次完成一个阶段的分析就把Seurat对象保存一次。尤其是跑完SCTransform、harmony这种比较耗时的步骤后一定要保存。为什么因为这些步骤的随机种子如果不固定每次跑出的结果可能略有不同保存好结果可以保证你论文里的图和之后补跑的图一致。saveRDS(expr, file Tumor01_sct_harmony.rds) # 下次加载直接 expr - readRDS(Tumor01_sct_harmony.rds)关于随机种子UMAP和聚类都带有随机性。如果希望结果可复现请在运行之前设置set.seed(42)。一个很尴尬的场景是你之前做的聚类结果和marker基因图都对得上但重跑一次后cluster编号全变了比如原来的cluster 3现在是cluster 5如果没保存结果论文里图和代码对不上会折腾很久。最后一个建议整个分析过程中一定要养成记录的习惯。哪些参数用了什么值、为什么这么选、结果怎么看都记下来。我早期做的时候嫌麻烦不记录后期写方法部分时经常要回忆当初到底用的哪个版本的Cell Ranger、哪个版本的Seurat、filter参数是什么非常痛苦。项目多了之后我发现这些记录不仅是对自己负责更是对一篇论文可复现性的基本要求。现在我把每个样本的参数记录在一个Markdown文件里跟分析代码放在同一个目录下方便随时回溯。这个方法看似笨拙但长期坚持下来收益非常大。