1. 肿瘤细胞推断这件事为什么不能只靠“猜”做单细胞测序分析的人几乎都会在某个阶段碰到同一个问题拿到一堆细胞注释完免疫细胞、基质细胞剩下的那一坨“未知细胞”到底是什么是不是肿瘤细胞如果是又该怎么证明给审稿人看我之前在处理一批结直肠癌样本的时候就卡在这道坎上。CellRanger跑完、Seurat标准流程走完、免疫细胞亚群也标好了然后盯着UMAP上一团既不像T细胞也不像B细胞、既没有典型髓系标记也没有内皮/成纤维标记的细胞群发呆——这团细胞到底是不是恶性克隆如果直接按“拷贝数变异推断”下结论又担心被质疑如果只靠已知标记基因很多实体瘤的肿瘤细胞根本不表达什么特异性的表面标记。后来把整个推断逻辑捋清楚之后才发现这件事本质上不是一个“找标记”的问题而是一个“整合证据链”的问题。肿瘤细胞之所以难认是因为它没有一个通用的、像CD3之于T细胞那样稳定的标记它的恶性身份来自基因组层面的异常而不是转录层面的某个分子。所以单细胞数据里推断肿瘤细胞的核心思路从根本上是绕过“已知标记”转向“基因组状态”的证据拷贝数变异也就是CNV。这篇文章我会完整复盘我当时用的推断流程从CNS级别的经典方法到实际操作中的坑再到最后如何把结果呈现到论文里。整个过程不依赖任何商业软件全部基于R和几个开源包适合正在做肿瘤单细胞数据分析、尤其是实体瘤样本的人参考。2. 为什么inferCNV是当前最稳的起点2.1 从原理上理解它在算什么inferCNV是目前单细胞肿瘤研究里使用最广泛的CNV推断工具至少在“第一步筛查”这个层面它依然是首选。它做的事情简单说就是利用整个转录组里的基因表达量去反推每个细胞在染色体水平上是否存在大范围的拷贝数异常。它的逻辑核心在于这样一个假设如果一条染色体臂或者某个染色体区段发生了整段的拷贝数缺失或扩增那么这个区段内所有基因的表达量会呈现一致性的下调或上调。这种信号非常微弱在单个基因的水平上完全看不出来但当你把整条染色体臂上的几十上百个基因的滑动平均表达量放在一起看时这种“区域化的表达偏差”就会浮现出来。可以类比成看一张被压缩过的图片单独看一个像素点噪点很多但把这一片像素整体做模糊和降采样之后原本被掩盖的大块色块就能看出来了。inferCNV做的就是这种“空间平滑”的活。2.2 输入数据该怎么准备很多人第一次跑inferCNV就挂在输入数据的格式上。这里我把当时实际可用的流程写出来。inferCNV需要两个核心输入单细胞表达矩阵可以是count矩阵也可以是标准化后的矩阵基因在染色体上的位置注释文件如果你的数据是从10X的CellRanger输出来的表达矩阵是filtered_feature_bc_matrix目录下的matrix.mtx直接用Seurat读进来然后转成普通矩阵就行。我当时是先用Seurat做了全套的标准降维聚类确认了主要细胞大类之后才把表达矩阵导出来给inferCNV用。基因注释文件如果不想自己写可以直接用inferCNV官方仓库里自带的gencode_v21_gen_pos.complete.txt人源的版本也比较新。如果你用的是小鼠数据也有对应的gencode_vM23版本。这个文件长这样chr1 11869 14409 DDX11L1每一行记录一个基因的染色体、起始位置、终止位置和基因名用tab分隔。注意inferCNV要求基因名这一列不要带ENSEMBLID最好是Symbol格式。如果你从Seurat里拿到的矩阵行名是ENSEMBL ID记得先转换。2.3 参考细胞群的选择直接决定结果质量选参考细胞reference cells是inferCNV最关键的一步但也是实际操作中最多人踩坑的地方。选择参考细胞的目的是提供一个“正常的基线”让那些拷贝数正常的细胞在结果里呈现为平坦的、无明显区域波动的状态而肿瘤细胞则相对地显现出CNV峰谷。什么样的细胞适合做参考两个要求一是确实来自非恶性细胞群二是数量上不能太少。我当时用的是同一份样本里注释出的T细胞和髓系细胞作为参考。为什么不用上皮细胞因为在很多上皮来源的肿瘤里所谓“正常上皮”和“肿瘤上皮”很难从转录组上完全切开用它们做参考会把肿瘤信号稀释掉。还需要注意一点参考细胞的数量最好在几十到几百的量级太少的话估计基线不稳太多则可能出现部分参考细胞自身含有弱CNV信号的情况导致肿瘤细胞的信号被相对压低。一般我建议每种参考细胞类型选50~200个细胞作为代表。2.4 跑通一次inferCNV的最小代码示例这里给出一个我认为最精简又能直接跑的inferCNV流程。前提是你已经用Seurat处理过数据对象叫seu里面已经有一个cell_type的metadata列且你定义好了哪些是参考细胞。library(infercnv) # 导出表达矩阵行是基因列是细胞 expr_mat - as.matrix(seuassays$RNAcounts) # 生成细胞分组注释 cell_annot - data.frame( cell colnames(seu), type as.character(seu$cell_type) ) # 选择参考细胞类型 ref_group - c(T_cell, Myeloid) # 运行inferCNV infercnv_obj - CreateInfercnvObject( raw_counts_matrix expr_mat, annotations_file cell_annot, gene_order_file gencode_v21_gen_pos.complete.txt, ref_group_names ref_group ) infercnv_obj - infercnv::run( infercnv_obj, cutoff 0.1, out_dir infercnv_output, cluster_by_groups TRUE, denoise TRUE, HMM TRUE )cutoff这个参数的值需要解释一下。对10X的UMI count数据官方推荐cutoff0.1对Smart-seq2这类全长转录组数据推荐cutoff1。原因在于UMI本身的定量方式会让低表达基因的计数稀疏性更强用一个较小的cutoff来过滤掉非常低表达的基因避免它们对滑动平均造成干扰。跑完之后inferCNV会输出一个热图infercnv.png和一系列HMM分析的结果文件。但请注意输出的热图颜色是“红蓝”风格红色代表扩增蓝色代表缺失看的时候别搞反。很多初学者第一次看到一片红蓝会觉得眼花缭乱其实只需要看每列细胞在染色体上的连续波动模式就能初步判断哪些细胞的CNV信号更集中。3. 拿到inferCNV结果后下一步不是直接下结论3.1 别把“热图颜色深浅”当证据很多人跑完inferCNV看一眼热图觉得某群细胞颜色偏红或者偏蓝就写进文章说“这群是肿瘤细胞”。这种事我见过太多最后被审稿人问倒的例子也不少。原因很简单热图是可视化的结果不是统计检验的结果不同样本、不同批次、不同参考细胞热图的色标范围都不一样。A样本里看起来“红得发紫”的细胞放到B样本里可能只是普通波动。所以拿到inferCNV结果之后第一件事不是看热图而是看HMM的输出。HMM隐马尔可夫模型分析是inferCNV里用来判断“哪些区段确实存在CNV”的统计模块会给出每个细胞在每个区域上的状态0表示正常1表示一个拷贝的缺失2表示拷贝数翻倍3表示更高倍数的扩增。这些状态才是可以量化和后续分析的基础。当时我处理样本时最有效的做法是把每个细胞的HMM状态统计一下计算每个细胞的“总CNV长度”或“CNV评分”。一群细胞如果整体上在多个染色体臂存在大范围的HMM非0状态那它大概率是恶性的相反如果你的参考细胞在HMM结果里也出现了大段的非0区域那就说明参考群选得不对或者denoise参数没调好需要回头重新跑。这个“计算每个细胞的CNV评分”的步骤在inferCNV输出结果里其实没有直接给一个现成的表需要自己写几行代码从HMM的结果矩阵里汇总。这是整个流程里我认为最值得投入时间做的一步。3.2 从HMM状态到每个细胞的CNV评分inferCNV跑完后在输出目录里会生成一个HMM_CNV_predictions.HMMi6.leiden.hmm_mode-subclusters.Pred_CNV_call.txt文件。这个文件里每一行是一个细胞每一列是一段基因组区域值就是前面说的0/1/2/3状态。我当时写了一段很简单的代码来汇总每个细胞的CNV负荷hmm - read.delim(HMM_CNV_predictions.HMMi6.leiden.hmm_mode-subclusters.Pred_CNV_call.txt, row.names 1, check.names FALSE) # 统计非0区域的数量和总长度 cnv_burden - data.frame( cell rownames(hmm), n_events apply(hmm, 1, function(x) sum(x ! 0)), total_gain apply(hmm, 1, function(x) sum(x[x 1])), total_loss apply(hmm, 1, function(x) sum(x[x 0])) ) # 简单评分非0区域占比 cnv_burden$score - cnv_burden$n_events / ncol(hmm)这里n_events是每个细胞发生非正常状态的区域数量score是区域占比。用这个分数把细胞按群分组画箱线图你会看到肿瘤细胞群普遍显著高于参考细胞群。这个分数是可以放到论文里的定量依据。3.3 每个样本单独跑不要混着跑这是我在实际分析中踩过的比较大的一次坑。最开始为了省时间想把自己手上的几个样本合并成一个Seurat对象然后用一套参考细胞跑一次inferCNV。结果发现不同样本的测序深度、细胞组成、参考细胞比例差异很大合并后跑出来的CNV模式混乱很难解释。原因在于inferCNV的CNV推断是一种相对推断——它比较的是每个细胞相对于参考细胞群体的表达偏离。不同样本的技术批次效应会干扰这种相对比较。正确做法是每个样本单独跑一次inferCNV各自选各自的参考细胞。这样出来的结果既干净也更容易解释。后面做跨样本比较时再基于每个细胞各自的CNV评分进行整合。4. 从CNV推“肿瘤身份”时最容易犯的三个错误4.1 把全部上皮细胞都当作肿瘤细胞如果你研究的是癌种是上皮来源肺癌、结直肠癌、乳腺癌等那么上皮细胞这个大类里通常同时混有正常上皮和恶性上皮。inferCNV跑完不能一看到“上皮细胞”这个注释就直接说它们是肿瘤细胞。必须结合CNV结果细分。这一步的实际操作是把inferCNV的CNV评分映射回Seurat的UMAP上用FeaturePlot看分数的分布。如果上皮细胞亚群内部出现明显的分数差异通常意味着有些上皮细胞是正常的另一些是恶性的。这时候需要根据CNV评分重新对上皮细胞进行亚聚类和注释。我当时在结直肠样本里就遇到过这种情况同一群上皮细胞一部分细胞的CNV评分接近参考细胞另一部分则明显偏高而且这种差异和它们在UMAP上的位置几乎完美对应。如果当时不加区分地全部注释为“肿瘤上皮”后续差异表达分析的结果会混入大量正常上皮的信号导致假阳性。4.2 把“高CNV评分”直接等同于“增殖能力强”这是个概念混淆的问题。CNV是基因组层面的静态变异增殖能力是转录层面的动态状态两者并不是一回事。一个肿瘤细胞可能携带大量CNV但处于静息状态G0期而一个正常组织干细胞可能没有CNV但增殖非常活跃。所以如果你在文章里需要讨论“肿瘤干性”或者“增殖亚群”不要用CNV评分来代替应该单独用细胞周期评分或者干性相关基因集打分。CNV的作用是“确立身份”不是“描述状态”。4.3 忽略低质量细胞对inferCNV的干扰inferCNV对输入数据中的低质量细胞非常敏感。具体表现是低质量细胞低UMI数、高线粒体比例在表达矩阵中表现出整体表达水平的随机波动这种波动经过滑动平均后看起来很像CNV信号。我通常会在跑inferCNV之前做一轮更严格的质控过滤比Seurat标准流程更严格一些。具体指标参考指标常规分析inferCNV前推荐nFeature_RNA 200 500nCount_RNA 500 1000percent.mt 20% 10%双细胞预测可选必须过滤如果细胞数量充足宁可多滤掉一些可疑细胞也要保证进入inferCNV的每个细胞都是高质量的。我自己的经验是这一步过滤做得越干净后面HMM结果的假阳性就越少。5. 除了inferCNV还有哪些方法可以做交叉验证5.1 CopyKAT用“阴阳”思路做快速判断CopyKAT是另一个比较常用的工具它和inferCNV的思路类似也是利用基因表达来推断CNV但有一个明显不同的地方它内置了一种叫“阴阳相关”的统计方法用来区分“二倍体细胞”和“非整倍体细胞”。运行速度比inferCNV快很多因为不需要做HMM迭代输出也比较友好直接给每个细胞一个标签aneuploid非整倍体或diploid二倍体。我的习惯是先用CopyKAT快速筛一遍再用inferCNV做精细确认。两个工具如果结论一致证据链就非常扎实如果不一致就需要单独排查。CopyKAT在R里用法也很简单library(CopyKAT) copykat - copykat( rawmat expr_mat, id.type S, ngene.chr 5, sam.name sample1, plot.genes TRUE, genome hg20 )解析copykat返回结果里的prediction列即可其中aneuploid就是候选肿瘤细胞。5.2 用已知的肿瘤特异突变做锚点如果你的样本有配套的肿瘤组织WES或WGS数据那么最直接的验证方式就是把单细胞数据里检测到的关键驱动基因突变和WES里的突变对比。常见做法是在单细胞表达矩阵里直接查看已知突变基因比如TP53、KRAS、PIK3CA等在候选肿瘤细胞群中的表达情况或者更严格地用varTrix之类的工具直接在单细胞BAM里检测突变。这个方法需要的前提是你的单细胞测序包含了全转录组信息10X 5或3都可以且已知的驱动突变位点在你测序覆盖的区域内。在结直肠癌里APC、TP53、KRAS、SMAD4是四大高频突变基因如果候选肿瘤细胞群在这些基因上存在表达异常或突变支持说服力会大幅上升。不过要注意的是RNA表达水平和基因组突变并不是一一对应的关系——有些突变会导致无义突变介导的mRNA降解反而在表达量上表现为下降。所以这个方法更适合做“正支持”而不是“负排除”。5.3 染色体整倍性变化的另一个视角inferCNV之外的“染色体拷贝数谱”最后一个交叉验证方法是直接从表达矩阵里提取“染色体表达程序”来做无监督聚类。这个做法不需要任何CNV推断工具只需要用所有基因的平均表达量按染色体位置重新排列就能得到每个细胞的一条“染色体表达曲线”。大体步骤如下把表达矩阵按基因在染色体上的位置排序计算每条染色体上每个基因的平均表达量或用滑动窗口平滑把每个细胞的染色体表达曲线做相关性聚类如果某个亚群的细胞在特定染色体区域整体上呈现出与参考细胞不同的表达模式那它很可能存在CNV。这个方法虽然粗糙但可以作为inferCNV结果之外的一个独立证据尤其在数据量大、计算资源有限的情况下非常好用。6. 如何把CNV推断结果写进论文图表设计与描述规范6.1 一张“堆叠式CNV热图”的标准展示逻辑论文里最常见的展示方式就是一张“细胞列 × 染色体位置行”的CNV热图上面用颜色表示拷贝数增益/缺失。这张图的读者是审稿人和同行他们的阅读习惯是把注意力集中在配色、细胞分群边界和染色体位置标注这三点上。我当时做图用的是inferCNV默认输出的热图但在放进论文前做了两点额外处理把细胞按“参考细胞”“肿瘤细胞亚群1”“肿瘤细胞亚群2”等分组顺序重新排好列使每个组的CNV模式一眼能看出来在热图上方加了一个简单的聚类树或分组色条标明每个细胞属于哪个注释类型这样做的好处是不需要读者自己去数细胞颜色分区直接说明了“肿瘤细胞群在多个染色体臂上存在一致的CNV模式”。6.2 用“CNV评分箱线图”量化恶性程度热图适合给人看直观印象但缺少统计感。我会在正文或补充材料里加一张箱线图横轴是细胞注释类型纵轴是前面计算的CNV评分每个点是一个细胞。然后做一个简单的统计检验Wilcoxon秩和检验比较肿瘤细胞群与非肿瘤细胞群的CNV评分是否有显著差异。这张图的优点在于它把“某群细胞有没有CNV”变成了“某群细胞相对参考细胞的CNV负荷显著更高”这个表述是非常标准的论文写作语言。我当时写的就是类似这样一句话通过inferCNV分析恶性上皮细胞亚群显示出显著的染色体拷贝数变异负荷与正常上皮细胞及免疫细胞参考群相比差异具有统计学意义Wilcoxon rank-sum test, P 0.001。这种写法不会引起审稿人“你怎么证明这是肿瘤细胞”的质疑因为结论是建立在可重复的、量化的数据之上。6.3 结合病理和公共数据库做最后的“确认”最后一个建议很多人会忽略如果你有条件拿到同一例患者的病理切片或者公共数据库里同类肿瘤的CNV信息尽量做一次交叉比对。比如GISTIC2.0输出的TCGA同类肿瘤的显著CNV区域列表和你从单细胞里推断出的CNV区域如果高度重叠那就不是巧合是生物学信息。我当时把自己样本中推断出的高频扩增区域和TCGA-COAD队列做了对比发现不少重叠的染色体区段这直接让审稿人对“肿瘤细胞身份判定”这一环节打了高分。7. 当数据质量不完美时如何调整分析策略7.1 参考细胞数量太少怎么办有些样本特别是肿瘤纯度特别高的样本里面正常细胞——T细胞、B细胞、髓系细胞——的数量都非常少可能只有几十个甚至十几个。这种情况下inferCNV的参考基线很容易不稳。一个可行的应急方案是从另一个检测平台或同一个患者的外周血单细胞数据里提取免疫细胞作为参考。但前提是这些免疫细胞的测序深度、平台类型要和肿瘤组织基本一致否则批次效应会干扰基线估计。另一个方案是降低对“精确CNV区域”的依赖只做“有没有CNV”这种二分类判断也就是只看CopyKAT的阴阳结论。7.2 冷冻样本和新鲜样本对CNV推断的影响不同样本处理方式对CNV推断的敏感度影响很大。我对比过同一患者的冷冻样本和新鲜样本新鲜样本的细胞质量明显更高inferCNV的信号也更干净冷冻样本由于细胞死亡率高、空液滴多整体表达水平偏低CNV信号被稀释。所以如果你手上只有冷冻样本的数据不要直接套用默认参数。可以考虑把cutoff调低比如从0.1降到0.05或者把denoise参数调得更激进一些让低质量细胞带来的噪声被更充分地抑制。7.3 如果不同工具之间结果冲突当inferCNV说某个细胞群是肿瘤而CopyKAT说它是二倍体时优先信任谁我的建议是以HMM结果为准以CopyKAT的阴阳分类为辅助参考。因为inferCNV的HMM模块考虑的是区域化的连续CNV信号而CopyKAT在某些情况下容易被单个染色体臂的弱信号迷惑导致漏判。如果两者真的冲突还有一种很实用的排查方法把冲突细胞单独拎出来做亚聚类看它们和周围的细胞在转录组上更接近谁。通常你会看到它们更接近肿瘤细胞群体而不是参考细胞群体——这说明CopyKAT的结果更可能是阈值敏感性造成的假阴性而不是真正常倍体。8. 我的最终分析流程与参数备忘单整套流程跑下来最终形成了一套相对稳定的分析管线。下面按顺序整理方便直接复用。这里面包含了我在多个项目中反复调整后确认比较稳妥的参数。第一步数据质控与预处理过滤标准nFeature_RNA 500nCount_RNA 1000percent.mt 10Doublet过滤用DoubletFinder或Scrublet保留score低于阈值的细胞标准化LogNormalizescale.factor 10000第二步细胞大类注释先用经典标记基因注释出所有细胞大类T/B/髓系/上皮/内皮/成纤维等确认参考细胞类型和候选肿瘤细胞群的大致位置第三步每个样本独立跑inferCNV参考细胞免疫细胞T细胞髓系细胞每类50~200个cutoff 0.1denoise TRUEHMM TRUE聚类方式cluster_by_groups TRUE这样每个注释类型在热图里是分开的方便看组间差异第四步从HMM结果中量化CNV负荷读取Pred_CNV_call.txt计算每个细胞的n_events、total_gain、total_loss计算每个细胞类型的平均CNV评分做统计检验第五步用CopyKAT做交叉验证同一个表达矩阵跑CopyKAT比较aneuploid/diploid标签和inferCNV结论的一致性不一致的细胞单独排查第六步整合可视化与论文输出UMAP上用CNV评分做FeaturePlot展示细胞群的恶性程度梯度CNV热图用于展示染色体层面的整体模式箱线图用于量化比较必要时叠加突变证据和公共数据库比对这套流程的时间成本大概是一次inferCNV跑30到60分钟取决于细胞数CopyKAT快一些10分钟左右。整个样本从原始数据到最终结论顺利的话一两天能跑完一轮。9. 关于肿瘤细胞推断一些更长期的思考肿瘤细胞身份推断这件事本质上是在做一件“用转录组的影子去推断基因组的实体”的事情。它永远不是完美的——因为转录组受到太多因素影响包括细胞周期、微环境信号、技术噪声。但它是目前单细胞水平上能够系统性区分肿瘤细胞和正常细胞的最实用路径。真正让我觉得这个流程有价值的不是某一次分析跑出了多么漂亮的热图而是它迫使你把“细胞身份”这个问题拆解到基因组层面去思考一个细胞之所以是肿瘤细胞并不是因为它表达某个基因而是因为它的基因组不再稳定携带了导致恶性行为的变异。所以如果你也在做类似的分析我的建议是不要把inferCNV当黑盒也不要只依赖一张热图下结论。把证据链建完整——CNV推断、参考细胞选择、突变锚点、公共数据库比对——每一个环节都经得起追问那么最终的结果无论是对你自己还是对审稿人都是可靠的。如果你在跑流程时遇到什么奇怪的结果欢迎回头来对比你踩的坑和我上面提到的这些是不是同一类。单细胞分析里奇怪的结果往往不是分析本身错了而是某个隐蔽的假设没有满足。