ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

单细胞肿瘤注释实战:copyKAT与inferCNV联合识别恶性细胞

单细胞肿瘤注释实战:copyKAT与inferCNV联合识别恶性细胞 做过肿瘤单细胞数据分析的朋友大概都有过这种经历降维聚类做得漂漂亮亮marker 也一个个验证过结果把注释结果拿给病理方向的合作者看对方一句“你以为的正常上皮其实才是肿瘤大本营”心态当场就崩了。我早期在实体瘤项目里就被这样怼过。后来我开始把体细胞拷贝数变异CNV当成肿瘤细胞的“身份证”来用最常用的就是 copyKAT 和 inferCNV 这一对工具。前阵子几个内部数据集跑下来两个方法互相印证准确率确实比我预想的高不少。这篇就把我的完整操作流程、关键参数、结果判读逻辑和踩过的坑一起写下来给做肿瘤微环境注释、肿瘤纯度估计和亚克隆分析的朋友做个参考。1. 为什么单细胞注释里“肿瘤细胞认不出”是个真问题1.1 传统 marker 注释为什么会在肿瘤细胞这里翻车如果你只用表达 marker 来判断恶性细胞很快会发现一个尴尬的事实肿瘤细胞压根没有一个通用的“身份证蛋白”。上皮来源的癌通常表达 EPCAM、KRT8、KRT18、KRT19但这些 marker 在正常上皮细胞里也同样高表达。你拿到一个组织样本里面可能混着正常上皮、肿瘤上皮、成纤维细胞、内皮细胞和一堆免疫细胞单看 EPCAM 阳性根本无法区分这群细胞到底是“正常的”还是“恶性的”。更麻烦的是一部分肿瘤细胞在发生上皮-间质转化后上皮 marker 的表达量会大幅下降同时开始表达 VIM 等间质基因。如果你机械地按“EPCAM 阳性 肿瘤”来注释这些已经“脱上皮化”的肿瘤细胞会被你莫名归到成纤维细胞或者未知细胞里。反过来正常组织中某些应激状态的上皮细胞也可能上调一系列增殖和转化相关基因看上去特别像恶性的。免疫细胞反而好认一些因为 CD3D、CD8A、CD68 这些谱系 marker 在肿瘤和正常组织里相对稳定。肿瘤细胞的问题在于它本身是“正常细胞突变后长歪了的版本”表达程序高度异质一个 marker 在不同亚克隆之间的差异可能比不同细胞类型之间都大。所以靠单个 marker 甚至一组 marker 来认定恶性细胞本身就是一件风险极高的事。1.2 拷贝数异常为什么能当“身份指纹”真正靠谱一点的思路是去看基因组层面的改变。实体瘤在发展过程中几乎都会积累大量体细胞拷贝数变异染色体臂的大段扩增、片段缺失、局部基因的高倍扩增等等。这些 CNV 是细胞在恶性转化过程中留下的“克隆印迹”同一个肿瘤亚克隆里的细胞基本共享同一套大片段改变。而正常体细胞在没有恶性转化的情况下绝大多数区域都保持二倍体状态。scRNA-seq 虽然不直接测 DNA但转录组数据里其实藏着拷贝数信息某个染色体区域平均拷贝数越高这个区域内所有基因检测到的 read count 通常也会成比例地偏高。这就是所谓的基因剂量效应。把一个细胞所有基因按染色体位置排开再用滑动窗口看平均表达量的高低就能大致还原出一条“染色体拷贝数轮廓”。肿瘤细胞的轮廓往往东缺一块、西多一块而正常细胞基本是一条平线。顺着这个思路社区里出现了两个主流的推断工具copyKAT 和 inferCNV。前者更像一个“肿瘤细胞分类器”直接给每个细胞打标签后者更像一个“可视化显微镜”输出全基因组的 CNV 热图。把两者结合着用是我目前觉得最稳的做法。2. copyKAT给每个细胞发一份“非整倍体体检报告”2.1 原理滑动窗口如何变成判断依据copyKAT 的核心思路不复杂三句话能说清先把基因按照染色体位置排列然后在每个染色体内部划分滑动窗口对每个窗口的 read count 信号做标准化和背景校正判断这个窗口对应的区域是正常拷贝、扩增还是缺失最后把一个细胞所有窗口的判断汇总起来计算它整体偏离二倍体的程度据此给出这个细胞是 “aneuploid” 还是 “diploid” 的预测。注意一个关键点copyKAT 直接吃原始 read count 矩阵而不是 log-normalized 后的表达矩阵。原因在于基因剂量效应是线性的raw count 的均值变化能天然反映拷贝数变化而 Seurat 里默认的 log 归一化会把这种线性关系扭曲掉导致信号失真。我见过有人贪方便直接拿GetAssayData(obj)的结果丢进 copyKAT出来的热图简直像撒了一把芝麻完全没法看。copyKAT 还有一个我很喜欢的特性——它不要求你预先提供一组已知的正常细胞作为参考。工具会从数据内部寻找二倍体基线再据此判断哪些细胞明显偏离。当然如果你确实知道某些细胞是正常的也可以显式传进去后面我会专门讲这个参数。2.2 一份可以直接改的 copyKAT 运行代码实际运行代码比较短不过参数背后的含义值得逐一说清楚。下面这段是我在 10x 数据上常用的配置library(copykat) # 矩阵要求行为基因列为细胞值是 read count # 基因名建议用 symbolid.type 对应 S res - copykat( rawmat as.matrix(counts_matrix), id.type S, ngene.chr 5, win.size 25, KS.cut 0.05, sam.name my_tumor_sample, distance.method euclidean, norm.cell.names normal_cells, # 可选传已知正常细胞的 barcode n.cores 8 )参数逐个说id.type表达矩阵里的基因名类型。S 表示 symbol如果你用的是 Ensembl ID换成 E。选错会直接导致基因坐标匹配不上。ngene.chr每个染色体最少需要多少个表达基因才参与计算。设置成 5 是为了过滤掉那些基因注释稀疏的染色体区域减少噪声。win.size滑动窗口大小窗口越小对局部 CNV 越敏感但也更容易被单个基因的表达波动带偏。25 是官方默认值在大多数实体瘤数据集上平衡性都不错。KS.cut判断窗口拷贝数状态时显著性检验的阈值默认 0.05。想更严格可以调到 0.01但有些微弱的 CNV 信号可能就检测不出来了。distance.method后续对细胞做层次聚类时用的距离度量默认 euclidean 在常见样本上表现稳定。norm.cell.names可选参数。如果通过免疫 marker 注释出了一批很明确的 T、B、NK、髓系正常细胞把它们的 barcode 传进去算法会以这些细胞作为强先验能明显降低误判。如果你的矩阵很大建议先用 Seurat 做一轮标准过滤再取子集跑 copyKAT。几万个细胞跑起来虽然不会爆内存但时间会拉得比较长。2.3 输出文件与 R 对象里到底有什么copyKAT 跑完以后会在当前工作目录生成以sam.name为前缀的几个文件my_tumor_sample_copykat_prediction.txt每个细胞的最终预测结果列里包括 barcode 和对应的 copykat.pred诊断标签是 aneuploid 或 diploid。my_tumor_sample_copykat_CNA_results.txtCNV 矩阵相当于每个细胞在每个窗口的拷贝数状态可以拿去做后续亚克隆分析或者自定义可视化。my_tumor_sample_copykat_heatmap.png默认的热图我建议自己用copykat.heatmap函数重新画因为默认图的比例和配色在论文里不一定合用。在 R 会话里最常用的是res$copykat.pred$prediction。这是一个两列的数据框你可以直接把它匹配回 Seurat 对象方便后续和聚类、marker 一起看pred - res$copykat.pred$prediction seu$copykat_pred - pred$copykat.pred[match(colnames(seu), pred$cell)] # 快速看一眼预测分布在哪些 cluster table(Idents(seu), seu$copykat_pred)读结果的时候有一点要提醒aneuploid在 copyKAT 的输出里代表“预测为非整倍体疑似肿瘤细胞”它不是一个绝对确凿的病理诊断diploid也不要直接等同于“正常细胞”因为某些低 CNV 肿瘤在单细胞分辨率下确实可能被归到 diploid。把预测当成一个强信号而不是最终结论。2.4 copyKAT 的“准”和“不那么准”我实际用下来copyKAT 对上皮来源的实体瘤类型特别有效比如头颈癌、肺癌、乳腺癌、结直肠癌这些。这些肿瘤通常携带大量大片段拷贝数改变信号足够强single cell 层面一眼就能看出谁不在二倍体状态。在我最近跑的一个样本里copyKAT 判出来的 aneuploid 细胞超过九成都落在同一个主要恶性上皮 cluster 里和病理注释高度吻合。但它也有明显的盲区。第一种情况是肿瘤整体 CNV 负荷很低比如一些分化比较好的甲状腺癌、部分血液肿瘤、低级别胶质瘤拷贝数改变很微弱copyKAT 很容易把真肿瘤细胞也判成 diploid。第二种情况是样本里免疫细胞占比极高T 细胞和 B 细胞在 T 细胞受体和 B 细胞受体区域的重组信号会被算法当成“缺失”个别克隆性扩增的 T 细胞甚至会在局部染色体上表现出类似 CNV 的信号。这个问题不只在 copyKAT 里有inferCNV 也会遇到后面专门讲。3. inferCNV用正常参照细胞把肿瘤的 CNV 轮廓照出来3.1 它在算什么相对表达 vs 相对拷贝数inferCNV 的思路最早可以追溯到 Tirosh 等人 2016 年在 Science 上发表的黑色素瘤单细胞工作。核心假设和 copyKAT 一样还是基因剂量效应但 inferCNV 更强调“相对比较”你把所有基因按染色体位置排序用滑动窗口计算每个细胞窗口内的平均表达量然后把肿瘤细胞组和一组正常参考细胞组放在一起比较看肿瘤组在哪些基因组区域系统性高于或低于参考组。正因为它是“相对”的所以 inferCNV 对参考组的选择极其敏感。你用 T 细胞做参考和用正常上皮做参考得到的信号轮廓会有实质差别。实体瘤样本里最常用的参照是肿瘤浸润的免疫细胞因为它们整体就是二倍体背景而且通常数量够多、表达程序相对稳定。需要注意是不要让参考组里混入太多有克隆扩增的免疫细胞否则受体区域会带进额外噪声。与 copyKAT 直接把每个细胞分到 two categories 不同inferCNV 输出的是一张连续信号的热图保留了大量空间信息。你看热图时能直观看到某个染色体臂是整体扩增了还是只有局部片段改变这对亚克隆分析特别有意义。3.2 三个输入文件怎么备inferCNV 要求三个输入表达矩阵、细胞注释文件、基因坐标文件。没有一个是能随便糊弄的依次来说。表达矩阵最好用文本格式行为基因、列为细胞值同样是 read count 或 UMI count。如果直接从 Seurat 导出用GetAssayData(obj, slot counts)不要用data那层。细胞注释文件是一个两列的 tab 分割文本第一列是细胞 barcode第二列是分组名。重点在于必须有一个组叫normal或你在代码里指定的参考组名inferCNV 会拿这个组的平均表达量做基线AAACCTGGTACGCAAT-1 normal AACACGTGTGTACGCT-1 tumor基因坐标文件是另一个容易踩雷的地方。常见格式是四列、无表头基因名、染色体、起始位置、终止位置。比如AKT1 chr14 104477349 104501767 TP53 chr17 7661779 7687538下载基因注释的时候尽量从 UCSC 或 Ensembl 导出一份蛋白编码基因列表过滤掉线粒体基因、核糖体蛋白基因、没有明确染色体定位的记录。还有一个容易忽略的点矩阵里的基因名如果有一部分是 symbol、一部分是 Ensembl IDinferCNV 匹配时会直接产生一堆 NA后面的结果就没法看了。我习惯跑之前用intersect统一基因范围。3.3 运行参数cutoff、denoise 和 HMM代码方面标准流程是这样library(infercnv) infercnv_obj - CreateInfercnvObject( raw_counts_matrix counts_matrix.txt, annotations_file cell_annotations.txt, gene_order_file gene_order.txt, ref_group_names c(normal) ) infercnv_obj - infercnv::run( infercnv_obj, cutoff 0.1, out_dir infercnv_output, cluster_by_groups TRUE, denoise TRUE, HMM FALSE )这里cutoff是有讲究的。官方建议 10x Genomics 数据用 0.1Smart-seq2 这种全长转录组方法用 1。它本质上是在过滤低表达基因低于阈值的基因会被丢到背景里避免太多 dropout 造成的虚假信号。如果你用 UMI 数据却把 cutoff 设成 1会过滤掉大量真实信号热图会变得很淡。denoise TRUE会在聚类后做一轮背景降噪把组内普遍存在的低水平波动抹平。HMM TRUE则是在降噪基础上进一步把连续信号转换成离散的拷贝数状态缺失、中性、扩增适合想定量比较不同样本 CNV 状态的场景但计算时间会明显增加。如果只是想先看个大概我建议第一次先设 FALSE跑通了再回头补 HMM。3.4 怎么看 inferCNV 热图inferCNV 默认输出里有一张infercnv.png横轴是基因组坐标从 1 号染色体排到 X 染色体纵轴是细胞所有细胞按注释组排列。颜色方面红色代表相对参考组表达更高、推测为扩增蓝色代表表达更低、推测为缺失。拿到热图后我的阅读顺序是固定的先看正常参考组的区域是不是一条相对平坦的条带。如果参考组自己都红一块蓝一块说明参考细胞选得有问题后面肿瘤组的信号也别急着信。再看肿瘤组区域有没有明显的“片段嵌合体”现象——不是单个基因的差异而是连续几十个窗口、整段染色体都偏红或偏蓝。这种连续区域才是真正的 CNV 信号。最后一定要看信号落在哪些染色体上。如果只在 14 号染色体、7 号染色体的受体区见信号别太兴奋后面避坑章节会具体讲为什么。4. 联合判读两个工具都说“是”才算稳4.1 我平时用的叠加流程copyKAT 和 inferCNV 各有优势copyKAT 快能给每个细胞独立的二倍体/非整倍体标签inferCNV 慢但能给出全基因组轮廓能看见大片段改变的边界。合在一起用的时候我有一套固定流程先用 Seurat 或 Scanpy 做标准 QC、聚类和细胞大类注释得到 T 细胞、B 细胞、髓系细胞、上皮细胞、成纤维细胞等大群。跑 copyKAT得到每个细胞的预测标签快速定位疑似非整倍体的细胞群。从注释结果里挑出靠谱的正常细胞通常是 T 细胞、NK 细胞、巨噬细胞等免疫细胞作为 inferCNV 的参考组。跑 inferCNV获得肿瘤细胞相对参考组的 CNV 热图。把 copyKAT 预测结果、inferCNV 分组信号和细胞类型 marker 放到同一个表格里按 cluster 逐群比对。只在两者结论一致、且 marker 支持的情况下才把细胞定义为高置信肿瘤细胞。用这组高置信肿瘤细胞做后续差异表达、拟时间分析或肿瘤纯度估计。先跑 copyKAT 再跑 inferCNV 的原因很实际copyKAT 速度快、判读简单先给你一个大方向inferCNV 算得慢但能验证方向同时补充 CNV 的具体位置和边界。如果一开始直接 inferCNV你可能面对一张巨大的热图非常迷茫不知道哪些群是重点。4.2 结果不一致时怎么拉架两个工具的结果不可能永远一致遇到不一致也没什么好慌的按下面这个思路处理就够了。我做了个简单的判断表基本覆盖了常见情况copyKAT 标签inferCNV 信号更可能的情况建议动作aneuploid有清晰 CNV 块高置信肿瘤候选结合 marker 和染色体位置做最终确认aneuploid无明显信号局部/微小 CNV或近整倍体肿瘤不轻易排除看 marker 并复核亚克隆diploid有明显 CNV 块参考组不干净或某群表达程序被当成了 CNV检查参考组细胞构成看信号位置diploid无明显信号正常细胞或低 CNV 肿瘤以其他证据辅助判断第二种情况在实体瘤里不算少见尤其是一些携带点突变驱动的肿瘤大片段 CNV 本来就不多。遇到时我会退回去看差异表达结果如果这个群同时高表达一系列肿瘤相关通路基因即使拷贝数层面没有强信号也不能草率地划成正常细胞。第三种情况更值得警惕。如果 inferCNV 的参考组里混进了非整倍体细胞比如你注释时把一群上皮来源的肿瘤细胞误标成了“正常上皮”那整个基线都会偏移肿瘤细胞的信号反而可能被抹平。此时 copyKAT 的 diploid 判断反而更可信。4.3 marker 是第三把尺子copyKAT 和 inferCNV 都是统计推断不是金标准。最终过 marker 这一关不能省。这里说的不是拿单个 EPCAM 就拍板而是用一组表征谱系和恶性状态的 marker 去交叉验证。举个例子如果一个细胞群被 copyKAT 判成 aneuploidinferCNV 热图上也有连续 CNV 信号但它的 CD3D、CD8A 表达和 T 细胞群完全重叠那这大概率不是实体瘤细胞而是发生克隆性扩增的 T 细胞。反过来一个群呈上皮 marker 阳性同时又有明显 CNV那基本就是恶性上皮细胞没跑了。我习惯用的 marker 组合至少包括 EPCAM/KRT8/KRT18、PTPRC、CD3D、MS4A1、CD68、COL1A1、PECAM1覆盖主要细胞大类再根据癌种追加组织特异 marker。这一步看着繁琐但能避免后面差异表达分析里出现系统性假象。你想象一下如果肿瘤细胞群里混了一大堆克隆扩增的 T 细胞后面所有“肿瘤特异性”的差异基因都会被免疫程序污染审稿人一复核就露馅。5. 最容易翻车的五个细节这五个坑我都踩过5.1 直接用 log-normalized 矩阵跑结果全乱这是新手上路最常见的坑。Seurat 的默认数据槽是data也就是 log-normalized 后的连续浮点值不是原始计数。copyKAT 和 inferCNV 都要求 read count 或 UMI count 这种整数计数。如果你把 log 值丢进去模型的离散计数假设直接失效窗口内的平均信号会被高表达基因严重主导出来全是假信号。正确做法是显式取 counts 槽counts_mat - GetAssayData(seu, slot counts) counts_mat - as.matrix(counts_mat)跑之前顺手做一次质量检查class(counts_mat)应该是matrixcounts_mat[1:5, 1:5]应该是一堆整数。如果看到小数点赶紧重来。5.2 基因坐标文件的版本和格式坑基因坐标这事看起来简单炸起来也是真的炸。最常见的是染色体编号不统一有的注释文件写chr1有的写1inferCNV 匹配时会把它们当成完全不同的字符串导致大量基因被丢掉。解决方式很简单预处理时统一加chr前缀。另一个坑是参考基因组版本。copyKAT 里有genome参数hg19 和 hg20 分别对应 GRCh37 和 GRCh38你提供的基因坐标必须和这个参数一致。很多教程里给的是 hg19 的坐标文件而你的比对又是 GRCh38 的结果坐标错位CNV 边界全弯掉。我的习惯是构建一套固定的 gene order 文件放在项目里反复用每次新样本都先验证版本。重复基因名也要处理。Ensembl 导出时同一个基因可能出现多行不同转录本如果你不先做aggregate或distinct下游会报“duplicate row names”。线粒体基因和核糖体基因同样是高频噪声源通常直接过滤掉因为它们的高表达与拷贝数状态关系不大。5.3 内存、时间与“跑一半崩了”inferCNV 是真吃内存这不是开玩笑。三万个细胞的全量运行加上 denoise 和聚类64GB 内存的机器很容易直接 OOM。我第一次跑一个两万多细胞的数据集整个工作站卡到鼠标都动不了最后只能分群跑。我的建议有三个第一输入前先过滤把低质量细胞和明显不在分析范围内的细胞去掉保留代表性细胞做 inferCNV比如每群随机抽一部分第二设一个合理的out_dirinferCNV 会把中间结果写盘崩了能续第三copyKAT 的n.cores不要贪多8 核和 32 核在超大矩阵上的差距不一定大但内存开销会成倍增加。如果你确实需要全量分析几万个细胞建议把数据切成几个亚群分组跑最后再按染色体位置合并信号。虽然麻烦但至少不会跑到一半被系统杀死。5.4 免疫受体区域和 X 染色体的假信号这个坑值得单独拿出来说因为几乎所有免疫细胞富集的样本都会遇到。T 细胞受体基因位于 7 号和 14 号染色体B 细胞受体基因主要位于 2 号、14 号和 22 号染色体这些区域在 V(D)J 重组过程中会被切割和重排在 scRNA-seq 里就会表现为持续的“缺失”或“扩增”信号和真 CNV 长得一模一样。MHC 区域也是一样6 号染色体 p21 区域的 HLA 基因多态性和表达异质性极高如果拿髓系细胞做参考这一区域经常会显示出片状信号。我的处理方法是先观察信号是否集中在上述区域如果 copyKAT 和 inferCNV 的信号只在受体区出现那就不能认定为肿瘤特异性 CNV更常见的是免疫细胞正常的受体编辑痕迹。还有一个容易忽略的是 X 染色体。男女样本混在一起跑的时候女性细胞的 X 染色体通常会因为剂量补偿机制表现出整体表达偏高看上去就像整条 X 染色体扩增。遇到 X 染色体上的信号先确认样本性别和参考细胞性别再决定要不要把它当作真实 CNV。5.5 没有正常参考组的时候怎么办inferCNV 对参考组是刚需没有 normal 组就没法算相对信号。最理想的情况是在实验设计阶段就预留正常组织或确定无疑的正常细胞。但很多回顾性数据分析里并没有这个条件尤其是你手里只有一份纯肿瘤样本时。这时候有个临时办法用 copyKAT 预测为 diploid 的那群细胞作为参考组。注意这只是一个“近似正常”参考因为前面说过低 CNV 的肿瘤细胞也可能被 copyKAT 判成 diploid。更稳的做法是从所有二倍体细胞里再筛一遍 marker只保留明确的 T 细胞、NK 细胞、巨噬细胞作为参考。不少人也试过用公共数据集里的正常组织单细胞做参考理论上可行但批次效应会带来额外噪声同一个基因在两个数据集的测序深度差异可能被误判成拷贝数差异。从我实测结果看还是样本内部挑选参考细胞最靠谱实在不行再考虑外部数据而且一定要做好归一化和批次校正。至于 copyKAT它虽然不强制要求参考但如果你有已知的正常细胞建议照样传进norm.cell.names。在一个肿瘤细胞占比特别高的样本里内部二倍体基线会被肿瘤信号拉偏显式指定正常细胞能显著减少假阴性。最后说一个我私人的工作习惯跑完 copyKAT 和 inferCNV 后我从来不会只看两个工具的标签就直接写结论而是把它们和 marker 表达、聚类结果放进同一张表按 cluster 逐群过一遍。是肿瘤的细胞群通常不会只在其中一个工具里有强信号也不会突然高表达一堆免疫谱系基因。看似多花了十几分钟做交叉验证但这一步能帮你省下后面无数次的反复质疑。做肿瘤细胞注释这件事慢一点、稳一点比什么都强。
返回列表