做了几年单细胞RNA测序数据分析,最深的感受是:这个领域看起来门槛高,但真正卡人的地方往往不是算法本身,而是对每一步为什么要这么做的理解。很多人拿到10X Genomics下机的数据就开始跑Seurat默认流程,跑完发现聚类结果一塌糊涂,也不知道该调哪里。这篇文章把我从基础原理到高级应用的完整思路写下来,包含我实际跑项目时的参数选择、踩坑记录和排查经验,希望能帮你少走弯路。
1. 单细胞测序到底是什么,为什么传统转录组不够用了
1.1 从“群体平均”到“单个细胞”
传统转录组测序(bulk RNA-seq)测的是组织块里成千上万细胞的平均表达量。一个肿瘤组织样本里的癌细胞、T细胞、巨噬细胞、成纤维细胞混在一起,最终拿到的是一个“混合口味的果汁配方”——你能知道里面大概有什么成分,但不知道每种水果原本的味道。
单细胞RNA测序把每一个细胞的转录组单独测出来。以10X Genomics平台为例,单个细胞被油滴包裹,带上独特的barcode标签,后续测序数据里每个reads都能溯源到它来自哪个细胞。这样你就可以回答一些传统转录组回答不了的问题:这个组织里到底有哪些细胞类型?每种类型占多少比例?某个亚群是否发生了状态转变?
我经常用一个类比来理解单细胞数据:bulk RNA-seq像一碗蔬菜汤打成泥,你知道里面有胡萝卜味也有芹菜味,但分不清比例;单细胞测序则是把每根菜叶单独夹出来看,你能数清楚这碗汤里到底有几种菜、各占多少、哪些菜已经开始蔫了。
1.2 一个完整的单细胞项目长什么样
在实际项目里,从拿到下机数据到最终发表图,大致要经过这样一个流程:
| 环节 | 核心任务 | 常用工具 |
|---|---|---|
| 数据产出 | 碱基序列 → 基因表达矩阵 | Cell Ranger / STARsolo |
| 数据质控 | 过滤低质量细胞和基因 | Seurat / Scanpy |
| 去双细胞 | 去除一个液滴里的多个细胞 | DoubletFinder / Scrublet |
| 归一化 | 消除测序深度差异 | LogNormalize / SCTransform |
| 降维 | PCA + UMAP / tSNE | Seurat / scanpy |
| 聚类 | 无监督识别细胞群 | Louvain / Leiden |
| 注释 | 标记基因鉴定细胞类型 | SingleR / CellTypist + 人工核对 |
| 下游分析 | 差异表达、轨迹、通讯等 | Seurat / Monocle / CellChat |
这个流程看着简单,但每一步都有不少需要决策和处理的地方。很多人拿到数据就直接跳到聚类,忽略了质控和归一化阶段的问题,后续分析就会很被动。预处理阶段的每一个决定都会传导到最终结果,所以前期的理解特别重要。
2. 入门前的准备工作:数据理解和环境搭建
2.1 表达矩阵里的三种数据类型
开始分析前,你先要把手里的数据搞清楚。Cell Ranger输出的filtered_feature_bc_matrix目录下有三个文件:
- matrix.mtx:稀疏矩阵,三列分别代表基因ID、细胞Barcode、表达量计数
- features.tsv:基因注释信息,包含基因ID、基因名、数据类型
- barcodes.tsv:细胞barcode列表
这三个文件本质上描述了一张二维矩阵:行是基因(大约2万个),列是细胞(几千到几万不等),矩阵里的值代表每个基因在每个细胞中检测到的UMI(Unique Molecular Identifier)计数。
UMI是10X平台的一个设计巧思。每个mRNA分子在反转录时会被加上一段随机的UMI序列,测序后通过UMI去重,可以让相同序列的reads归并为一个转录本分子,大大降低PCR扩增带来的偏差。理解了这一点,你就知道为什么单细胞表达量通常用UMI count而不是reads count来做定量。
2.2 计算环境和工具链怎么选
单细胞分析主流的两个框架是R语言生态下的Seurat和Python生态下的Scanpy。我个人的建议是:入门用Seurat,因为参考资料多,社区活跃,画图方便;进阶可以学Scanpy,处理超大数据的性能更好,机器学习的生态也更完整。
如果是个人电脑做练习,建议优先用RStudio配合Seurat。10万细胞以下的数据集,内存16GB也能勉强跑,但做UMAP和聚类时可能卡顿,建议64GB内存或者用服务器。做真实的临床样本项目,一个样本通常有5万个细胞以上,我一般直接上服务器跑,避免本地内存不足浪费时间。
安装Seurat时容易忽视的是依赖包的版本问题。直接用install.packages("Seurat")往往会被安装到比较旧的版本,建议用:
install.packages("Seurat") packageVersion("Seurat") # 确认版本号如果你用的是Seurat v5,需要额外注意数据对象的存储方式从原来的data变成了layers结构,很多网上的旧代码会报错。比如GetAssayData(object, slot = "data")在v5中变成了LayerData(object, assay = "RNA", layer = "data"),这些细节在跑旧教程时很常见,遇到报错不妨先看版本号。
3. 核心细节解析与实操要点:从矩阵到可分析对象
3.1 创建Seurat对象——不是套一个函数那么简单
读入数据创建Seurat对象,看起来就是一行代码,但里面的逻辑值得讲清楚:
library(Seurat) data_dir <- "filtered_feature_bc_matrix" counts <- Read10X(data_dir) obj <- CreateSeuratObject(counts = counts, project = "demo", min.cells = 3, min.features = 200)关键参数是min.cells和min.features。min.cells = 3表示这个基因至少在3个细胞中有表达才保留,主要是过滤掉那些只在极少数细胞中零星出现的基因,这些基因大概率是测序噪音。min.features = 200表示一个细胞至少检测到200个基因才保留,这个阈值可以过滤掉大部分破裂的、空的液滴。
你的项目如果样本很小,比如只有几千个细胞,min.cells可以适当降低,否则大量低表达但真实的基因会被过滤掉。但也不要设得太低,否则下游分析会明显变慢且噪音高。实际项目中我通常看一个“先过滤后统计”的策略:先用比较宽松的标准建对象,然后根据质控指标再做精细过滤。
3.2 质控指标解读——线粒体比例为什么这么重要
质控是单细胞分析中最关键的一步,核心看三个指标:
- 每个细胞的基因数(nFeature_RNA):代表复杂度和细胞活力
- 每个细胞的UMI数(nCount_RNA):代表测序深度
- 线粒体基因比例(percent.mt):代表细胞健康状况
线粒体基因比例的底层逻辑是:如果细胞膜破损,细胞质的mRNA会从破裂处泄漏出去,而线粒体因为相对稳定,它的转录本会相对富集。所以percent.mt高的细胞通常说明这个细胞已经受损或者正在凋亡。
计算线粒体比例在Seurat里很容易,但有个容易出错的细节——物种不同,线粒体基因前缀不同:
# 人: MT- # 小鼠: mt- # 斑马鱼: mt- (注意大小写规则不完全相同) obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")一个小坑:你在网页数据库下载的公共数据可能已经做了质控,也可能没做。拿到数据后先用VlnPlot和FeatureScatter快速看一下分布,再决定是否需要进一步过滤。公共数据集的质控标准未必适合你的下游分析目标,比如你要做罕见细胞亚群,就不能用太严格的阈值把低表达基因的细胞全滤掉了。
3.3 质控阈值到底怎么定——先看分布再定门坎
很多教程直接告诉你nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10,但实际项目中这个硬编码不总是适用。我通常的做法是:
VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")先看分布图。如果大部分细胞的基因数集中在中间,两侧有明显的离群尾巴,那就在尾巴拐弯的地方设阈值。比如大部分外周血单核细胞的nFeature_RNA在800到3500之间,一个肿瘤组织的细胞可能因为类型多样,范围更宽。
线粒体阈值在10%到20%之间调整是常态。一些特殊的组织样本,比如肝脏、肌肉因为代谢旺盛,baseline的线粒体比例天然会高一些,这时10%可能过于严格,可以放宽到20%。关键是阈值设定要有可视化依据,不要拍脑袋。我一般会把过滤前后的细胞数、基因数变化记下来,形成项目日志,方便后续写方法部分。
3.4 去双细胞(Doublet)——最容易被忽略却影响深远的一步
双细胞是两个细胞被一个液滴捕获,同一个barcode下混着两个不同细胞的转录本。双细胞占的比例通常在0.4%到8%之间,在细胞数量多的样本中绝对数量不小。如果没有去除,它们会形成独立的中间状态聚类群,污染你后续的注释。
常见的去双细胞工具包括DoubletFinder(R包,基于人工构造的双细胞模拟)和Scrublet(Python,同理的原理)。
DoubletFinder的原理可以简化理解为:人为把已知的细胞两两合并成模拟双细胞,然后训练一个分类器去识别真实数据中与模拟双细胞相似的细胞。使用时有几个参数值得注意:
library(DoubletFinder) ## 先完成Seurat标准流程的前几步 obj <- NormalizeData(obj) %>% FindVariableFeatures(nfeatures = 2000) %>% ScaleData() %>% RunPCA() ## 根据细胞数估算Doublet rate # 10X官方:1000 cells约0.4%,10000 cells约7.6%实际经验是:DoubletFinder需要先做PCA,且pN = 0.25是默认值,pK建议用paramSweep来寻找最优值,不要直接套默认参数。我在一个8万细胞的数据集上跑默认参数得到1万多个“假双细胞”,调整pK后降到合理的5000左右。这个步骤没有绝对正确参数,建议跑完看被标记细胞的marker表达——如果标记的细胞里有很多同时表达两种明确不共存的marker,比如T细胞和B细胞的marker,那么去双大概率是合理的。
3.5 归一化和SCTransform——消除测序深度影响的两种思路
归一化的核心是解决“同一个基因在两个细胞里表达量不同,可能是因为两个细胞本身的RNA总量不同,而不是真正的表达差异”。
Seurat默认的LogNormalize方法:
log1p( count / total_counts_per_cell * 10000 )相当于把每个细胞的UMI总数标准化到1万,再取log。这个方法是有效的,但它不能完全消除测序深度与技术变异的相关性。
更好用的替代是SCTransform,它构建了一个正则化负二项回归模型来去除技术噪音,效果确实更好。如果样本量不大、计算资源充足,我通常直接用SCTransform代替NormalizeData + FindVariableFeatures + ScaleData三步:
obj <- SCTransform(obj, vars.to.regress = "percent.mt", verbose = FALSE)但要注意,SCTransform生成的SCT assay在后续找marker基因时,和RNA assay的用法不太一样,FindAllMarkers在SCT assay下建议用logfc.threshold适当调低,因为SCT的表达值是残差而非实际count,logFC的解读方式有差异。实际项目中我习惯用RNA assay跑差异表达,用SCT assay跑聚类。
4. 降维、聚类与可视化——从计算到图形的每一步决策
4.1 高变基因筛选和PCA——为什么不用全部2万个基因
做完归一化后,Seurat会默认寻找高变基因(variable features),默认选2000个。这一步的意义在于:大部分基因在不同细胞间的表达变化很小,用这些基因做下游分析只会引入噪音;而真正决定细胞身份差异的,往往是那批在不同细胞中表达差异大的基因。
PCA降维的本质是把高维表达矩阵(几千个高变基因)转换为一组线性不相关的“主成分”,其中前几个主成分往往能解释数据中的主要差异来源。在单细胞项目中,PCA还有提高信噪比的作用——后续聚类的输入不是原始的2万个基因维度,而是那几十个主成分。
PCA要跑多少个主成分?Seurat里可以用ElbowPlot看拐点:
ElbowPlot(obj, ndims = 50)但说实话,拐点法只是一个参考。我现在的习惯是用JackStraw的p值显著性结合拐点图,同时参考自己的经验值——一般10万细胞以下的数据集,dims设置在20到30之间比较常见。主成分数量太少会丢失罕见的细胞亚群,太多又会引入噪音。聚类和UMAP图对dims的选择比较敏感,建议同一个数据集分别跑dims = 15和dims = 30,如果聚类结果差别很大,说明你的数据结构和参数选择都需要重新审视。
4.2 聚类算法背后的直觉——Louvain与Leiden
单细胞聚类现在最主流的是图聚类算法,Louvain在很长时间是默认选项,而Leiden在2023年后的Seurat版本里也逐渐成为替代选择。它们的基本流程是:先构建细胞间KNN图(每细胞连接最近的K个邻居),然后在这个图上寻找模块化的分区。
Leiden比Louvain的改进在于解决了已知的社区检测“断连”问题,在单细胞数据的测试中大多表现更好。实际操作上,你只需要关心分辨率参数resolution——它控制聚类的“颗粒度”。同一组细胞,resolution越高,分出来的类就越多。
选分辨率没有绝对标准,我的做法是:
- 从
resolution = 0.5起跑,结合marker热图看每个群的生物学意义; - 如果发现某些群里混多种细胞类型(比如一个群里T细胞、B细胞marker都高),适当调大分辨率;
- 如果已知的某个细胞类型被拆散了,调低分辨率再观察。
这两天我跑一个肿瘤样本数据,用了resolution = 0.8,分出了27个cluster,其中有5个cluster都表达T细胞marker(CD3D、CD3E),但是单独看这几个cluster,有的高表达CD8A,有的高表达CD4,有的高表达GZMB(细胞毒性或耗竭),这说明这个精度下T细胞亚群可以被进一步区分,是有意义的拆分,而不是过聚类。
4.3 UMAP vs tSNE——选哪个、参数怎么调
UMAP几乎是现在单细胞可视化的默认选择。相比tSNE,UMAP在保持全局结构方面更好,而且计算速度快得多,更重要的是细胞间的距离有一定意义——UMAP中靠得近的细胞不一定相似,但离得远的细胞大概率不相似。这个特性对于解释下游聚类结果很有帮助。
RunUMAP的关键参数是n.neighbors和min.dist:
n.neighbors越大,UMAP越关注全局结构,得到的图越“松散”,适合看大的发育轨迹min.dist越小,Cluster间的簇越紧凑,适合区分关系相近的细胞亚群
我常用的是min.dist = 0.3,这个数值能让不同来源的细胞在图上分开,又不会完全把同一类型拆散。如果你发现UMAP图过度拥挤,每个群都挤成一团看不清楚结构,试试调大min.dist到0.5。
另外强烈建议在注释之前不要反复用不同随机种子跑UMAP——这会消耗你大量时间去“碰运气”找一个视觉效果最好的图。UMAP的分布有随机性,小范围的重排是正常的,只要聚类结果是稳定的,不必执着于某个具体布局。
4.4 细胞类型注释——用marker还是用自动注释工具
注释是单细胞分析中最耗费精力的环节,也是最能体现经验的部分。我的工作流分两步:
第一步:自动注释粗筛。常用工具包括SingleR、CellTypist、scCATCH等。对PBMC样本,SingleR搭配celldex里的参考数据效果不错;对组织样本,CellTypist的免疫细胞和上皮细胞模型表现更好。
library(SingleR) ref <- celldex::HumanPrimaryCellAtlasData() pred <- SingleR(test = GetAssayData(obj, assay = "RNA", layer = "data"), ref = ref, labels = ref$label.main) obj$singleR_labels <- pred$labels但注意,自动注释的本质是基于已知参考数据的模式匹配,遇到罕见细胞、疾病状态下的异常细胞或者新细胞亚型,参考数据库可能完全没有覆盖。所以自动注释只能用来“缩小范围”,不能作为最终答案。
第二步:marker gene人工核对。这一步是不能跳过的。我把常用的marker整理成表格方便随时查阅:
| 细胞类型 | 经典marker(人) |
|---|---|
| T细胞 | CD3D, CD3E |
| CD8+ T细胞 | CD8A, CD8B |
| CD4+ T细胞 | CD4, IL7R |
| NK细胞 | NKG7, GNLY, KLRD1 |
| B细胞 | MS4A1, CD79A |
| 单核细胞 | CD14, LYZ |
| 树突状细胞 | FCER1A, CST3 |
| 内皮细胞 | PECAM1, VWF |
| 成纤维细胞 | COL1A1, DCN |
| 上皮细胞 | EPCAM, KRT8, KRT18 |
用FeaturePlot和DotPlot在cluster水平检查marker表达,确认每个cluster的注释是否符合预期。实际项目中常见的问题是某些cluster表达多个细胞类型的marker——这种情况要么是双细胞污染残留,要么是某种中间态,要么是注释分辨率不够需要细分。我的原则是:宁可多分出一个小群,也不要强行合并表达模式不同的细胞。
5. 从聚类走向故事:差异表达、轨迹分析与细胞通讯
5.1 差异表达和marker基因鉴定——统计检验的坑要知道几个
FindAllMarkers是注释后的常用操作,它比较每个cluster与其他所有cluster的基因表达差异。默认使用的Wilcoxon秩和检验在这里有一些值得注意的局限:单细胞数据存在大量零值(dropout事件),而且细胞与细胞在同一个样本中并不完全独立,这会导致检验统计量偏乐观——p值往往极小,但生物学意义未必同样突出。
所以我更看重两个维度:效应量(average log2 fold change)和表达比例(pct.1 / pct.2)。一个基因在某个cluster里90%的细胞表达、其他簇只有5%表达,比一个基因虽然p值极低但两边的表达比例都是30%更值得关注。
markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) top_markers <- markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)一个常见的坑是:多组间比较时建议用FindConservedMarkers分别对不同样本组(比如疾病组vs对照组)做分组内的marker检测,再取交集或合并结果,这样可以避免某些“marker”完全是由某个样本主导的表达差异。
5.2 轨迹分析(拟时序分析)——什么时候该用Monocle,什么时候不该用
轨迹分析的目的是推断细胞群之间的发育、分化或状态转变关系。常用工具包括Monocle3、Slingshot、scVelo(RNA velocity)等。
但这条分析需要克制。很多人在数据里看到一种细胞类型的大群和另一种类型的亚群,就想着跑轨迹做分化推断,但轨迹推断有两个前提:1)数据确实存在一个连续的生物学过程;2)你对起点细胞有合理的生物学判断或证据。
我的建议是:先用UMAP观察是否存在明显连续过渡的结构,如果多种细胞类型边界清晰,没有中间态,那跑轨迹很可能只会得到人为拼接的假轨迹。如果UMAP里细胞沿一条路径连续分布,且已知marker沿路径渐变表达,那轨迹分析才是有意义的。
Monocle3在Seurat对象上的衔接比较顺畅:
library(monocle3) cds <- as.cell_data_set(obj) cds <- cluster_cells(cds) cds <- learn_graph(cds) cds <- order_cells(cds, root_pr_nodes = get_earliest_principal_node(cds))但要注意,as.cell_data_set()后基因名和细胞ID的顺序、重复名字等都可能报错,通常需要重新构建CDS对象并手动传入gene_metadata和cell_metadata。我遇到过几次从Seurat转换时基因重复导致的报错,通常是Ensembl ID和symbol混用导致的,保持一致能避免不少麻烦。
5.3 细胞通讯分析——CellChat和CellPhoneDB怎么选
细胞通讯分析利用配体-受体数据库推测不同细胞类型之间的互作关系。CellChat(R包)使用简单、结果丰富;CellPhoneDB(Python)数据库质量高;NicheNet则可以做配体-靶基因调控推断。
实用建议是:先用CellChat跑一个整体图,看大方向上哪种细胞与哪种细胞的互作最强、涉及哪些信号通路,然后聚焦到与项目背景相关的信号通路做细看。CellChat在运行前需要把表达矩阵归一化到CellChat自身的数据结构:
library(CellChat) data.input <- GetAssayData(obj, assay = "RNA", layer = "data") meta <- data.frame(labels = obj$celltype, row.names = colnames(obj)) cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels") cellchat <- addCellChatDB(cellchat, db = "CellChatDB.human")一个常见问题是:细胞数太多时CellChat运行会非常慢。可以对每个细胞类型先取一部分细胞做subsampling,保证每个类型500到1000个细胞即可稳定估计互作强度,不需要把全部细胞放进去。
5.4 批量效应整合——多个样本合并时的必修课
真实项目大部分都包含多个样本(不同患者、不同处理组),直接合并分析会带来明显的批次效应:聚类结果很容易按样本来源而非生物学差异分开。这时需要用整合算法。
常见选择:
- Harmony:速度快,在大型数据集上表现稳,是目前最常用的选择
- Seurat CCA整合:适合不同细胞类型差异很大的情况(比如跨组织比较)
- scVI:深度学习框架,适合大数据集,但需要配置Python环境
Harmony在Seurat里用起来很方便:
obj <- NormalizeData(obj) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA() obj <- RunHarmony(obj, group.by.vars = "sample") obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30)整合后要检查两个问题:1)不同样本是否在UMAP上混合分布(对比整合前后的图);2)已知的细胞类型marker是否仍然清晰分离。如果某个细胞群在整合后消失了,可能说明整合过度,需要调整lambda参数或改用别的算法。
6. 常见问题与排查技巧实录
6.1 质控后细胞数少得可怜怎么办
先确认是不是原始数据就很少。Cell Ranger的web summary里会报告estimated number of cells,如果这个数本来就只有预期的一半,问题出在文库构建或测序环节。如果Cell Ranger报告的细胞数正常,但过滤后少了很多,多半是质控标准太严格了。这时候可以逐项放宽并观察细胞数的变化,比如把nFeature_RNA的下限从200降到100,percent.mt从10放宽到15,但要警惕这可能是样本质量本身不佳的信号——遇到这种情况最好重新评估是否该补做样本。
6.2 聚类结果里某个cluster的marker不明确
这大概是单细胞分析里最常遇到的情况。处理思路有三条路径:1)调大分辨率,看这个cluster是不是可以细分;2)跑FindMarkers看这个cluster前30差异表达基因是什么——有时候marker是你之前没查过的基因,需要去文献里确认;3)如果无论如何都找不到明确的细胞身份,可以考虑标注为”Unknown”或”TBD”,后续用更下游的分析(如拷贝数变异推断)来辅助判断。
某个cluster的marker表达不高也可能是混合了多种细胞类型,此时重新检查是否存在双细胞残留、是否整合过度抹掉了差异,或者直接提高分辨率后再看。
6.3 数据量太大,内存溢出怎么办
我曾在单台服务器上处理过一个40万细胞的数据集,Seurat的默认流程跑到ScaleData就崩了。解决办法是按顺序尝试:1)用Plan参数改为"multicore"或"multisession"并行计算;2)用上BPCells(Seurat v5支持的磁盘存储格式)把表达矩阵放在磁盘而不是内存里;3)先用Harmony整合后降采样到10万细胞做初探,确认主要细胞类型后,在需要特别关注的亚群上重新聚类全量细胞。对于真实的项目场景,很多时候你不需要在完整数据集上调参,先跑通子集、确认流程、再全量运行是更高效的策略。
6.4 常见报错速查
| 报错信息 | 常见原因 | 解决方案 |
|---|---|---|
Cannot find a "counts" layer | 对象没有RNA assay的counts层 | 检查数据是否用Read10X正确读入;SCTransform后RNA assay的counts还在 |
Error in validObject(.Object) | 元数据或assay不兼容 | 检查细胞barcode是否有重复;用make.unique(colnames(obj)) |
Cannot find variable features | 没有跑FindVariableFeatures | 先执行FindVariableFeatures再跑RunPCA |
Error: matrix has non-numeric entry | 表达矩阵里有NA或字符 | 用as.matrix检查并转换数据类型 |
6.5 一张高质量可视化的沟通要点
单细胞可视化不仅是给自己看的,也是给合作者或审稿人看的。我习惯产出的标准图包括:UMAP分群图(按cluster着色)、细胞类型注释后的UMAP图、关键marker的FeaturePlot叠加图、不同样本/条件在UMAP上的分布对比图、marker热图。这些图组合在一起,基本可以讲清楚一个单细胞项目的故事。
绘制UMAP时统一主题和配色,保存时用高分辨率(ggsavedpi=300,或pdf格式保存矢量图)。文字说明里标注清楚用了什么软件和版本,是方法部分的基本要求。
7. 高级应用的进阶思路——从一个“合格分析”走向一个有深度的生物学故事
7.1 细胞比例分析和差异丰度检验
细胞类型注释完成后,比较不同样本组之间每种细胞类型的比例变化,是单细胞分析从“描述是什么”走向“某某处理改变了什么”的第一步。但细胞比例分析有个统计陷阱:每个样本中细胞比例不是独立的,一个类型比例升高会拉低其他类型的比例,构成成分数据(compositional data)问题。用普通t检验那是不合适的,建议至少使用propeller(speckle包)这类专门针对单细胞比例差异的方法,或者在样本量足够时用广义线性模型。
7.2 CNV推断——从转录组看基因组变异
对于肿瘤单细胞数据,inferCNV或CopyKAT可以利用转录组的基因表达模式推断大范围的染色体拷贝数变异(CNV)。这个分析可以帮助你区分肿瘤细胞和正常细胞——肿瘤细胞往往会有大量的CNV,正常的免疫细胞一般没有。
inferCNV的运行提醒:需要提供明确的正常细胞作为参考(通常用T细胞或髓系细胞),且运行时间较长。还有一个常见的误区:肿瘤细胞和正常细胞在细胞周期上存在差异,高增殖的肿瘤细胞会被错误推断出染色体范围表达异常,建议先做细胞周期评分(CellCycleScoring)并考虑回归掉周期差异再跑inferCNV。
7.3 转录因子调控网络推断
SCENIC是单细胞转录因子调控网络分析的主流工具,它的核心思想是利用转录因子结合位点信息推断每个细胞中的“调节子”(regulon)活性。SCENIC在R里也有pyscenic的Python版本。这个分析比较吃时间和内存,通常只在关注某条调控通路时使用,比如在T细胞耗竭或巨噬细胞极化研究中,想看哪些转录因子驱动了状态转变。
7.4 空间转录组数据整合
空间转录组(如10X Visium)提供了一个组织切片的转录组空间分布信息。把单细胞数据的细胞类型映射到空间位置上是如今很热门的分析思路,常用工具包括Seurat的FindTransferAnchors、Cell2location等。
这类分析弥补了单细胞数据“消化道”的信息——单细胞把组织打碎成单个细胞,空间转录组则保留位置信息,两者结合可以回答“某种细胞在组织里分布在哪里、旁边是什么细胞”一类的问题。分析方法还不像单细胞流程那么固定,但大方向是有参考价值的。
8. 写在最后的一些实操体会
做了这么多单细胞项目,最大的体会是:分析流程本身并不是壁垒,网上教程太丰富了,真正拉开差距的是对生物学问题的理解和对数据质量的判断力。你有再漂亮的UMAP图,如果质控时metadata搞错了样本分组,或者注释时把T细胞标成了B细胞,后面所有的下游分析都会白费。
最后分享一个我自己坚持的工作习惯:每跑一个项目,我都会在项目根目录建一个analysis_notes.md,记录下每一步用的参数、为什么这样改、中间产物的细胞数变化、以及临时的想法。单细胞分析是一个高度迭代的过程,很多时候你回头看一个月前的选择记不清当初的原因,而这偏偏是写文章方法部分、回复审稿人意见时最需要的东西。这个习惯已经帮我避免了好几次返工。
如果你正卡在某个环节报错,或者对某个参数选择拿不准,可以照着这篇文章的步骤前四个章节排查一遍,大部分问题都在数据处理前几步能找到源头。