本教程基于Linux中的Rstudio环境演示，计算资源不足的同学可参考：生信分析为什么要使用服务器？足够支持你完成硕博生涯的生信环境 配置一个心仪的工作站(硬件+环境配置) 独享服务器，生信分析不求人 为实验室准备一份生物信息学不动产

二、可视化集锦 三、绪论 1、定义 拷贝数变异（Copy number alterations，以下简称CNA）是一种重要的基因组变异，在癌症的发生和发展过程中起着至关重要的作用。确定肿瘤细胞中CNA的特征对早期肿瘤检测、划分肿瘤异质性、了解肿瘤进展模式和揭示耐药机制具有重要意义。随着单细胞RNA测序（scRNA-seq）的蓬勃发展，研究人员开发出了多种计算方法来从scRNA-seq研究中推断CNA。

目前从scRNA-seq数据中推断CNA的计算方法大致可分为两类。一类是基于表达的，即直接从基因表达推断CNA。这类工具包括inferCNV、copykat和SCEVAN。这些工具的工作原理是，CN扩增（AMPs）或缺失（DELs）将分别导致受影响基因组区域内基因的上调或下调。InferCNV利用基因窗口的校正移动平均数；copykat采用综合贝叶斯分割方法；SCEVAN采用多通道分割算法。

第二类将基因表达与等位基因信息相结合，相关工具有Numbat和CaSpER等。2、主要研究成果 这里我们参考了一篇2025年发表Briefings in Bioinformatics杂志上面的文献，文献的标题为“Benchmarking copy number aberrations inference tools usingle-cell multi-omics datasets”。

作者对五种拷贝数变异（CNA）推断工具进行了基准测试。这些工具的选择基于两个主要标准：（i）需单细胞RNA测序（scRNA-Seq）数据作为输入；（ii）在该领域内被公认为流行工具。这些工具大致可分为两类：仅需表达矩阵的（inferCNV、CopyKAT和SCEVAN），以及需要表达矩阵和等位基因频率数据的（Numbat和CaSpER）。其性能通过肿瘤细胞中的三个常见应用进行评估：（i）肿瘤细胞与正常细胞的分类；

（ii）拷贝数变异事件的准确性；（iii）肿瘤亚克隆推断。2.1 区分肿瘤细胞与正常细胞 为了根据推断出的拷贝数变异（CNA）谱区分肿瘤细胞和正常细胞，作者将上述工具应用于包括8例结直肠癌（CRC）、1例神经内分泌肿瘤和1例胶质瘤在内的实体瘤。所有细胞的细胞类型注释均来自相应的已发表研究。仅保留结直肠癌中的上皮细胞、神经内分泌肿瘤中的内分泌细胞以及胶质瘤中的神经胶质细胞用于肿瘤细胞与正常细胞的分类。

总体而言，Numbat在区分肿瘤细胞与正常细胞方面表现最佳。在仅使用表达矩阵的三个工具中，CopyKAT的整体性能最佳。SCEVAN在肿瘤细胞频率较低（<40%）的样本中表现不佳（图2B）。InferCNV除在CRC02中表现不可靠（其中大部分细胞为肿瘤细胞，约88%）外，其余情况均可靠。CaSpER也犯了同样的错误。

为解决这一问题，作者研究了在肿瘤与正常细胞比例失衡的情况下，纳入肿瘤微环境（TME）细胞（如免疫细胞、内皮细胞和成纤维细胞）是否能提高inferCNV、CopyKAT和SCEVAN的性能（图2E）。结果表明，纳入TME细胞显著提高了SCEVAN对肿瘤细胞预测的准确性（图2E）。这种改进是合理的，因为SCEVAN算法整合了来自公共数据库的一系列基因特征，包括来自TME的细胞（如基质细胞、免疫细胞），以识别高可信度的正常细胞。

因此，添加TME细胞增强了SCEVAN区分肿瘤细胞和正常细胞的能力。尽管在inferCNV中未观察到整体的改进，但某些肿瘤细胞数量较多的样本，如CRC02（88%）、CRC22（91%）和CRC03（67%），其F1分数却有了显著提升，分别从0、0.5和0.55上升至1、0.9和0.99（图2E）。

为了系统地探究肿瘤纯度对拷贝数变异（CNA）推断准确性的影响，作者生成了一个涵盖广泛肿瘤/正常细胞比例（从1:100到100:1）的合成数据集，以使用最大的胶质瘤数据集来考察其对性能的影响（图2H）。即使在肿瘤/正常细胞比例极低的极端情况下，Numbat也始终优于其他工具（图2I）。

为了探究测序深度对拷贝数变异（CNA）推断准确性的影响，作者对几个样本进行了下采样，使其测序深度的中位数分别降至约10,000、3,000和1,000个UMI/细胞。随着测序深度的降低，所有工具的整体肿瘤-正常细胞分类F1分数均有所下降，尤其是在Numbat中（图2J）。然而，CopyKAT的性能基本不受影响。

备注：F1 score是一种用于评估分类模型性能的指标，它综合了精确率（Precision）和召回率（Recall）这两个指标，能够更全面地反映模型的分类效果。其中，精确率是指预测为正例的样本中实际为正例的比例，召回率是指实际为正例的样本中被正确预测为正例的比例。F1 score的取值范围在0到1之间。F1 score越高，说明模型的性能越好。

当F1 score为1时，表示模型的精确率和召回率都达到了100%，即所有的预测结果都是正确的，这是一种理想的完美状态；当F1 score为0时，则表示模型的预测结果完全错误，没有一个正例被正确预测出来。2.2 推断的拷贝数变异（CNA）谱的准确性 为了评估推断出的拷贝数变异（CNA）谱的准确性，确保每个软件工具的最佳性能至关重要。

从单细胞RNA测序（scRNA-seq）数据推断CNA的过程中，正常参考细胞的指定对最终的CNA结果有显著影响。Numbat和CaSpER软件需要输入参考细胞，最好来自相同的生物学和技术条件。对于InferCNV、CopyKAT和SCEVAN软件，使用参考细胞是可选的。

在9个案例中的8个案例中，指定了参考细胞之后，从单细胞DNA测序（scDNA-Seq）得出的拷贝数变异（CNA）图谱与从单细胞RNA测序（scRNA-seq）推断出的CNA图谱之间的相似性显著提高（图3B-C）。相比之下，CopyKAT和SCEVAN在其算法中已经采用了类似两步走的方法，即先预测正常参考细胞，然后将其作为基线来校正肿瘤细胞的CNA。

Numbat和CaSpER生成了整数倍的拷贝数（CN）图谱，这些图谱更易于解读且看起来更清晰，而inferCNV、CopyKAT和SCEVAN输出的是连续的CN。在所有情况下，没有一种工具在单细胞RNA测序（scRNA-Seq）推断图谱与单细胞DNA测序（scDNA-Seq）图谱的一致性方面始终优于其他工具（图3D-E）。在肿瘤细胞和正常细胞数量都足够的情况下，例如胶质瘤、CRC13、CRC11和CRC12，所有工具的表现都很好。

在肿瘤细胞和正常细胞数量不平衡的样本中，例如CRC02、CRC22和CRC15，这五种工具预测的拷贝数变异（CNAs）图谱的整体准确性大幅下降。作者评估了在血液系统癌症ALL中的表现。但由于原始FASTQ数据不可用，作者只能评估inferCNV、CopyKAT和SCEVAN的性能。其中，CopyKAT表现最佳。此外，作者还评估了这些工具检测断点的能力。

除inferCNV外，其他四个工具均能在单细胞水平上输出拷贝数变异（CNV）片段（图3F）。这些工具在每个细胞中检测到的片段数量存在较大差异。Numbat和CaSpER检测到的片段数量为数十到数百个，而CopyKAT和SCEVAN则检测到数千个片段，这显然超出了癌症中预期的染色体断点数量。在单细胞水平上评估断点检测的技术难度仍然很大。因此，作者重点评估了inferCNV、Numbat和SCEVAN在肿瘤克隆断点检测方面的准确性。

结果显示，SCEVAN在F1分数和灵敏度方面总体表现最佳，而inferCNV在精度方面表现更优。与SCEVAN和来自DNA-Seq的真实情况相比，inferCNV和Numbat检测到的断点要少得多。更具体地说，inferCNV和Numbat较低的F1分数断点可能是由于在识别复杂拷贝数变异（CNAs）时分辨率降低所致。

例如，在chr10上，SCEVAN检测到两个断点，而Numbat和inferCNV将这两个相邻的缺失合并为一个缺失（图3H）。这也可能是由于断点检测遗漏造成的，例如在chr14上（图3I）。最后，在Numbat中偶尔也会观察到假阳性CNAs调用（图3J）。由于整合了AB等位基因分析，Numbat和CaSpER均能够检测到杂合性缺失（LOH）。

由于CaSpER在恶性细胞和非恶性细胞中均做出了过多明显错误的LOH调用，因此未对其LOH调用性能进行进一步评估。Numbat在多个样本中总共检测到17个LOH事件。为了评估其准确性，作者使用配对的单细胞DNA测序（scDNA-seq）数据对这些LOH事件的实际拷贝数变异（CNV）状态进行了分析。分析表明，Numbat检测到的大多数LOH事件（17个中的14个，约82%）实际上是CNV扩增（AMP）或缺失（DEL）（图3K）。

例如，在CRC13样本中，Numbat将chr8q识别为LOH事件；然而，对DNA数据的CNA图谱分析显示，chr8q实际上经历了CNV扩增（图3L）。Numbat的这一错误LOH调用是由于8q区域A等位基因显著增加，而B等位基因保持不变，正如DNAseq分析所观察到的（图3M）。但Numbat在识别LOH方面表现出高灵敏度。

2.3 亚克隆结构推断准确性的评估 人类癌症表现出广泛的瘤内异质性，通过不同亚克隆中不断出现的突变和拷贝数变异（CNAs）持续进化，这会影响其表型并赋予适应性优势。单细胞拷贝数（CN）谱可以推断出亚克隆以及肿瘤的进化史。在胶质瘤中，先前的研究已经确定了亚克隆结构，只有Numbat正确地将较小的C1亚克隆且具有较低的拷贝数变异负荷的细胞归类为肿瘤细胞（图4A-C）。

对于结直肠癌（CRC）病例，作者通过平均轮廓系数来确定最优的k值，估计出的结直肠癌样本的最优克隆数为2。如果通过不同方法预测出的亚克隆数大于2，则根据其拷贝数变异（CNV）谱的相似性合并这些簇。经人工检查，两个样本（CRC03和CRC11）显示出两个较大的亚克隆。在CRC03的情况下，只有inferCNV出现了类似的肿瘤识别错误（图4D）。inferCNV错误地将肿瘤C1亚克隆识别为正常细胞，从而未能正确检测出两个亚克隆。

但是，将肿瘤微环境（TME）细胞添加到inferCNV中纠正了错误分类，并改善了亚克隆结构的推断。其他四个工具在亚克隆分配方面表现良好（图4E-F），尽管单细胞DNA测序（scDNA-Seq）获得的亚克隆CNAs谱并未完全重现。在CRC11中，所有工具都正确地对肿瘤细胞进行了分类，所有工具都取得了良好的性能（所有工具的ARI值均大于0.8，图4G-I）。

另外，在一个极端案例中，CRC12中定义的一个亚克隆仅由一个单细胞组成，而inferCNV和SCEVAN正确检测到了这种单细胞亚克隆。因此，在肿瘤细胞分类正确的前提下，这五种工具在描绘亚克隆结构方面均表现出色。2.4 计算效率 为了评估各种方法的计算速度，作者追踪了所用的CPU时间。结果表明，在运行时间方面，CopyKAT和SCEVAN表现最佳。InferCNV和CaSpER的计算速度中等，而Numbat则需要最多的CPU时间。

此外，随着数据集规模的增大，运行时间略有增加。不过，对于所测试的所有工具，运行时间都在可接受范围内。从BAM文件计算B等位基因频率通常是耗时最长且资源消耗最大的步骤。根据作者的测试，配备普通配置的笔记本电脑能够轻松处理所有五种工具的少于1000个细胞的数据集。然而，对于包含数千个细胞的数据集，建议在服务器或高性能计算平台上进行分析，尤其是对于Numbat和InferCNV。

3、小结 总体而言，以上工具都能在一定程度上区分肿瘤细胞和正常细胞，并准确推断出拷贝数变异（CNA）的特征。在这些方法中，Numbat在各种评估标准中表现最佳。对于仅提供表达矩阵的情况，建议优先使用CopyKAT方法。基于基因表达的方法更容易受到肿瘤纯度的影响。在没有参考设置的情况下，inferCNV依靠输入的细胞群体来生成背景。当肿瘤纯度较高时，实际的拷贝数变异（CNA）信号可能会被错误地视为背景并被移除。

相反，当肿瘤纯度较低时，SCEVAN倾向于错误地将非恶性细胞归类为恶性细胞。对于这两种情况，将肿瘤微环境（TME）细胞纳入分析都能明显提高性能。整合B等位基因信息使Numbat和CaSpER能够识别出LOH事件，这些事件在癌症进展中起着关键作用。Numbat在检测LOH方面具有较高的灵敏度，但特异性较低。因此，用户在解读Numbat的杂合性丧失结果时应谨慎。建议使用独立工具确认这些杂合性丧失区域不存在拷贝数改变。

另一个观察结果是，拷贝数变异（CNAs）的复杂性和负担对每种方法的影响是相同的。不同工具在实体瘤上的表现通常优于在血液系统肿瘤或细胞系上的表现。这可能是由于血液系统肿瘤中的拷贝数变异较少且染色体拷贝数变异的异质性较低。鉴于没有一种工具能在所有任务中都表现出色，特殊情况下我们也可以结合两种及以上工具进行联合分析。

四、infercnv 1、简介 inferCNV是 TrinityCTAT 的一个组件，被用于在单细胞数据中鉴定体细胞大规模的染色体拷贝数变异，例如染色体部分区域的拷贝数增加与缺失。在软件执行过程中，infercnv会与参考的正常细胞之间进行肿瘤基因组中基因密度的评估，最终输出热图展示各染色体的基因相对密度，从而反映待测细胞类型的拷贝数变异情况。

InferCNV 提供了多个残差表达过滤器，以探索最小化噪音并进一步揭示支持拷贝数变异(copy number alteration, CNA)的信号。此外，InferCNV还包括预测CNA区域的方法，并根据异质性模式定义细胞聚类。

2、安装 if ( ! requireNamespace ( "BiocManager" , quietly = TRUE )) install.packages ( "BiocManager" ) if ( ! requireNamespace ( "infercnv" , quietly = TRUE )) BiocManager :: install ( "infercnv" ) 3、输入文件 一共需要准备三个输入数据：raw count矩阵、细胞注释、基因注释，这里我们以 infercnv 内置的文件展示数据结构：# 清空环境变量 rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 836583 44.7 1578222 84.3 1277907 68.3 # Vcells 1551361 11.9 8388608 64.0 2518594 19.3 # 展示系统时间 Sys.time # [1] "2025-08-29 17:53:28 CST" # 设置工作目录 setwd ( "/data/09_scRNA-seq_cnv/" ) # 如果你成功安装了infercnv，那么你的系统中自带以下文件，无需下载：# 获取测matrix的路径：matrix_path <- system.file ( "extdata" , "oligodendroglioma_expression_downsampled.counts.matrix.gz" , package = "infercnv" ) matrix_path # [1] "/usr/local/lib/R/site-library/infercnv/extdata/oligodendroglioma_expression_downsampled.counts.matrix.gz" # 比较诡异的是作者声称这里输入的是raw count(如果你是用Seurat对象生成的矩阵，GetAssayData(Seurat_obj_V5, assay="RNA", slot='counts'))，但实际上作者提供的是normalization后的矩阵(包含小数，如果你是用Seurat对象生成的矩阵，那么可以这么获取 GetAssayData(Seurat_obj_V5, assay="RNA", slot='data')) # 原文对counts的描述：InferCNV is compatible with both smart-seq2 and 10x single cell transcriptome data, and presumably other methods (not tested). The counts matrix can be generated using any conventional single cell transcriptome quantification pipeline, yielding a matrix of genes (rows) vs. cells (columns) containing assigned read counts. my_matrix <- read.table (matrix_path, header = T, row.names = 1 , sep = ' \t ' ) # 读取这个压缩文件 my_matrix[ 1 : 4 , 1 : 4 ] # 可以看到行名是基因，列名是细胞名 # MGH54_P16_F12 MGH54_P12_C10 MGH54_P11_C11 MGH54_P15_D06 # A2M 0 0.000 0.000 0.000 # A4GALT 0 0.000 0.000 0.000 # AAAS 0 37.008 30.935 21.011 # AACS 0 0.000 0.000 0.000 # 获取细胞注释的路径：anno_path <- system.file ( "extdata" , "oligodendroglioma_annotations_downsampled.txt" , package = "infercnv" ) anno_path # [1] "/usr/local/lib/R/site-library/infercnv/extdata/oligodendroglioma_annotations_downsampled.txt" # 其中包含的是细胞的注释信息，即Seurat中的meta.data anno_data <- read.table (anno_path, header = F, sep = ' \t ' ) # 读取txt文件 head (anno_data) # 查看表头，第一列包含的是细胞名/barcode名称，第二列包含的是细胞对应的注释 # V1 V2 # 1 MGH54_P2_C12 Microglia/Macrophage # 2 MGH36_P6_F03 Microglia/Macrophage # 3 MGH53_P4_H08 Microglia/Macrophage # 4 MGH53_P2_E09 Microglia/Macrophage # 5 MGH36_P5_E12 Microglia/Macrophage # 6 MGH54_P2_H07 Microglia/Macrophage print ( c ( '数据集中包含的细胞类型有：' , paste ( unique (anno_data[, 2 ])))) # [1] "数据集中包含的细胞类型有：" "Microglia/Macrophage" # [3] "Oligodendrocytes (non-malignant)" "malignant_MGH36" # [5] "malignant_MGH53" "malignant_97" # [7] "malignant_93" # 获取基因注释的路径：gene_ann_path <- system.file ( "extdata" , "gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt" , package = "infercnv" ) gene_ann <- read.table (gene_ann_path) # 包含的是基因在染色体中的坐标信息，通过这个信息才可以进行对应基因组区域的拷贝数计算 # 四列分别为基因名、染色体名、基因在染色体上的起始位点、基因在染色体上的终止位点 head (gene_ann) # V1 V2 V3 V4 # 1 WASH7P chr1 14363 29806 # 2 LINC00115 chr1 761586 762902 # 3 NOC2L chr1 879584 894689 # 4 MIR200A chr1 1103243 1103332 # 5 SDF4 chr1 1152288 1167411 # 6 UBE2J2 chr1 1189289 1209265 4、infercnv workflow 准备好以上测试文件后，两步法即可完成infercnv的分析流程：library (infercnv) library (Seurat) # 创建递归目录 dir.create ( "/data/09_scRNA-seq_cnv/result_infercnv/test1" , recursive = T) # （1）、构建infercnv对象 # 依次填入上述文件路径即可：infercnv_obj <- CreateInfercnvObject ( raw_counts_matrix= matrix_path, # 矩阵文件路径 annotations_file= anno_path, # 细胞注释文件路径 delim= " \t " , # 读取文件时的分隔符，如果是csv文件，则用',' gene_order_file= gene_ann_path, # 基因注释文件路径 ref_group_names= c ( "Microglia/Macrophage" , "Oligodendrocytes (non-malignant)" )) # 正常参考细胞类型的名称，需要根据大家的生物学知识来判断，设定一定未发生肿瘤变异的细胞类型 # 如果出现以下报错：# Error in `ScaleData`: # ! No layer matching pattern 'data' found. Please run NormalizeData and retry # Backtrace: # 1. infercnv::run(...) # 可以设置：options ( "Seurat.object.assay.version" = "v3" ) # options 函数是 R 语言中用于管理全局选项的函数 # 这行代码的主要功能是将 Seurat 对象的分析版本设置为版本V3，以便在后续的分析中使用该版本的功能和数据处理方式 # 参考https://github.com/broadinstitute/infercnv/issues/664 # 运行infercnv_obj <- infercnv :: run (infercnv_obj, # 创建的infercnv对象 cutoff= 0.1 , # 使用Smart-seq2这类全长数据时填cutoff=1 , 使用10x Genomics这类3'测序策略数据时填cutoff=0.1 out_dir= '/data/09_scRNA-seq_cnv/result_infercnv/test1' , # 输出路径 cluster_by_groups= TRUE , # 决定是否按提供的细胞分组进行聚类。

如果 TRUE，细胞会基于之前分配的组信息进行聚类，这有助于分析不同组间的 CNV 变化。如果 FALSE，则不使用组信息进行聚类。denoise= TRUE , # 决定是否对数据进行去噪处理 HMM= TRUE , # 是否启用隐马尔科夫模型（Hidden Markov Model, HMM）来进一步细化拷贝数变异的推断。HMM 会根据推断出的 CNV 模式，进行更精确的基因组段拷贝数变异的检测，帮助识别潜在的基因组异常区域。

num_threads= 4 , # 并行计算设置，以4线程为例 write_expr_matrix = TRUE ) 经过22个内置step的运行后程序运行完毕，得到如下结构的输出文件夹：system ( 'tree /data/09_scRNA-seq_cnv/result_infercnv/test1' , intern = T) # [1] "\033[01;34m/data/09_scRNA-seq_cnv/result_infercnv/test1\033[00m" # [2] "├── 01_incoming_data.infercnv_obj" # [3] "├── 02_reduced_by_cutoff.infercnv_obj" # [4] "├── 03_normalized_by_depthHMMi6.infercnv_obj" # [5] "├── 04_logtransformedHMMi6.infercnv_obj" # [6] "├── 08_remove_ref_avg_from_obs_logFCHMMi6.infercnv_obj" # [7] "├── 09_apply_max_centered_expr_thresholdHMMi6.infercnv_obj" # [8] "├── 10_smoothed_by_chrHMMi6.infercnv_obj" # [9] "├── 11_recentered_cells_by_chrHMMi6.infercnv_obj" # [10] "├── 12_remove_ref_avg_from_obs_adjustHMMi6.infercnv_obj" # [11] "├── 14_invert_log_transformHMMi6.infercnv_obj" # [12] "├── 15_tumor_subclustersHMMi6.leiden.infercnv_obj" # [13] "├── 17_HMM_predHMMi6.leiden.hmm_mode-subclusters.cell_groupings" # [14] "├── 17_HMM_predHMMi6.leiden.hmm_mode-subclusters.genes_used.dat" # [15] "├── 17_HMM_predHMMi6.leiden.hmm_mode-subclusters.infercnv_obj" # [16] "├── 17_HMM_predHMMi6.leiden.hmm_mode-subclusters.pred_cnv_genes.dat" # [17] "├── 17_HMM_predHMMi6.leiden.hmm_mode-subclusters.pred_cnv_regions.dat" # [18] "├── 18_HMM_pred.Bayes_NetHMMi6.leiden.hmm_mode-subclusters.mcmc_obj" # [19] "├── 19_HMM_pred.Bayes_NetHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.infercnv_obj" # [20] "├── 20_HMM_pred.repr_intensitiesHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.infercnv_obj" # [21] "├── 22_denoiseHMMi6.leiden.NF_NA.SD_1.5.NL_FALSE.infercnv_obj" # [22] "├── \033[01;34mBayesNetOutput.HMMi6.leiden.hmm_mode-subclusters\033[00m" # [23] "│ ├── cellProbs.0.5.pdf" # [24] "│ ├── cnvProbs.0.5.pdf" # [25] "│ ├── CNV_State_Probabilities.dat" # [26] "│ ├── infercnv.NormalProbabilities.PostFiltering.heatmap_thresholds.txt" # [27] "│ ├── infercnv.NormalProbabilities.PostFiltering.observation_groupings.txt" # [28] "│ ├── infercnv.NormalProbabilities.PostFiltering.observations_dendrogram.txt" # [29] "│ ├── infercnv.NormalProbabilities.PostFiltering.png" # [30] "│ ├── infercnv.NormalProbabilities.PreFiltering.heatmap_thresholds.txt" # [31] "│ ├── infercnv.NormalProbabilities.PreFiltering.observation_groupings.txt" # [32] "│ ├── infercnv.NormalProbabilities.PreFiltering.observations_dendrogram.txt" # [33] "│ ├── infercnv.NormalProbabilities.PreFiltering.png" # [34] "│ └── MCMC_inferCNV_obj.rds" # [35] "├── expr.infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.dat" # [36] "├── expr.infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.dat" # [37] "├── expr.infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.dat" # [38] "├── expr.infercnv.dat" # [39] "├── expr.infercnv.preliminary.dat" # [40] "├── HMM_CNV_predictions.HMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_genes.dat" # [41] "├── HMM_CNV_predictions.HMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_regions.dat" # [42] "├── infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.heatmap_thresholds.txt" # [43] "├── infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.observation_groupings.txt" # [44] "├── infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.observations_dendrogram.txt" # [45] "├── infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.observations.txt" # [46] "├── infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.png" # [47] "├── infercnv.17_HMM_predHMMi6.leiden.hmm_mode-subclusters.references.txt" # [48] "├── infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.heatmap_thresholds.txt" # [49] "├── infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.observation_groupings.txt" # [50] "├── infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.observations_dendrogram.txt" # [51] "├── infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.observations.txt" # [52] "├── infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.png" # [53] "├── infercnv.19_HMM_pred.Bayes_Net.Pnorm_0.5.references.txt" # [54] "├── infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.heatmap_thresholds.txt" # [55] "├── infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.observation_groupings.txt" # [56] "├── infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.observations_dendrogram.txt" # [57] "├── infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.observations.txt" # [58] "├── infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.png" # [59] "├── infercnv.20_HMM_predHMMi6.leiden.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.references.txt" # [60] "├── infercnv.heatmap_thresholds.txt" # [61] "├── infercnv.observation_groupings.txt" # [62] "├── infercnv.observations_dendrogram.txt" # [63] "├── infercnv.observations.txt" # [64] "├── infercnv.png" # [65] "├── infercnv.preliminary.heatmap_thresholds.txt" # [66] "├── infercnv.preliminary.observation_groupings.txt" # [67] "├── infercnv.preliminary.observations_dendrogram.txt" # [68] "├── infercnv.preliminary.observations.txt" # [69] "├── infercnv.preliminary.png" # [70] "├── infercnv.preliminary.references.txt" # [71] "├── infercnv.references.txt" # [72] "├── infercnv_subclusters.heatmap_thresholds.txt" # [73] "├── infercnv_subclusters.observation_groupings.txt" # [74] "├── infercnv_subclusters.observations_dendrogram.txt" # [75] "├── infercnv_subclusters.png" # [76] "├── preliminary.infercnv_obj" # [77] "└── run.final.infercnv_obj" # [78] "" # [79] "1 directory, 75 files" 文件比较多，它们收录了分析过程中保存的信息：01_incoming_data.infercnv_obj：原始创建的infercnv对象，包含了三个输入文件的信息。

02_reduced_by_cutoff.infercnv_obj : 特定阈值过滤后的数据。03_normalized_by_depthHMMi6.infercnv_obj : 基于HMM深度归一化后的表达矩阵。04_logtransformedHMMi6.infercnv_obj : 对基于HMM深度归一化后的表达矩阵进行对数转化的结果。

08_remove_ref_avg_from_obs_logFCHMMi6.infercnv_obj : 从观测数据中去除参考细胞平均值后的表达数据，可以突出观测细胞的异常表达情况。

09_apply_max_centered_expr_thresholdHMMi6.infercnv_obj 中心化表达阈值后的数据，减少离群值的干扰 10_smoothed_by_chrHMMi6.infercnv_obj : 按染色体坐标平滑化后的数据，减少局部波动对整体结果的影响。

11_recentered_cells_by_chrHMMi6.infercnv_obj : 以染色体为单位重新中心化的细胞表达数据，减少染色体的影响，提升结果稳定性。12_remove_ref_avg_from_obs_adjustHMMi6.infercnv_obj : 移除参考平均值并做调整后的数据，可以提升细胞异常表达的观测效果。

14_invert_log_transformHMMi6.infercnv_obj : 对数变换的反向操作，将之前的对数转换数据还原到原来的表达值范围。15_no_subclusteringHMMi6.infercnv_obj : 未进行亚群聚类的对象。17-22_HMM_pred* 系列文件: 基于HMM（隐藏马尔科夫模型）方法对拷贝数变异（CNV）的预测结果。

包括细胞分组、基因使用、CNV预测的基因和区域、以及贝叶斯网络的MCMC（马尔科夫链蒙特卡罗）推断结果。这些文件是预测结果的核心，包含细胞的拷贝数状态，帮助识别潜在的基因组异常。BayesNetOutput.HMMi6.hmm_mode-samples 目录：贝叶斯网络方法的输出，包括细胞和CNV状态的概率、热图、分组、树状图和细胞的拷贝数状态等。用于展示拷贝数变异的概率估计，帮助可视化细胞群体中的变异分布。

expr.infercnv.* 系列文件：不同分析阶段下的表达数据，用于对拷贝数变异进行进一步分析，以及生成热图和分组信息。HMM_CNV_predictions.* 系列文件：HMM预测的CNV区域和基因信息，用于预测哪些基因或基因区域发生了拷贝数变异。

heatmap_thresholds , observation_groupings , observations_dendrogram：包含生成热图的阈值、细胞分组信息和层次聚类结果。MCMC_inferCNV_obj.rds : 包含MCMC推断结果的rds对象。preliminary.* 系列文件：初步分析阶段生成的数据文件，包括初步的热图、分组和树状图信息。最终分析前用于检查数据处理的效果和初步的拷贝数变异结果。

run.final.infercnv_obj : 最终的inferCNV对象，包含了所有分析步骤的最终结果。我们看一下几个比较重要的输出文件：infercnv.png : 以热图的形式展示 infercnv score，图中可以观察到，作为参考细胞的”Microglia/Macrophage”与”Oligodendrocytes (non-malignant)“，图中几乎没有蓝色(cnv loss)与红色(cnv gain)。

而受试组细胞类型”malignant_MGH36”、“malignant_MGH53”、“malignant_97”、“malignant_93”则存在大量密集的红色(CNV gain)、蓝色区域(CNV loss)。当然，这里作者已经注释出了这些恶性细胞群作为受试细胞，因此它们的CNV事件均比较明显，实际计算下，大家可以通过infercnv的结果来推断未命名亚群的恶性程度。

infercnv.references.txt : 参考细胞的CNV分值矩阵，行名是基因名，列名是细胞名：# 读取数据：ref_score <- read.table ( '/data/09_scRNA-seq_cnv/result_infercnv/test1/infercnv.references.txt' ) ref_score[ 1 : 4 , 1 : 4 ] # 可以看出大部分都在1附近 # MGH53_P2_G04 MGH54_P16_F12 MGH53_P1_D07 MGH53_P11_F08 # WASH7P 1.003744 1.003744 1.264298 1.003744 # LINC00115 1.003744 1.003744 1.261337 1.003744 # NOC2L 1.003744 1.003744 1.258332 1.003744 # SDF4 1.003744 1.003744 1.255518 1.003744 infercnv.observations.txt : 受试细胞的CNV分值矩阵，行名是基因名，列名是细胞名：# 读取数据：obs_score <- read.table ( '/data/09_scRNA-seq_cnv/result_infercnv/test1/infercnv.observations.txt' ) obs_score[ 1 : 4 , 1 : 4 ] # 可以看出相较于ref_score，obs_score与"1"偏离程度更大 # X93_P8_G11 X93_P8_C11 X93_P8_A12 X93_P8_B06 # WASH7P 1.003744 1.003744 0.7897687 0.8224919 # LINC00115 1.003744 1.003744 0.7923169 0.8239459 # NOC2L 1.003744 1.003744 0.7960422 0.8228446 # SDF4 1.003744 1.003744 0.8010698 0.8231211 5、infercnv实战 5.1 数据概况 这里的演示数据来源于转移性肺腺癌的单细胞转录组谱，共包含从44例正常组织或早期到转移期癌症患者的208,506个细胞。

该数据集的原文发表在《Nature communication》中：https://pmc.ncbi.nlm.nih.gov/articles/PMC7210975/。数据可以从GEO数据库中直接下载 GSE131907，当然我也已经给大家下载到本地。

# 读取原始矩阵 mycount <- readRDS ( '/data/09_scRNA-seq_cnv/data/GSE131907_Lung_Cancer_raw_UMI_matrix.rds' ) # 读取注释信息 mymetada <- read.table ( '/data/09_scRNA-seq_cnv/data/GSE131907_Lung_Cancer_cell_annotation.txt' , header = T, row.names = 1 , sep = ' \t ' ) # 创建Seurat对象 scRNA <- CreateSeuratObject ( counts = mycount) # 添加注释信息到Seurat对象中 scRNA <- AddMetaData (scRNA, metadata = mymetada) # 创建文件夹 dir.create ( "/data/09_scRNA-seq_cnv/result_infercnv/test2" , recursive = T) # 保存rds文件 saveRDS (scRNA, "/data/09_scRNA-seq_cnv/result_infercnv/test2/GSE131907_Lung_Cancer_anno.rds" ) library (Seurat) # 看看细胞数：scRNA <- readRDS ( "/data/09_scRNA-seq_cnv/result_infercnv/test2/GSE131907_Lung_Cancer_anno.rds" ) ncol (scRNA) # [1] 208506 # 看看有哪些注释信息：names (scRNA @ meta.data) # [1] "orig.ident" "nCount_RNA" "nFeature_RNA" # [4] "Barcode" "Sample" "Sample_Origin" # [7] "Cell_type" "Cell_type.refined" "Cell_subtype" # 一共有这些样本：unique (scRNA $ Sample_Origin) # [1] "nLN" "mBrain" "nLung" "tLung" "PE" "mLN" "tL/B" # 一共有这些亚群：unique (scRNA $ Cell_type) # [1] "B lymphocytes" "Epithelial cells" "T lymphocytes" # [4] "Myeloid cells" "NK cells" "Undetermined" # [7] "MAST cells" "Fibroblasts" "Endothelial cells" # [10] "Oligodendrocytes" # 取出Epithelial cells用于infercnv分析：Idents (scRNA) <- 'Cell_type' Epi <- subset (scRNA, idents = 'Epithelial cells' ) # 查看亚群分布：table (Epi $ Cell_subtype) # # AT1 AT2 Ciliated Club Malignant cells # 530 2020 654 439 24784 # tS1 tS2 tS3 Undetermined # 3270 3018 64 60 # 取出特定亚群参与分析 Idents (Epi) <- 'Cell_subtype' Epi <- subset (Epi, idents = c ( 'AT1' , 'tS1' , 'Malignant cells' )) # transcriptional states 1(tS1)、tS2、tS3、alveolar types (AT1)、AT2是normal Epithelial的亚群 # 细胞数量还是太多，为了加快演示速度，每个细胞亚群抽500 Epi <- subset (Epi, downsample= 500 ) # 查看当前留存的细胞数：table ( Idents (Epi)) # # Malignant cells AT1 tS1 # 500 接下来就是标准的单细胞分析流程，这部分内容看不懂的同学可以参考scRNA-Seq学习手册Seurat V4修订版 # 标准化 Epi <- NormalizeData (Epi, normalization.method = "LogNormalize" , scale.factor = 10000 ) # 寻找高变基因 Epi <- FindVariableFeatures (Epi, selection.method = "vst" , nfeatures = 2000 ) # 归一化数据 Epi <- ScaleData (Epi, features = VariableFeatures (Epi)) # 计算线粒体含量 Epi[[ "percent.mt" ]] <- PercentageFeatureSet (Epi, pattern = "^MT-" ) Epi <- ScaleData (Epi, vars.to.regress = "percent.mt" ) Epi <- RunPCA (Epi, features = VariableFeatures ( object = Epi)) ElbowPlot (Epi) # 寻找邻近值 Epi <- FindNeighbors (Epi, dims = 1 : 20 ) Epi <- RunUMAP (Epi, dims = 1 : 20 ) Epi <- RunTSNE (Epi, dims = 1 : 20 ) # 展示细胞类型的降维图 DimPlot (Epi, reduction = "tsne" , label = TRUE ) 这里为了演示起来方便，我直接展示的是注释后的细胞标签，以显著的展示出正常细胞与癌症细胞的拷贝数变异事件差异。

实际的分析过程中，你手里拿到的可能是亚群分析的结果，需要通过infercnv的分析才能够完成亚群的注释，如果这部分内容大家理解不了，可以参考下面两个教程：如何分析细胞亚群 Seurat中如何让细胞听你指挥 5.2 infercnv流程 这里我们再走一遍infercnv的基础流程 # 加载包 library (infercnv) library (tidyverse) library (ComplexHeatmap) library (circlize) library (RColorBrewer) library (infercnv) # 整理数据 # 获取count矩阵 exprMatrix <- as.matrix ( GetAssayData (Epi, slot= 'counts' )) # 取出分组注释信息 cellAnnota <- subset (Epi @ meta.data, select= 'Cell_subtype' ) # 把分组注释信息写到本地 write.table (cellAnnota, file = "/data/09_scRNA-seq_cnv/result_infercnv/test2/groupFiles.txt" , sep = ' \t ' , col.names = F) # 制作基因与染色体位置的文件 # 加载注释基因信息的包 if ( ! require (AnnoProbe)) install.packages ( 'AnnoProbe' ) # 获取基因注释信息，这里物种是人，如果是小鼠，把human换成mouse geneInfor = annoGene ( rownames (Epi), "SYMBOL" , 'human' ) # 包含以下信息 colnames (geneInfor) # [1] "SYMBOL" "biotypes" "ENSEMBL" "chr" "start" "end" # 选取基因的 geneInfor = geneInfor[ with (geneInfor, order (chr, start)), c ( 1 , 4 : 6 )] # 去重 geneInfor = geneInfor[ ! duplicated (geneInfor[, 1 ]),] # 查看一共有多少基因 length ( unique (geneInfor[, 1 ])) # [1] 20363 # 查看最终的注释信息 head (geneInfor) # SYMBOL chr start end # 23 AP006222.2 chr1 266855 268655 # 51 FAM87B chr1 817371 819837 # 53 LINC00115 chr1 826206 827522 # 55 FAM41C chr1 868071 876903 # 62 SAMD11 chr1 923928 944581 # 63 NOC2L chr1 944203 959309 # 注释信息需要写到本地供infercnv读取 write.table (geneInfor, '/data/09_scRNA-seq_cnv/result_infercnv/test2/geneLocate.txt' , row.names= F, col.names= F, sep= ' \t ' ) # 创建infercnv对象：infercnv_obj <- CreateInfercnvObject ( raw_counts_matrix= exprMatrix, # 矩阵文件路径 annotations_file= "/data/09_scRNA-seq_cnv/result_infercnv/test2/groupFiles.txt" , # 细胞注释文件路径 delim= " \t " , # 读取文件时的分隔符，如果是csv文件，则用',' gene_order_file= "/data/09_scRNA-seq_cnv/result_infercnv/test2/geneLocate.txt" , # 基因注释文件路径 ref_group_names= "AT1" ) # 正常参考细胞类型的名称，需要根据大家的生物学知识来判断，设定一定未发生肿瘤变异的细胞类型 # 运行infercnv：infercnv_obj <- infercnv :: run (infercnv_obj, # 创建的infercnv对象 cutoff= 0.1 , # 使用Smart-seq2这类全长数据时填cutoff=1, 使用10x Genomics这类3'测序策略数据时填cutoff=0.1 out_dir= '/data/09_scRNA-seq_cnv/result_infercnv/test2' , # 结果输出目录 hclust_method= "ward.D2" , # 指定层次聚类方法，ward.D2最常用的方法，试图最小化类间的方差（推荐用于 scRNA-seq） plot_steps= F, # 如果为真，则保存 infercnv 对象并在中间步骤绘制数据 denoise= TRUE , # 决定是否对数据进行去噪处理 HMM= T, Markov Model, HMM）来进一步细化拷贝数变异的推断。

HMM 会根据推断出的 CNV 模式，进行更精确的基因组段拷贝数变异的检测，帮助识别潜在的基因组异常区域。num_threads= 4 , # 并行计算设置，以4线程为例 cluster_by_groups= TRUE # 决定是否按提供的细胞分组进行聚类。如果 TRUE，细胞会基于之前分配的组信息进行聚类，这有助于分析不同组间的 CNV 变化。如果 FALSE，则不使用组信息进行聚类。

) # 如果出现Error in seq_len(max(obs_annotations_groups)) : argument must be coercible to non-negative integer # 可以参考https://github.com/broadinstitute/infercnv/issues/526进行调整,设置cluster_by_groups=TRUE 可以看到作为参考的 AT1 (上方热图)发生的拷贝数变异事件比较少。

而 Malignant cells 出现了大量且密集的拷贝数变异。tS1 似乎也存在着少量的拷贝数变异，所以单凭热图还不行，还得继续分析。

5.3 热图优化 # 继续可视化 infercnv_obj = readRDS ( "/data/09_scRNA-seq_cnv/result_infercnv/test2/run.final.infercnv_obj" ) # 读取运行好的infercnv对象 expr <- infercnv_obj @ expr.data # 取出infercnv评分数据 # 获取参考/待测细胞的索引 normal_loc <- infercnv_obj @ reference_grouped_cell_indices # 正常参考细胞的索引 normal_loc <- normal_loc $ AT1 test_loc <- infercnv_obj @ observation_grouped_cell_indices # 待测细胞的索引 test_loc <- c (test_loc $ tS1,test_loc $ ` Malignant cells ` ) test_loc[ 1 : 4 ] # [1] 3 6 8 10 # 把之前写到本地的分组信息读回来 anno.df = read.table ( '/data/09_scRNA-seq_cnv/result_infercnv/test2/groupFiles.txt' , header = F, sep = ' \t ' ) colnames (anno.df)[ 1 ] <- 'CB' # 读取并整理此前保存的基因坐标信息：geneFile <- read.table ( "/data/09_scRNA-seq_cnv/result_infercnv/test2/geneLocate.txt" , header = F, sep = " \t " , stringsAsFactors = F) rownames (geneFile) = geneFile $ V1 # 将基因名设置为行名 co_gene <- intersect ( rownames (expr),geneFile $ V1) # 表达矩阵和基因注释文件都存在的基因 sub_geneFile <- geneFile[co_gene,] # 取出基因注释文件的交集 sub_geneFile $ num_chr <- gsub ( 'chr' , '' ,sub_geneFile $ V2) %>% as.numeric # 以数值型进行排序：sub_geneFile <- sub_geneFile[ order (sub_geneFile $ num_chr),] expr = expr[co_gene,] # 取出表达矩阵的交集 set.seed ( 2024 ) kmeans.result <- kmeans ( t (expr), 7 ) # 完成kmeans聚类 kmeans_df <- data.frame ( kmeans_class= kmeans.result $ cluster) # 取出聚类的cluster信息并转化为数据框 kmeans_df $ CB = rownames (kmeans_df) # 设置行名为新的一列，即cell barcode kmeans_df = kmeans_df %>% inner_join (anno.df, by= "CB" ) # 按照barcode合并分组注释文件与聚类结果 kmeans_df_s = arrange (kmeans_df,kmeans_class) # 按照聚类结果进行排序 rownames (kmeans_df_s) = kmeans_df_s $ CB # 设置行名为cell barcode kmeans_df_s $ CB = NULL # 后面用不到cellbarcode了，将这行删除 kmeans_df_s $ kmeans_class = as.factor (kmeans_df_s $ kmeans_class) # 聚类结果是数字，需要转换为因子，否则后面热图注释会报错 colnames (kmeans_df_s)[ 2 ] <- 'Sub_type' # 改变第二列的列名 head (kmeans_df_s) # kmeans_class Sub_type # AAACGGGAGGCGCTCT_NS_12 1 Malignant cells # AATCGGTAGCGTAGTG_NS_12 1 Malignant cells # ACATACGTCGTGGGAA_NS_12 1 Malignant cells # ACCTTTAAGATCACGG_NS_12 1 Malignant cells # ACGAGGAGTAGGGACT_NS_12 1 Malignant cells # ACGAGGATCTGCCCTA_NS_12 1 Malignant cells # 绘制热图 # 热图这部分代码看不懂的同学请先行学习：相关链接见专题页 # 顶部注释 top_anno <- HeatmapAnnotation ( foo = anno_block ( gp = gpar ( fill = 1 : 22 , col = "NA" ), labels = paste0 ( 'Chr' , 1 : 22 ), labels_gp = gpar ( cex = 0.8 , rot = 90 )), annotation_height = unit ( 0.3 , "cm" )) # 调整注释的高度 color_v = RColorBrewer :: brewer.pal ( 8 , "Dark2" )[ 1 : 7 ] names (color_v) = as.character ( 1 : 7 ) # 提取两个子集 malignant_tS1_data <- expr[, rownames (kmeans_df_s[kmeans_df_s $ Sub_type %in% c ( "Malignant cells" , "tS1" ), ])] AT1_data <- expr[, rownames (kmeans_df_s[kmeans_df_s $ Sub_type == "AT1" , ]) ] # 配色方案：my_fn <- fivenum (expr) # 返回最低值、下四分位数、中位数，上四分位数、最大值 col_fun <- colorRamp2 ( c (my_fn[ 1 ], my_fn[ 3 ], my_fn[ 5 ]), c ( " , " , " )) # 确定最低值、中位数、最高值的颜色分配方案 # 为 Malignant cells 和 tS1 生成注释：malignant_tS1_anno <- rowAnnotation ( df = kmeans_df_s[kmeans_df_s $ Sub_type %in% c ( "Malignant cells" , "tS1" ), ], col = list ( Sub_type = c ( "tS1" = " , 'Malignant cells' = ' ), kmeans_class = color_v)) # 为 AT1 生成注释：AT1_anno <- rowAnnotation ( df = kmeans_df_s[kmeans_df_s $ Sub_type == "AT1" , ], col = list ( Sub_type = c ( "AT1" = " ), kmeans_class = color_v)) # 热图1：Malignant cells + tS1的cnvscore热图 ht1 <- Heatmap ( t (malignant_tS1_data)[,sub_geneFile $ V1], col = col_fun, cluster_rows = FALSE , cluster_columns = FALSE , show_column_names = FALSE , show_row_names = FALSE , top_annotation = top_anno, left_annotation = malignant_tS1_anno, column_split = factor (sub_geneFile $ num_chr), border = TRUE , # 添加边框 column_title = NULL ) # `use_raster` is automatically set to TRUE for a matrix with more than # 2000 columns You can control `use_raster` argument by explicitly # setting TRUE/FALSE to it. # # Set `ht_opt$message = FALSE` to turn off this message. # 热图2：AT1的cnv score热图 ht2 <- Heatmap ( t (AT1_data)[,sub_geneFile $ V1], col = col_fun, cluster_rows = FALSE , cluster_columns = FALSE , show_column_names = FALSE , show_row_names = FALSE , left_annotation = AT1_anno, column_split = factor (sub_geneFile $ num_chr), border = TRUE , # 添加边框 show_column_dend = FALSE # 关闭column_split标签 ) # `use_raster` is automatically set to TRUE for a matrix with more than # 2000 columns You can control `use_raster` argument by explicitly # setting TRUE/FALSE to it. # # Set `ht_opt$message = FALSE` to turn off this message. # 把两张热图拼接起来 draw (ht1 %v% ht2, heatmap_legend_side = "right" ) 这样还是可以明显看出Malignant cells与tS1在聚类中分离的比较开，并且发生了更多的拷贝数变异事件。

5.4 小提琴图 5.4.1 从细胞类群的角度看 # 转化cnvscore数据 # fivenum (expr) # 原始的cov_score是以1为中心分布的数据 # [1] 0.7623885 1.0011799 1.0011799 1.0011799 1.5786998 CNV_score = expr -1 # 转换为以0为中心分布的数据 CNV_score = CNV_score ^ 2 # 平方，防止累加的过程中cnv lose和cnv gain相互抵消，造成没有cnv事件发生的假象 CNV_score = as.data.frame ( colMeans (CNV_score)) # 以细胞为单位，获得cnv score colnames (CNV_score) = "CNV_score" # 更改列名 CNV_score $ CB = rownames (CNV_score) # 加入cell barcode列 kmeans_df_s $ CB = rownames (kmeans_df_s) # 把kmeans_df_s中的cell barcode也加上 CNV_score = CNV_score %>% inner_join (kmeans_df_s, by= "CB" ) # 合并cnvscore和聚类结果 library (ggplot2) library (ggsignif) ggplot (CNV_score, aes (Sub_type,CNV_score, fill= Sub_type)) + geom_violin ( alpha= 0.4 ) + stat_boxplot ( geom= "errorbar" , position = position_dodge ( width = 0.1 ), width= 0.1 ) + # 加上误差棒 geom_boxplot ( alpha= 0.5 , outlier.size= 0 , size= 0.3 , width= 0.3 ) + # 加上箱线图 scale_fill_manual ( values = ggsci :: pal_aaas ( 3 )) + # 填充颜色 theme_bw + # 主题为空白 geom_signif ( comparisons = list ( c ( "AT1" , "tS1" ), c ( "tS1" , "Malignant cells" ), c ( "AT1" , "Malignant cells" )), step_increase = 0.1 , map_signif_level = F) 可以看到结果非常完美，作为正常细胞的 AT1 与 tS1 与 Malignant cells之间的统计结果呈现出极显著的关系，而AT1与tS1之间则没有显著性差异(0.086) 5.4.2 从聚类的角度来看 ggplot (CNV_score, aes (kmeans_class,CNV_score, fill= kmeans_class)) + geom_violin ( alpha= 0.4 ) + stat_boxplot ( geom= "errorbar" , position = position_dodge ( width = 0.1 ), width= 0.1 ) + # 加上误差棒 geom_boxplot ( alpha= 0.5 , outlier.size= 0 , size= 0.3 , width= 0.3 ) + # 加上箱线图 scale_fill_manual ( values = ggsci :: pal_nejm ( 8 )) + # 填充颜色 theme_bw # 主题为空白 可以看出1、3、5、6、7明显cnv score要高于其它几个亚群，我们可以通过cluster在Sub_type中的亚群来辅助判断它们的身份：cell.prop <- as.data.frame ( prop.table ( table (CNV_score $ kmeans_class, CNV_score $ Sub_type))) colnames (cell.prop) <- c ( "kmeans_class" , "Sub_type" , "proportion" ) ggplot (cell.prop, aes (kmeans_class,proportion, fill= Sub_type)) + geom_bar ( stat= "identity" , position= "fill" ) + ggtitle ( "" ) + theme_bw + theme ( axis.ticks.length= unit ( 0.5 , 'cm' )) + guides ( fill= guide_legend ( title= NULL )) + theme ( axis.text.x.bottom = element_text ( angle = 45 , hjust = 1 , vjust = 1 )) + scale_fill_manual ( values = ggsci :: pal_nejm ( 8 )) # 填充颜色 果然不出所料，1、3、6、7中均存在很高比例的Malignant cells，而5含有一定比例的Malignant cell，非常的迷惑，假设这是一个未经注释的数据集，我一定会这么命名：CNV_score $ Type <- recode (CNV_score $ kmeans_class, '1' = 'Malignant cells' , '2' = 'Normal Epi' , '3' = 'Malignant cells' , '4' = 'Normal Epi' , '5' = 'Malignant-transition cells' , '6' = 'Malignant cells' , '7' = 'Malignant cells' ) unique (CNV_score $ Type) # [1] Normal Epi Malignant-transition cells # [3] Malignant cells # Levels: Malignant cells Normal Epi Malignant-transition cells # 看一下重命名过的cnv score怎么分布：ggplot (CNV_score, aes (Type,CNV_score, fill= Type)) + geom_violin ( alpha= 0.4 ) + stat_boxplot ( geom= "errorbar" , position = position_dodge ( width = 0.1 ), width= 0.1 ) + # 加上误差棒 geom_boxplot ( alpha= 0.5 , outlier.size= 0 , size= 0.3 , width= 0.3 ) + # 加上箱线图 scale_fill_manual ( values = ggsci :: pal_nejm ( 8 )) + # 填充颜色 theme_bw + # 主题为空白 geom_signif ( comparisons = list ( c ( 'Malignant cells' , 'Normal Epi' ), c ( 'Normal Epi' , 'Malignant-transition cells' ), c ( 'Malignant cells' , 'Malignant-transition cells' )), step_increase = 0.1 , map_signif_level = F) 猜的没错，正常组和另外两组cnv score间均有显著性差异，我们可以将这个注释导回原来的上皮细胞Seurat对象中：rownames (CNV_score) <- CNV_score $ CB # 设置Cell barcode行名，AddMeta会按照行名匹配 Epi <- AddMetaData (Epi, metadata = CNV_score) # 添加CNV score注释信息 Idents (Epi) <- 'Type' # 改变Seurat对象默认标签 DimPlot (Epi, reduction = "tsne" , cols = ggsci :: pal_nejm ( 8 )) # 查看不同type的降维情况 # 下面的工作就是生物学方面的意义啦，比如做个差异分析：DEG <- FindMarkers (Epi, ident.1 = 'Malignant cells' , ident.2 = 'Normal Epi' ) head (DEG) # p_val avg_log2FC pct.1 pct.2 p_val_adj # IGFBP3 1.443336e-146 6.379162 0.604 0.021 4.277182e-142 # GAPDH 1.958204e-116 1.945720 0.994 0.830 5.802941e-112 # LY6D 1.324926e-111 5.484267 0.457 0.011 3.926287e-107 # S100P 3.648096e-110 3.440433 0.569 0.047 1.081077e-105 # FAM83A 2.220502e-106 4.188541 0.569 0.055 6.580235e-102 # STEAP1 4.566941e-102 4.089125 0.466 0.023 1.353367e-97 DEG $ Gene <- row.names (DEG) # 画个火山图看一下：ggplot (DEG, aes (avg_log2FC , - log10 (p_val))) + # 横向水平参考线：geom_hline ( yintercept = - log10 ( 0.05 ), linetype = "dashed" , color = " ) + # 纵向垂直参考线：geom_vline ( xintercept = c ( - 1.2 , 1.2 ), linetype = "dashed" , color = " ) + # 散点图: geom_point ( aes ( size= - log10 (p_val), color= - log10 (p_val))) + # 指定颜色渐变模式：scale_color_gradientn ( values = seq ( 0 , 1 , 0.2 ), colors = c ( " , " , ' , " , " )) + # 指定散点大小渐变模式：scale_size_continuous ( range = c ( 0.5 , 4 )) + # 主题调整：theme_bw + # 调整主题和图例位置：# theme(panel.grid = element_blank, # legend.position = c(0.01,0.7), # legend.justification = c(0,1) # )+ # 设置部分图例不显示：guides ( col = guide_colourbar ( title = "-Log10_q-value" ), size = "none" ) + # 添加标签：geom_text ( aes ( label= Gene, color = - log10 (p_val)), size = 3 , vjust = 1.5 , hjust= 1 ) + # 修改坐标轴：xlab ( "Log2FC" ) + ylab ( "-Log10(FDR q-value)" ) 常规套路就是富集分析之流的，大家可以参考：富集分析 Sys.time # [1] "2025-08-29 17:57:39 CST" save.image ( "result_infercnv/infercnv.rdata" )

五、copykat 1、copykat介绍 1.1 简介 copykat 工具自带的文章在2021年发表在《Nature biotechnology》上，题为Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes。

copykat 是一款基于贝叶斯的工具，作者也在文章对比了其与 infercnv 相比的优势。以前的方法（如 inferCNV和HoneyBadger）可以从足够大的基因组区域的RNA读取计数中估计基因组拷贝数谱。然而，这些方法旨在分析具有较低细胞通量和较高覆盖深度的第一代scRNA-seq技术。

这些方法不适用于分析新开发的高通量scRNA-seq平台（微滴和纳米孔平台），这些平台进行全转录组扩增并在非常稀疏的覆盖深度对mRNA的3’或5’端进行测序。此外，以前的方法无法准确解析特定染色体断点的基因组位置，也无法根据非整倍体拷贝数谱对肿瘤和正常细胞进行分类。区分肿瘤微环境中的正常细胞类型与恶性细胞以及解析肿瘤内的亚克隆在当时比较有挑战性。

为了应对这些挑战，作者开发了copykat并将其应用于各种人类肿瘤，以识别非整倍体肿瘤细胞并描绘肿瘤块内共存的不同亚群的克隆亚结。作者从包含21种肿瘤(包括三阴性乳腺癌、胰腺导管腺癌、甲状腺未分化癌、浸润性导管癌和胶质母细胞瘤)、共 46501 个单细胞的 scRNA-seq 数据中、以5Mb的平均基因组分辨率估计基因组拷贝数图谱。

1.2 copykat工作流程描述 copykat的工作流程将贝叶斯方法与分层聚类相结合，计算得到单个细胞的基因组拷贝数谱。接着从高通量3’scRNA-seq数据中定义亚克隆（Figure 1）。该工作流程将唯一分子标识符的基因表达矩阵作为输入数据。分析从注释基因开始，按基因的坐标进行排序。执行 Freeman-Tukey变换以使方差更加稳定，然后进行多项式动态线性建模让单细胞UMI计数中的异常值更加平滑（Figure 1a）。

接下来通过检测具有高置信度的二倍体细胞亚群，以推断正常 2N细胞的拷贝数基线值。为达到这一目的，作者将单个细胞汇集到几个小的分层集群中，并使用高斯混合模型（GMM）估计每个集群的方差（Figure 1b）。通过遵循严格的分类标准，具有最小估计方差的簇被定义为“置信的二倍体细胞”。当数据只有少数正常细胞，或者当肿瘤细胞具有近二倍体基因组且CNA事件有限时，可能会发生潜在的错误分类。

在这种情况下，copykat提供了一种GMM定义”模式来逐个识别二倍体正常细胞，其中假设单个细胞中基因表达的三个高斯模型的混合物代表基因组增益、损失和中性状态。当处于中性状态的基因占表达基因的至少99%时，单个细胞被定义为“具有置信度的二倍体细胞”。

Figure 1 1.3 性能评估 为了评估copykat的性能，作者通过高通量3’ scRNA-seq (10X Genomics)对来自癌前乳腺肿瘤（DCIS1）的1,480个单肿瘤细胞进行了测序(Figure 2)。作者使用copykat从scRNA-seq数据中计算了全基因组拷贝数谱(Figure 2a-b)，并将结果与inferCNV对相同数据进行的分析结果进行了比较(Figure 2c-d)。

为了获得DNA拷贝数谱的真实实验参考数据，作者从DCIS1中对数百万个非整倍体肿瘤细胞进行了流式分选，用于全基因组bulk DNA测序(Figure 2e-f)。结果表明，copykat在220Kb基因组分辨率下与批量参考DNA 拷贝数谱实现了高度一致性(Pearson 相关性 = 0.82)。

在bulk DNA-seq数据中检测到的大多数主要CNA都在scRNA-seq数据中被鉴定出来，包括chr4、6、8、12、17、20、X染色体出现的拷贝数增加和chr2、3、9、10、15染色体出现的拷贝数丢失(Figure 2a-b)。作者用 inferCNV 处理了同一数据集,并根据成纤维细胞marker 基因ACTA2和FN1手动鉴定基质细胞，这些基因用于提供所需的内部基线参考。

接下来，作者通过使用220Kb可变的分bin将inferCNV的结果转换为与copykat相同的基因组分辨率来比较数据。尽管inferCNV的信号较低，但通过相关性分析，这些数据也与bulk DNA-seq数据参考实现了高度一致(Pearson 相关性 = 0.79)(Figure 2c-d）。

然而，inferCNV的一个主要局限性是它只能报告基因窗口的平滑平均值，并且不检测染色体断点或拷贝数片段的特定坐标，这些内容可以通过copykat来实现。作者通过重复采样不同基因大小区间的相邻局部区域，进一步计算了从两种方法推断的拷贝数状态到参考bulk DNA-seq拷贝数谱的相对距离。

该分析表明，与inferCNV输出的滑动窗口平均值相比，copykat分割结果显著(p 值< 0.001，T检验)更接近参考的DNA拷贝数状态(Figure 2g)。此外，作者发现copykat在5到500个基因的不同大小的基因区间中表现出更稳定的性能(Figure 2h)。

Figure 2 接下来，作者将scRNA-seq结果与’background’即bulk DNA-seq拷贝数谱(Figure S1d)进行比较，评估了copykat从DCIS1数据中检测染色体断点的灵敏度和效率。进行 Bootstrapping 以对单个细胞进行重采样并估计检测灵敏度的变化。

平均而言，我们估计在1Mb分辨率下检测到19%的CNA，在5Mb分辨率下检测到56%，在20Mb分辨率下检测到88%，这与作者对 5Mb平均基因组分辨率的理论计算一致。作者还确定了copykat对分割参数(KS.cut)的敏感程度，该参数被设置为临时修剪截止参数以连接相邻的染色体片段。随着KS.cut值的增加，断点检测变得更加严格，导致检测到的断点更少(Figure S1e)。

我们观察到，当KS.cut超过 0.3(范围0~1)时，分割值的准确性显着下降(Figure S1f)。这些数据表明，copykat可以从高通量3’scRNA-seq数据中以中等基因组分辨率（5Mb）准确推断DNA拷贝数谱。

1.4 实体瘤中肿瘤细胞和正常细胞的分类 作者将copykat应用于先前发表的来自5名胰腺癌(PDAC)患者的3’scRNA-seq数据，以及我们从5名三阴性乳腺癌(TNBC)患者和5名甲状腺未分化癌(ATC)患者生成的3’scRNA-seq数据，以根据肿瘤细胞和正常细胞的拷贝数差异来区分肿瘤细胞和正常细胞(Figure 3)。

作者通过copykat分析了来自5名PDAC患者的9,717个单细胞转录组，并成功鉴定了所有患者的非整倍体肿瘤细胞亚群(Figure 3a & Figure S2a）。预测的肿瘤细胞具有全基因组CNA，包括染色染1q、3q、7p、8q、17、19、20的频繁扩增和染色体3p、6、8p的缺失，这些在 PDAC14-16肿瘤中很常见，而具有二倍体特征的正常细胞没有复发性 CNA。

作者分类的非整倍体肿瘤细胞的UMAP存在共定位现象，这些表达簇使用一组四种已建立的肿瘤上皮标志物(EPCAM、KRT19、KRT18和KRT8)检测到高上皮基因评分。值得注意的是，在所有5例PDAC患者中，仅在两个上皮簇中的一个(c1)中检测到全基因组CNA，表明另一个簇可能是正常的二倍体上皮细胞(c2)，这不能仅通过基因表达来判断。

作者将c1簇指定为肿瘤细胞，因为它们对应于非整倍体拷贝数谱，并且含有较高的KRT19表达水平，KRT19是识别PDAC肿瘤中癌细胞的广泛使用的标志物。作者认为copykat通过计算肿瘤细胞与c1表达簇的共定位，在鉴定肿瘤细胞方面实现了高准确性(98.5% 一致性)(Figure S2)。

根据这些数据，我们估计了肿瘤纯度，范围为6-18%(Figure 3d)，并且与以前的组织病理学数据一致，表明PDAC患者的肿瘤纯度通常较低，因为基质细胞群高18-20。Figure 3 作者还使用10X Genomics平台对来自5个ATC肿瘤的19,568个细胞进行了3’高通量scRNA-seq测序。并使用copykat分析scRNA-seq数据，发现所有5种ATC肿瘤中均存在非整倍体肿瘤细胞(Figure 3b)。

在ATC肿瘤中经常报道的常见CNA：包括在推断的拷贝数数据中鉴定的1p、2p、5p、7、8q、11p、12、18p、20的扩增和在1q、6p、13、17、22中出现的缺失(Figure S2b)。在UMAP表达数据中，预测的非整倍体细胞对应于所有患者的1-2个簇，这些簇显示出较高的上皮细胞评分，包括以前用于组织病理学识别ATC的KRT8。

基于预测的非整倍体肿瘤细胞对c1 (ATC1 中的 c1 和 c2)的共定位，作者估计识别肿瘤细胞的平均预测准确率为97%。作者在ATC患者中观察到广泛的肿瘤纯度范围（2-80%）(Figure 3e)，其中两个肿瘤(ATC2、ATC3)纯度低(2，12%)，三个患者(ATC1、ATC4、ATC5)肿瘤纯度高(42-80%)，与ATC肿瘤病理数据中报告的范围一致。

作者进一步对5名未经治疗的TNBC患者的8,944个单细胞进行了3’scRNA-seq数据进行了分析。在所有5个TNBC样本中，copykat鉴定了具有非整倍体和二倍体拷贝数谱的不同单细胞组。从估计的单细胞拷贝数谱的聚类热图中，我们确定了先前在27-29岁的TNBC患者中报道的频繁 CNA，包括 1q、6p、8、10p、18p、Xq的拷贝数增加和1p、5q、17p、Xp的拷贝数损失(Figure 2c)。

降维显示，预测的非整倍体肿瘤细胞对应于上皮基因评分呈阳性的基因表达簇，其中包括通常用于识别TNBC患者肿瘤细胞的EPCAM (Figure 3c)。在4种肿瘤(TNBC1、TNBC2、TNBC4、TNBC5)中，上皮评分阳性的簇(c1-c2)与copykat预测的非整倍体簇直接对应，然而在一个肿瘤(TNBC3)中，上皮基因组(c1)阳性的3个簇中只有1个显示全基因组CNA。

作者将预测的非整倍体肿瘤细胞与上皮细胞表达簇(TNBC3)中的c1)进行了比较，并估计了(TNBC)肿瘤的高预测准确性(98%)。值得注意的是，在三种病例(TNBC1、TNBC2和TNBC5)中，推断的拷贝数谱定位于两个肿瘤特异性表达簇，这两个簇的上皮基因面板表达均呈阳性，表明肿瘤块中存在多个非整倍体克隆。与PDAC和ATC肿瘤相比，TNBC 样本在患者的肿瘤纯度都很高(34-83%)(Figure 3f）。

总的来说，这些数据表明，copykat可以根据仅从scRNA-seq数据推断的非整倍体拷贝数谱准确(98% ± 3% S.D.)区分各种实体瘤中的肿瘤细胞和正常细胞，而无需特异性基因表达标记(事实上通过基因表达标记准确率堪忧)。

1.5 应用于其他 scRNA-seq 技术 虽然上面的数据表明copykat可以从3’ scRNA-seq数据中准确估计拷贝数数据，但作者进一步研究了copykat是否可以广泛应用于第一代scRNA-seq数据(SMART-Seq2)以及5’ scRNA-seq数据(10X Genomics)。

作者对两个雌激素受体阳性浸润性导管癌(ER+ IDC)肿瘤 (IDC1、IDC2)进行了5’ scRNA-seq，并分析了来自先前发表的SMART-Seq2数据，该研究涉及2 名多形性胶质母细胞瘤(GBM)患者(GBM1、GBM2)。使用5’ scRNA-seq(10X Genomics)从IDC1和IDC2中总共对7,780个单细胞进行了测序。

copykat分析在每个肿瘤中确定了两个簇，代表正常细胞(N) 和肿瘤细胞(T)，其中两个肿瘤都显示出染色体8p的大量扩增(MYC)(Figure 4a、c)。scRNA-seq表达数据的聚类表明，在两种肿瘤中推断的非整倍体肿瘤细胞与显示高上皮评分的簇共定位(Figure 4b、d)，验证了copykat预测的准确性。与TNBC样本类似，ER+ IDC肿瘤均表现出较高的肿瘤细胞结构 (分别为 97% 和 87%)。

接下来，作者分析了先前发表的 scRNA-seq 数据集，该数据集使用来自两名 GBM 患者的第一代全长 scRNA-seq （SMART-seq2） 方法30 进行测序：GBM1 （MGH125）和GBM2 （MGH128）、（GSE131928）。与仅具有部分基因体覆盖的高通量3’或5’scRNA-seq方法（10X Genomics）相比，全长SMART-seq2数据具有覆盖整个基因转录本的序列读数。

然而，SMART-seq2的细胞通量要低得多（每个患者分别为332和184个细胞），并且没有可以减轻扩增偏倚的UMI条形码。为了进行copykat 分析，作者使用了表示归一化基因表达计数数据的TPM/10矩阵。在这两个样本中，我们观察到非整倍体肿瘤细胞和二倍体正常细胞簇之间的明显分离（图 4e、g）。copykat推断的非整倍体肿瘤细胞簇表达高水平的 EGFR （图 4f，h），这是 GBM患者中已确定的肿瘤细胞标志物。

总的来说，这些数据表明 copykat与广泛的scRNA-seq技术兼容。Figure 4 1.6 推断乳腺肿瘤的克隆亚结构 这部分算是作者对copykat的实战检验，为了描绘克隆亚结构并将癌症基因型与表型联系起来，作者将copykat应用于从三名TNBC患者生成的scRNA-seq数据（Figure 5）。对推断的拷贝数谱进行聚类，以根据拷贝数差异识别亚群，并从单细胞谱簇中计算共有拷贝数谱，以识别具有拷贝数差异的基因组区域。

根据亚克隆的共有图谱，作者进行了差异表达（DE）和基因特征分析，以确定亚克隆之间的表型差异。在一个 TNBC 样本（TNBC1）中，797个非整倍体拷贝数谱的聚类确定了两个主要亚克隆（A，B），分别占肿瘤质量的44%和28%（Figure 5a），并在它们的相邻连接（NJ）树中被两个不同的谱系分开。剩下的两个患者也均能通过copykat区分出肿瘤亚克隆。

Figure 5 1.7 小结 作者开发了基于一种综合贝叶斯分割的方法，并封装为成熟的R包copykat，其最大的作用便是在无偏scRNA-seq数据中鉴定肿瘤细胞。正常上皮细胞通常最难仅通过表达谱与恶性肿瘤细胞区分开来，因为它们可以表达许多与癌细胞相同的上皮标志物，而copykat则可以解决这一问题。copykat 的另一个应用是根据拷贝数改变的差异描绘实体瘤中的克隆亚结构，这一应用上一段文字刚描述过，大家可以回看一下。

需要强调的是，copykat能够应用于广泛用于单细胞基因组学领域的3’和5’scRNA-seq 方法（例如10X Genomics、Drop-Seq、InDrop平台）。不可避免的，copykat也存在一定的局限性：并非所有癌症类型都具有可用于区分正常细胞和肿瘤细胞的非整倍体拷贝数事件。特别是，已知儿科癌症和造血癌（例如 AML、CLL）几乎没有拷贝数改变，因此可能不适用于copykat 分析。

另一个限制是copykat主要限于根据整个基因组读取深度的变化检测 CNA事件，不能用于检测其他导致基因组多样性的基因组事件，包括染色体结构重排、插入缺失和体细胞突变。此外，由于3’scRNA-seq数据的技术可变性，copykat无法提供有关具有独特基因型的单个细胞基因组的可靠拷贝数信息。这使得copykat更适合分析许多细胞扩增并具有相似基因型的肿瘤中的亚克隆，而不是复制细胞或极其罕见的亚群的分析。

我们使用copykat注意到的一个潜在问题是，当scRNA-seq数据集没有任何肿瘤细胞时，copykat可能会尝试错误地检测基因表达水平最高的簇中的CNA事件。

2、copykat流程 2.1 安装 copykat 是一个基于R语言的工具包：if ( ! require (devtools)) install.packages (devtools) if ( ! require (copykat)) install_github ( "navinlabcode/copykat" ) 2.2 counts矩阵预处理 # 清空环境变量 rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 9118572 487.0 14930285 797.4 14930285 797.4 # Vcells 16723864 127.6 805612615 6146.4 1006972481 7682.6 # 展示系统时间 Sys.time # [1] "2025-08-29 18:02:58 CST" # 设置工作目录 setwd ( "/data/09_scRNA-seq_cnv/" ) library (Seurat) load ( 'data/exp.rawdata.rda' ) # 作者已经准备好了rdata文件以载入 exp.rawdata[ 1 : 4 , 1 : 4 ] # 这是一个10X数据输出的counts矩阵 # AAACCTGCACCTTGTC AAACGGGAGTCCTCCT AAACGGGTCCAGAGGA # RP11-34P13.3 0 0 0 # FAM138A 0 0 0 # OR4F5 0 0 0 # RP11-34P13.7 0 0 0 # AAAGATGCAGTTTACG # RP11-34P13.3 0 # FAM138A 0 # OR4F5 0 # RP11-34P13.7 0 # 构建Seurat对象 raw <- CreateSeuratObject ( counts = exp.rawdata, project = "copycat.test" , min.cells = 0 , min.features = 0 ) # 可以保存到本地供后续分析 dir.create ( "/data/09_scRNA-seq_cnv/result_copykat/test1" , recursive = TRUE ) write.table (exp.rawdata, file= "result_copykat/test1/exp.rawdata.txt" , sep= " \t " , quote = FALSE , row.names = TRUE ) 2.3 常用参数 上面的矩阵需要保证行名为基因名称而不是基因ID。

copykat的流程同样被封装的很好，所以我们需要了解到分析过程中的一些具体参数：- ngenes.chr 初步过滤细胞时，细胞的矩阵需要满足需要每个染色体上存在至少五个基因，才能够进行后续分析。当然，你也可以通过调整 ngenes.chr=1 去设置这个参数以保留更多的细胞(不推荐)。LOW.DR与与UP.DR 过滤基因通过LOW.DR(默认0.05)与UP.DR(默认0.2)两个比例来控制。LOW.DR：基因在细胞中表达的最小比例。

只有在至少这个比例的细胞中表达的基因才会被保留下来。默认值是0.05，意味着基因至少需要在5%的细胞中有表达才能被保留。这样就保证了过表达的管家基因被剔除，而表达比例较低的特异性基因也被剔除，保证了拷贝数评估的准确性。UP.DR：基因在细胞中表达的最大比例。只有在最多这个比例的细胞中表达的基因才会被保留下来。默认值是0.2，意味着基因最多只能在20%的细胞中有表达才能被保留。

win.size 评分滑窗的基因数 KS.cut 这是一个用于基因组分段的敏感性调节参数，通常用于决定分段（segments）或断点（breakpoints）的数量。它基于K-S检验（Kolmogorov-Smirnov test）来判断是否在某些位置进行分段或分割。KS.cut范围的值通常在0到1之间，常见的设置范围是0.05到0.15。数字越小敏感程度越高。n.cores 并行计算的核心数，默认为1，大于1的数字则开启并行计算。

norm.cell.names 可以通过一个向量输入供 copykat 参考的正常细胞的barcode名称(不是细胞类型)，但是与 infercnv 不同的是 copykat 可以不需要参考细胞。cell.line 如果你的数据是只包含二倍体和非整倍体的细胞系，选”yes”。如果你是组织样本，选”no”。

output.seg 输入布尔值，是否输出seg文件用于IGV浏览器读入 plot.genes 布尔值，是否在热图底部绘制基因名 distance 聚类方法，很难做到一种聚类方法适用所有的数据。因此这里支持距离聚类——“euclidean” 相似性聚类——1-“pearson”、“spearman”等方法。运行起来很简单：注意，作者建议每个样本单独运行，多样本一同运行可能存在批次效应，从而干扰 copykat 的判断。

2.4 copykat library (copykat) # copykat 通过自动生成结果文件并保存在当前工作目录下 # 如果想要指定输出目录，需要先修改工作目录 # setwd("/data/09_scRNA-seq_cnv/result_copykat/test1") copykat.test <- copykat ( rawmat= exp.rawdata, # 表达矩阵 id.type= "S" , # ID类型为symbol ngene.chr= 5 , # 对应细胞的基因至少在每个染色体上包含5个基因，这个细胞才会被包咯iu win.size= 25 , # 评分滑窗的基因数 KS.cut= 0.1 , # 控制敏感度 sam.name= "test1" , # 生成文件名称前缀 distance= "euclidean" , # 聚类的方法 norm.cell.names= "" , # 正常细胞对照，这里为空，即不需要对照 output.seg= "FLASE" , # 是否输出seg文件 plot.genes= "TRUE" , # 是否在热图中绘制基因 genome= "hg20" , # hg20为人类数据，mm10为小鼠数据 n.cores= 4 ) # 开启4核心的并行计算 # [1] "running copykat v1.1.0" # [1] "step1: read and filter data ..." # [1] "33694 genes, 302 cells in raw data" # [1] "12156 genes past LOW.DR filtering" # [1] "step 2: annotations gene coordinates ..." # [1] "start annotation ..." # [1] "step 3: smoothing data with dlm ..." # [1] "step 4: measuring baselines ..." # number of iterations= 315 # number of iterations= 679 # number of iterations= 214 # number of iterations= 823 # number of iterations= 431 # number of iterations= 823 # [1] "step 5: segmentation..." # [1] "step 6: convert to genomic bins..." # [1] "step 7: adjust baseline ..." # [1] "step 8: final prediction ..." # [1] "step 9: saving results..." # [1] "step 10: ploting heatmap ..." # Time difference of 2.919732 mins # 筛选出文件名包含test1的文件 files <- list.files ( "/data/09_scRNA-seq_cnv/" ) test1_files <- files[ grepl ( "test1" , files)] # 移动文件 file.rename ( from = file.path ( "/data/09_scRNA-seq_cnv/" , test1_files), to = file.path ( "/data/09_scRNA-seq_cnv/result_copykat/test1" , test1_files)) # [1] TRUE 这个数据包含300个细胞，计算时间约为3min，大家可以参考一下。

2.5 整理数据 # 拷贝数变异标签：pred.test <- data.frame (copykat.test $ prediction) # 取出确定为非整倍体、二倍体的细胞：pred.test <- pred.test[ which (pred.test $ copykat.pred %in% c ( "aneuploid" , "diploid" )),] defined cells head (pred.test) # cell.names copykat.pred # AAACCTGCACCTTGTC AAACCTGCACCTTGTC aneuploid # AAACGGGAGTCCTCCT AAACGGGAGTCCTCCT diploid # AAACGGGTCCAGAGGA AAACGGGTCCAGAGGA aneuploid # AAAGATGCAGTTTACG AAAGATGCAGTTTACG aneuploid # AAAGCAACAGGAATGC AAAGCAACAGGAATGC aneuploid # AAAGCAATCGGAATCT AAAGCAATCGGAATCT aneuploid # 拷贝数变异矩阵：CNA.test <- data.frame (copykat.test $ CNAmat) CNA.test[ 1 : 4 , 1 : 4 ] # chrompos abspos AAACCTGCACCTTGTC # 1 1 1042457 0.0001093239 # 2 1 1265484 0.0001093239 # 3 1 1519859 0.0001093239 # 4 1 1826619 -0.0249666042 # CNV矩阵的前三行是染色体名称、染色体位置、绝对位置。

每行代表一个220KB大小的bin。

后面则是对应细胞的CNA score 2.6 绘制热图 整体热图绘制：# 热图score的配色：my_palette <- colorRampPalette ( rev (RColorBrewer :: brewer.pal ( n = 3 , name = "RdBu" )))( n = 999 ) # 给染色体分配颜色：chr <- as.numeric (CNA.test $ chrom) %% 2 + 1 # 颜色对应的数字 rbPal1 <- colorRampPalette ( c ( 'black' , 'grey' )) # 2为黑色、1为灰色 CHR <- rbPal1 ( 2 )[ as.numeric (chr)] # 分配颜色 chr1 <- cbind (CHR,CHR) # 结合颜色向量 rbPal5 <- colorRampPalette (RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 2 : 1 ]) com.preN <- pred.test $ copykat.pred # 非整倍体/二倍体标签 pred <- rbPal5 ( 2 )[ as.numeric ( factor (com.preN))] # 按照标签分配注释颜色 cells <- rbind (pred,pred) # 细胞配色，非整倍体为" col_breaks = c ( seq ( - 1 , - 0.4 , length= 50 ), seq ( - 0.4 , - 0.2 , length= 150 ), seq ( - 0.2 , 0.2 , length= 600 ), seq ( 0.2 , 0.4 , length= 150 ), seq ( 0.4 , 1 , length= 50 )) # 绘制热图：heatmap.3 ( t (CNA.test[, 4 : ncol (CNA.test)]), dendrogram= "r" , distfun = function (x) parallelDist :: parDist (x, threads = 4 , # 并行计算距离 method = "euclidean" ), # 指定计算欧式距离 hclustfun = function (x) hclust (x, method= "ward.D2" ), # hclust的方法 ColSideColors= chr1, # 染色体配色 RowSideColors= cells, # 细胞配色 Colv= NA , Rowv= TRUE , notecol= "black" , col= my_palette, breaks= col_breaks, key= TRUE , keysize= 1 , density.info= "none" , trace= "none" , cexRow= 0.1 , cexCol= 0.1 , cex.main= 1 , cex.lab= 0.1 , symm= F, symkey= F, symbreaks= T, cex= 1 , cex.main= 4 , margins= c ( 10 , 10 )) legend ( "topright" , paste ( "pred." , names ( table (com.preN)), sep= "" ), pch= 15 , col= RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 2 : 1 ], cex= 0.6 , bty= "n" ) # copykat的默认配色似乎要更好看一些 当然也可以提取非整倍体单独绘制热图并分出亚群：# 提取aneuploid非整倍体标签：tumor.cells <- pred.test $ cell.names[ which (pred.test $ copykat.pred == "aneuploid" )] # 提取aneuploid非整倍体CNA score：tumor.mat <- CNA.test[, which ( colnames (CNA.test) %in% tumor.cells)] # 设定hclust聚类方法，这里算法与上一张热图相同 hcc <- hclust (parallelDist :: parDist ( t (tumor.mat), threads = 4 , method = "euclidean" ), method = "ward.D2" ) hc.umap <- cutree (hcc, 2 ) # 分配各亚群颜色 rbPal6 <- colorRampPalette (RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 3 : 4 ]) subpop <- rbPal6 ( 2 )[ as.numeric ( factor (hc.umap))] cells <- rbind (subpop,subpop) # 绘制热图：heatmap.3 ( t (tumor.mat), dendrogram= "r" , distfun = function (x) parallelDist :: parDist (x, threads = 4 , method = "euclidean" ), hclustfun = function (x) hclust (x, method= "ward.D2" ), ColSideColors= chr1, RowSideColors= cells, Colv= NA , Rowv= TRUE , notecol= "black" , col= my_palette, breaks= col_breaks, key= TRUE , keysize= 1 , density.info= "none" , trace= "none" , cexRow= 0.1 , cexCol= 0.1 , cex.main= 1 , cex.lab= 0.1 , symm= F, symkey= F, symbreaks= T, cex= 1 , cex.main= 4 , margins= c ( 10 , 10 )) legend ( "topright" , c ( "c1" , "c2" ), pch= 15 , col= RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 3 : 4 ], cex= 0.9 , bty= 'n' ) 3、copykat实战 3.1 数据概况 这里的演示数据来源于转移性肺腺癌的单细胞转录组谱，共包含从44例正常组织或早期到转移期癌症患者的208,506个细胞。

该数据集的原文发表在《Nature communication》中：https://pmc.ncbi.nlm.nih.gov/articles/PMC7210975/。数据可以从GEO数据库中直接下载[GSE131907](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131907 )，当然我也已经给大家下载到本地。

由于copykat建每个样本单独运行，因此这里我们只选取一个样本。

# 创建文件夹 dir.create ( '/data/09_scRNA-seq_cnv/result_copykat/test2/' , recursive = TRUE ) # 读取原始矩阵 mycount <- readRDS ( '/data/09_scRNA-seq_cnv/data/GSE131907_Lung_Cancer_raw_UMI_matrix.rds' ) # 读取注释信息：mymetada <- read.table ( '/data/09_scRNA-seq_cnv/data/GSE131907_Lung_Cancer_cell_annotation.txt' , header = T, row.names = 1 , sep = ' \t ' ) library (Seurat) # 创建Seurat对象：scRNA <- CreateSeuratObject ( counts = mycount) # 添加注释信息到Seurat对象中：scRNA <- AddMetaData (scRNA, metadata = mymetada) # saveRDS(scRNA,"/data/09_scRNA-seq_cnv/result_copykat/test2/GSE131907_Lung_Cancer_anno.rds") # 看看细胞数：ncol (scRNA) # 看看有哪些注释信息：names (scRNA @ meta.data) # 一共有这些样本：unique (scRNA $ Sample_Origin) # 一共有这些亚群：unique (scRNA $ Cell_type) # 取出特定亚群参与分析 Idents (scRNA) <- 'Cell_subtype' scRNA <- subset (scRNA, idents = c ( 'Naive CD4+ T' , 'Malignant cells' )) # 一般情况下Naive CD4+ T不会发生拷贝数变异，这里用作与Malignant cells做对比 # 保留一个样本：scRNA <- scRNA[,scRNA $ Sample == 'BRONCHO_11' ] # saveRDS(scRNA,"/data/09_scRNA-seq_cnv/result_copykat/test2/GSE131907_Lung_Cancer_anno_sub.rds") # 看看细胞数：scRNA <- readRDS ( "/data/09_scRNA-seq_cnv/result_copykat/test2/GSE131907_Lung_Cancer_anno_sub.rds" ) # 查看当前留存的细胞数：table ( Idents (scRNA)) # # Malignant cells Naive CD4+ T # 203 323 3.2 copykat library (copykat) # 如果想要指定输出目录，需要先修改工作目录 # setwd("/data/09_scRNA-seq_cnv/result_copykat/test2") # 注意seurat v4和v5不同的获取counts矩阵的方式 count <- GetAssayData (scRNA, slot = "counts" ) copykat.test <- copykat ( rawmat= count, # 表达矩阵 id.type= "S" , # ID类型为symbol ngene.chr= 5 , # 对应细胞的基因至少在每个染色体上包含5个基因，这个细胞才会被包咯iu win.size= 25 , # 评分滑窗的基因数 KS.cut= 0.1 , # 控制敏感度 sam.name= "test2" , # sample name distance= "euclidean" , # 聚类的方法 norm.cell.names= "" , # 正常细胞对照，这里为空，即不需要对照 output.seg= "FLASE" , # 是否输出seg文件 plot.genes= "TRUE" , # 是否在热图中绘制基因 genome= "hg20" , # hg20为人类数据，mm10为小鼠数据 n.cores= 4 ) # 开启4核心的并行计算 # [1] "running copykat v1.1.0" # [1] "step1: read and filter data ..." # [1] "29634 genes, 526 cells in raw data" # [1] "9495 genes past LOW.DR filtering" # [1] "step 2: annotations gene coordinates ..." # [1] "start annotation ..." # [1] "step 3: smoothing data with dlm ..." # [1] "step 4: measuring baselines ..." # number of iterations= 132 # number of iterations= 111 # number of iterations= 321 # number of iterations= 156 # number of iterations= 368 # [1] "step 5: segmentation..." # [1] "step 6: convert to genomic bins..." # [1] "step 7: adjust baseline ..." # [1] "step 8: final prediction ..." # [1] "step 9: saving results..." # [1] "step 10: ploting heatmap ..." # Time difference of 3.448345 mins # 筛选出文件名包含test2的文件 files <- list.files ( "/data/09_scRNA-seq_cnv/" ) test2_files <- files[ grepl ( "test2" , files)] # 移动文件 file.rename ( from = file.path ( "/data/09_scRNA-seq_cnv/" , test2_files), to = file.path ( "/data/09_scRNA-seq_cnv/result_copykat/test2" , test2_files)) # [1] TRUE 3.3 整理数据 # 拷贝数变异标签：pred.test <- data.frame (copykat.test $ prediction) # 取出确定为非整倍体、二倍体的细胞：pred.test <- pred.test[ which (pred.test $ copykat.pred %in% c ( "aneuploid" , "diploid" )),] defined cells head (pred.test) # cell.names copykat.pred # AAACGGGCAATGGATA_BRONCHO_11 AAACGGGCAATGGATA_BRONCHO_11 diploid # AAATGCCAGGCAGTCA_BRONCHO_11 AAATGCCAGGCAGTCA_BRONCHO_11 diploid # AAATGCCGTTATGCGT_BRONCHO_11 AAATGCCGTTATGCGT_BRONCHO_11 diploid # AACACGTAGCAGATCG_BRONCHO_11 AACACGTAGCAGATCG_BRONCHO_11 diploid # AACCATGCAACGATCT_BRONCHO_11 AACCATGCAACGATCT_BRONCHO_11 aneuploid # AACGTTGCACTTCGAA_BRONCHO_11 AACGTTGCACTTCGAA_BRONCHO_11 diploid # 拷贝数变异矩阵：CNA.test <- data.frame (copykat.test $ CNAmat) CNA.test[ 1 : 4 , 1 : 4 ] # chrompos abspos AAACGGGCAATGGATA_BRONCHO_11 # 1 1 1042457 -0.03289865 # 2 1 1265484 -0.03289865 # 3 1 1519859 -0.03289865 # 4 1 1826619 -0.03289865 # CNV矩阵的前三行是染色体名称、染色体位置、绝对位置。

每行代表一个220KB大小的bin。

后面则是对应细胞的CNA score 3.4 热图绘制 整体热图绘制：# 热图score的配色：my_palette <- colorRampPalette ( rev (RColorBrewer :: brewer.pal ( n = 3 , name = "RdBu" )))( n = 999 ) # 给染色体分配颜色：chr <- as.numeric (CNA.test $ chrom) %% 2 + 1 # 颜色对应的数字 rbPal1 <- colorRampPalette ( c ( 'black' , 'grey' )) # 2为黑色、1为灰色 CHR <- rbPal1 ( 2 )[ as.numeric (chr)] # 分配颜色 chr1 <- cbind (CHR,CHR) # 结合颜色向量 rbPal5 <- colorRampPalette (RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 2 : 1 ]) com.preN <- pred.test $ copykat.pred # 非整倍体/二倍体标签 pred <- rbPal5 ( 2 )[ as.numeric ( factor (com.preN))] # 按照标签分配注释颜色 cells <- rbind (pred,pred) # 细胞配色，非整倍体为" col_breaks = c ( seq ( - 1 , - 0.4 , length= 50 ), seq ( - 0.4 , - 0.2 , length= 150 ), seq ( - 0.2 , 0.2 , length= 600 ), seq ( 0.2 , 0.4 , length= 150 ), seq ( 0.4 , 1 , length= 50 )) # 绘制热图：heatmap.3 ( t (CNA.test[, 4 : ncol (CNA.test)]), dendrogram= "r" , distfun = function (x) parallelDist :: parDist (x, threads = 4 , # 并行计算距离 method = "euclidean" ), # 指定计算欧式距离 hclustfun = function (x) hclust (x, method= "ward.D2" ), # hclust的方法 ColSideColors= chr1, # 染色体配色 RowSideColors= cells, # 细胞配色 Colv= NA , Rowv= TRUE , notecol= "black" , col= my_palette, breaks= col_breaks, key= TRUE , keysize= 1 , density.info= "none" , trace= "none" , cexRow= 0.1 , cexCol= 0.1 , cex.main= 1 , cex.lab= 0.1 , symm= F, symkey= F, symbreaks= T, cex= 1 , cex.main= 4 , margins= c ( 10 , 10 )) legend ( "topright" , paste ( "pred." , names ( table (com.preN)), sep= "" ), pch= 15 , col= RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 2 : 1 ], cex= 0.6 , bty= "n" ) # 可以看到明显的非整倍体和二倍体的热图有所差别 当然也可以提取非整倍体单独绘制热图并分出亚群：# 提取aneuploid非整倍体标签：tumor.cells <- pred.test $ cell.names[ which (pred.test $ copykat.pred == "aneuploid" )] # 提取aneuploid非整倍体CNA score：tumor.mat <- CNA.test[, which ( colnames (CNA.test) %in% tumor.cells)] # 设定hclust聚类方法，这里算法与上一张热图相同 hcc <- hclust (parallelDist :: parDist ( t (tumor.mat), threads = 4 , method = "euclidean" ), method = "ward.D2" ) hc.umap <- cutree (hcc, 2 ) # 分为两个群 # 将分群结果添加回Seurat对象中：scRNA <- AddMetaData (scRNA, as.data.frame (hc.umap)) # 分配各亚群颜色 rbPal6 <- colorRampPalette (RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 3 : 4 ]) subpop <- rbPal6 ( 2 )[ as.numeric ( factor (hc.umap))] cells <- rbind (subpop,subpop) # 绘制热图：heatmap.3 ( t (tumor.mat), dendrogram= "r" , distfun = function (x) parallelDist :: parDist (x, threads = 4 , method = "euclidean" ), hclustfun = function (x) hclust (x, method= "ward.D2" ), ColSideColors= chr1, RowSideColors= cells, Colv= NA , Rowv= TRUE , notecol= "black" , col= my_palette, breaks= col_breaks, key= TRUE , keysize= 1 , density.info= "none" , trace= "none" , cexRow= 0.1 , cexCol= 0.1 , cex.main= 1 , cex.lab= 0.1 , symm= F, symkey= F, symbreaks= T, cex= 1 , cex.main= 4 , margins= c ( 10 , 10 )) legend ( "topright" , c ( "c1" , "c2" ), pch= 15 , col= RColorBrewer :: brewer.pal ( n = 8 , name = "Dark2" )[ 3 : 4 ], cex= 0.9 , bty= 'n' ) 3.5 探究一下拷贝数变异与原来亚群之间的关系 # 把拷贝数变异的结果作为注释文件输入回Seurat对象：scRNA <- AddMetaData (scRNA, metadata = pred.test) library (ggplot2) cell.prop <- as.data.frame ( prop.table ( table (scRNA $ Cell_subtype, scRNA $ copykat.pred))) colnames (cell.prop) <- c ( "Cell_subtype" , "CNV" , "proportion" ) ggplot (cell.prop, aes (Cell_subtype,proportion, fill= CNV)) + geom_bar ( stat= "identity" , position= "fill" ) + ggtitle ( "" ) + theme_bw + theme ( axis.ticks.length= unit ( 0.5 , 'cm' )) + guides ( fill= guide_legend ( title= NULL )) + theme ( axis.text.x.bottom = element_text ( angle = 45 , hjust = 1 , vjust = 1 )) 非常得明显，只有Malignant cells中存在拷贝数变异的细胞，而Naive CD4+ T中一个出现拷贝数变异的细胞都没有，这是在没有设置参考细胞下得到的哦！

3.6 CNA score小提琴图 # 转化cnascore数据 # expr <- CNA.test[, 4 : ncol (CNA.test)] # 前三列是坐标数据 CNV_score = expr # 以0为中心分布的数据 CNV_score = CNV_score ^ 2 # 平方，防止累加的过程中cnv lose和cnv gain相互抵消，造成没有cnv事件发生的假象 CNV_score = as.data.frame ( colMeans (CNV_score)) # 以细胞为单位，获得cnv score colnames (CNV_score) = "CNV_score" # 更改列名 CNV_score $ CB = rownames (CNV_score) # 加入cell barcode列 scRNA $ CB <- colnames (scRNA) # 统一cell barcode名称 CNV_score = CNV_score %>% inner_join (scRNA @ meta.data, by= "CB" ) # 合并cnvscore和聚类结果 library (ggplot2) library (ggsignif) # 按照CNA变异情况分类：ggplot (CNV_score, aes (copykat.pred,CNV_score, fill= copykat.pred)) + geom_violin ( alpha= 0.4 ) + stat_boxplot ( geom= "errorbar" , position = position_dodge ( width = 0.1 ), width= 0.1 ) + # 加上误差棒 geom_boxplot ( alpha= 0.5 , outlier.size= 0 , size= 0.3 , width= 0.3 ) + # 加上箱线图 scale_fill_manual ( values = ggsci :: pal_aaas ( 3 )) + # 填充颜色 theme_bw + # 主题为空白 geom_signif ( comparisons = list ( c ( "aneuploid" , "diploid" ) ), step_increase = 0.1 , map_signif_level = F) # 按照细胞亚群分类：ggplot (CNV_score, aes (Cell_subtype,CNV_score, fill= Cell_subtype)) + geom_violin ( alpha= 0.4 ) + stat_boxplot ( geom= "errorbar" , position = position_dodge ( width = 0.1 ), width= 0.1 ) + # 加上误差棒 geom_boxplot ( alpha= 0.5 , outlier.size= 0 , size= 0.3 , width= 0.3 ) + # 加上箱线图 scale_fill_manual ( values = ggsci :: pal_aaas ( 3 )) + # 填充颜色 theme_bw + # 主题为空白 geom_signif ( comparisons = list ( c ( "Naive CD4+ T" , "Malignant cells" )), step_increase = 0.1 , map_signif_level = F) Sys.time # [1] "2025-08-29 18:10:31 CST" save.image ( "/data/09_scRNA-seq_cnv/result_copykat/copykat.rdata" )

六、SCEVAN是一个R包，它从scRNA-seq数据的原始计数矩阵开始，通过将肿瘤微环境的非恶性细胞与恶性细胞分离来自动对活检中存在的细胞进行分类，并表征这些恶性细胞的克隆结构。它识别具有不同拷贝数结构的细胞亚群，并报告每个亚群的特异性和共享改变。该工具的目的是通过允许以非常简单且完全无监督的方式执行整个分析来自动化整个分析。

1、简介 DOI: 10.1038/s41467-023-36790-9 链接：https://pubmed.ncbi.nlm.nih.gov/36841879/ 1.1 工作流程 SCEVAN算法是一个旨在从原始scRNA-seq数据中自动解析肿瘤克隆亚结构的“端到端”计算管道（图1）。整个流程始于输入的原始基因表达计数矩阵。

在对低质量细胞和基因进行初步过滤后，算法首先会自动识别出一小部分可信度高的非恶性细胞，并将它们的表达谱作为基准，计算出所有细胞相对于这个正常基线的基因表达差异矩阵。为了在降噪的同时保留真实的拷贝数变化边界，该差异矩阵会经过一个边缘保留非线性扩散滤波器（edge-preserving nonlinear diffusion filter）进行平滑处理。处理后的平滑矩阵是后续分析的核心。

SCEVAN会采用其关键创新(一个源于计算机视觉领域的联合分割算法(joint segmentation algorithm))对所有细胞的表达谱进行同步分割。该算法基于一个核心假设：同一克隆内的所有细胞共享相同的拷贝数断点，因此可以整合所有细胞的信号来更准确地识别这些断点。分割完成后，算法会自动将所有细胞分为恶性与非恶性两大类，其依据是包含初始可信正常细胞最多的细胞簇被定义为非恶性群体。

随后，仅针对恶性细胞群体，SCEVAN会进一步通过聚类来识别不同的亚克隆。每个亚克隆的拷贝数谱会再经过一次独立的精细分割，最终将各基因组区域划分为“缺失”、“丢失”、“中性”、“增加”或“扩增”五种拷贝数状态。流程的最后，SCEVAN会进行全面的下游分析，包括区分主干、共享和克隆特异性的拷贝数变异，并构建克隆演化树，从而完整地描绘出肿瘤的克隆结构和演化关系。

1.2 合成数据上进行恶性细胞分类 研究团队希望在可控的条件下定量评估SCEVAN算法的准确性，于是他们在模拟数据（synthetic data）上进行了一系列严格的基准测试。利用scDesign2工具，一个基于真实的多发性骨髓瘤数据集，生成了500个具有已知细胞身份的模拟scRNA-seq数据集。

这些数据集被设计用来模拟真实的肿瘤样本，包含了不同的肿瘤纯度和细胞组成，并且为了进一步增加挑战性，研究者还使用SPLATTER工具向数据中人为添加了三种不同水平的dropout噪声。在对这总计超过32万个模拟细胞的测试中，SCEVAN的表现全面优于当前先进的CopyKAT算法。研究人员使用F1分数作为准确性的衡量标准，结果显示SCEVAN在所有噪声水平下都取得了显著更高的分数。

例如，在无额外噪声的情况下，SCEVAN的F1分数为0.948，而CopyKAT为0.798；在最高噪声水平下，SCEVAN的F1分数仍有0.824，而CopyKAT则降至0.726。这一系列在模拟数据上的测试结果，有力地证明了SCEVAN在自动识别恶性细胞方面具有卓越的准确性和鲁棒性。1.3 真实数据中的恶性细胞分类性能验证 作者将SCEVAN和CopyKAT应用到了一个包含106个真实肿瘤样本的大型公开数据集中进行现实场景验证。

这些样本覆盖了胶质母细胞瘤（GBM）、头颈癌（HNSCC）和结直肠癌三种不同的癌症类型和多种测序技术，总细胞数超过93,000个。在此次测试中，作者以原始研究中经过专家手动注释的细胞身份作为“金标准”，来评判两种自动化算法的准确性。结果再次显示了SCEVAN的优越性（图2）：它在63%的样本中取得了比CopyKAT更高的分类分数，而CopyKAT仅在23%的样本中占优。

从总体性能上看，SCEVAN在所有样本中的总平均F1分数达到了0.90，远高于CopyKAT的0.63。作者也指出，在肿瘤细胞数量极少（1-15个）的样本中，SCEVAN的准确性会有所下降，但这一系列在多样化真实数据集上的测试结果，仍旧有力地证明了SCEVAN能够准确地从scRNA-seq数据中区分肿瘤细胞和正常细胞。

1.4 合成数据中的拷贝数断点分割性能验证 在证明其卓越的细胞分类能力后，作者进一步对SCEVAN的核心技术之变分联合分割算法（variational joint segmentation algorithm）进行验证。在识别拷贝数变异（CNV）断点上的准确性进行了深入验证。研究团队在模拟数据（synthetic data）上进行了测试以便于进行可控的定量评估。

他们构建了包含克隆性和亚克隆性拷贝数变异的复杂数据集，其中变异的位置和强度都是已知的。在与CopyKAT及其他分割算法的直接比较中，SCEVAN在所有测试场景下都展现出显著更高的准确性。无论是在F1分数还是在精确率-召回率曲线下面积（AUC）这两个关键指标上，SCEVAN的表现都全面胜出，这表明其联合分割策略能更有效地从嘈杂的单细胞表达数据中识别出真实的拷贝数边界。

1.5 参考数据的分割准确性验证 那么如何在真实世界数据中验证其核心分割算法的准确性？作者使用了一批同时拥有scRNA-seq数据和作为“金标准”的群体DNA测序（WGS）数据的肿瘤样本。通过比较各种算法从scRNA-seq中推断出的CNV谱与DNA测序得到的真实CNV谱之间的皮尔逊相关性（Pearson correlation），来评估算法的准确性。

结果再次证实了SCEVAN的优越性（图3）：在胶质母细胞瘤（GBM）和多发性骨髓瘤（MM）两个数据集中，SCEVAN推断出的CNV谱与“金标准”的相关性均显著高于inferCNV和CopyKAT。例如，在GBM数据集中，SCEVAN的平均相关性为0.57，而inferCNV为0.44，CopyKAT仅为-0.03。此外，作者还通过一系列对照实验证明，SCEVAN的这种优势主要来源于其更优秀的分割算法，而非仅仅是细胞分类的准确性。

并且，该算法对于参考正常细胞的轻微错误分类表现出高度的鲁棒性（robustness），这在处理复杂的真实世界样本时是一个巨大的优势。1.6. 计算效率评估 除了准确性，计算效率也是衡量一个算法实用性的关键指标，尤其是在处理日益庞大的单细胞数据集时。研究团队对SCEVAN算法的运行速度进行了评估，并与CopyKAT和inferCNV进行了比较。

他们指出，SCEVAN之所以高效，关键在于其核心的分割步骤采用了一种贪心区域增长算法（greedy region-growing algorithm）。在对算法不同阶段的计时测试中，SCEVAN在初始的恶性细胞分类步骤中，速度是竞争对手的2至7倍。在更为耗时的分割步骤中，SCEVAN的优势更为显著：在一个胶质母细胞瘤（GBM）数据集上，它比CopyKAT快2倍，比inferCNV快5倍。

特别是在处理细胞数量庞大的10x Genomics数据集（如多发性骨髓瘤样本）时，SCEVAN的效率优势被进一步放大，其速度分别比inferCNV和CopyKAT快了11倍和19倍。这些结果有力地证明，SCEVAN不仅在精度上超越了同类工具，在计算效率上也表现卓越，使其特别适合应用于当前的大规模单细胞数据分析。

1.7 胶质母细胞瘤中的应用：揭示肿瘤内部的克隆异质性 作者将其应用于分析胶质母细胞瘤（Glioblastoma, GBM）这一具有高度瘤内异质性（intratumoral heterogeneity）的恶性脑瘤，以便于展示SCEVAN在真实生物学问题中的应用价值。

在一个名为MGH105的样本中，SCEVAN成功地从scRNA-seq数据中识别出了四个亚克隆（four subclones），而这些亚克隆是标准scRNA-seq分析流程无法分辨的，其存在此前仅通过单细胞DNA甲基化这一更直接的基因组学技术才得以证实，这突显了SCEVAN算法的超高分辨率。在另一个名为BT1160的样本中，SCEVAN不仅识别出三个亚克隆并构建了它们的演化树（phylogenetic tree）（图4）。

更重要的是，它还将这种由拷贝数变异（CNA）定义的遗传异质性与功能异质性关联了起来：分析显示，这三个亚克隆分别对应了GBM的不同代谢亚型（神经元型、线粒体型和增殖/祖细胞型）。此外，通过对增殖型亚克隆特有的基因扩增区域进行分析，算法还锁定了一个潜在的驱动基因UBE2T。最后，分析也揭示了关键抑癌基因（如PTEN和CDKN2A）的缺失是亚克隆级别的，即只存在于部分而非全部的肿瘤细胞中（图5）。

这些应用实例充分证明，SCEVAN不仅能精确地识别肿瘤亚克隆，还能将遗传结构与细胞功能状态联系起来，从而从scRNA-seq数据中发掘出具有重要生物学和临床意义的见解。1.8 多区域胶质母细胞瘤的克隆演化 作者还将其应用于一个多区域胶质母细胞瘤（GBM）样本的分析中，这样可以进一步展示SCEVAN在解析复杂肿瘤生物学问题上的能力。

由于单个肿瘤活检可能无法完全代表整个肿瘤的异质性，因此对来自同一肿瘤不同空间位置的样本进行分析，这对于理解肿瘤演化至关重要。在本文中，作者分析了一位患者的七个活检样本，其中两个来自肿瘤外周（tumor periphery），五个来自肿瘤核心（tumor core）。通过对每个活检样本进行独立的克隆分析，SCEVAN成功地重构了该肿瘤内部复杂的空间克隆结构，并推断出了一个演化树（evolutionary tree）（图6）。

该演化树清晰地显示，来自肿瘤外周的克隆形成了一个独立的演化分支，它们缺少存在于核心区域克隆中的4号和8号染色体扩增。此外，一个位于2号染色体上的扩增事件在外周区域是克隆性的（所有肿瘤细胞都存在），但在部分核心区域则转变为亚克隆性的（仅部分肿瘤细胞存在），这揭示了肿瘤在空间扩张过程中的动态演化过程。这一应用实例表明，SCEVAN能够有效整合多区域scRNA-seq数据，从而描绘出肿瘤的地理克隆图谱并推断其演化路径。

1.9 原发灶与转移淋巴结的克隆结构 作者利用该算法比较了头颈癌（HNSCC）患者的原发肿瘤（primary tumor）与其对应的淋巴结转移灶（lymph node metastasis）之间的克隆结构关系。在分析的多个病例中，大部分情况下转移灶的克隆结构与原发灶高度相似，显示出很高的相关性，这表明转移细胞直接来源于原发灶的主要克隆。

然而，在一个名为HNSCC5的特定病例中，SCEVAN揭示了一个重要的演化差异：原发灶中存在一个7号染色体的扩增，而这个扩增在淋巴结转移灶中却消失了（图7）。进一步分析发现，这个区域包含一个名为GPNMB的关键基因，该基因在转移灶中也相应地表现出下调。作者指出，GPNMB在其他癌症中通常与促进肿瘤生长和转移有关，因此它在转移过程中的丢失可能是一个值得深入研究的演化事件。

这个案例充分说明，SCEVAN可以被用来精确比较原发灶和转移灶的克隆差异，从而为研究肿瘤转移过程中的克隆演化和筛选过程提供有力的工具。1.10 小结 综上所述，SCEVAN算法通过变分方法成功从scRNA-seq数据中解卷肿瘤克隆子结构，准确区分恶性细胞，并揭示了肿瘤异质性、克隆演化和微环境交互的关键机制。在脑肿瘤应用中，它识别了如UBE2T等潜在治疗靶点，并阐明了PTEN和CDKN2A等抑癌基因的克隆状态。

然而，该算法的局限在于依赖非整倍体假设，不适合如白血病或儿童癌症等基因组变异较少的肿瘤类型。未来研究可扩展到更多技术平台和数据集，结合多组学整合，进一步验证其在临床转化中的鲁棒性。总体而言，SCEVAN为肿瘤异质性研究注入了新活力，有望推动精准医学的进步。

2、安装R包 library (devtools) install_github ( "miccec/yaGST" ) install_github ( "AntonioDeFalco/SCEVAN" ) 3、加载R包 # 清空环境变量 rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 9781130 522.4 14930285 797.4 14930285 797.4 # Vcells 17998057 137.4 644490092 4917.1 1006972481 7682.6 # 展示系统时间 Sys.time # [1] "2025-08-29 18:10:40 CST" library (SCEVAN) # 设置工作目录 setwd ( "/data/09_scRNA-seq_cnv/" ) 4、单样本分析 从胶质母细胞公共数据集 (GSE131928) 加载MGH106样本的scRNA-seq数据。

单个调用（pipelineCNA）能够执行整个关于克隆结构分类和表征的分析。+ count_mtx：具有基因行（允许使用基因符号或 Ensembl ID）和细胞列（列中包含细胞名称）的计数矩阵。

+ sample：保存结果的样本名称（可选） + par_cores：运行管道所需的内核数量（可选 - 默认值为 20） + norm_cells：用于作为可靠正常细胞的可能已知正常细胞向量（可选） + FIXED_NORMAL_CELLS：如果 norm_cell 向量用于作为固定参考，则为 TRUE，如果您仅对克隆结构感兴趣而不关心正常/肿瘤分类，则为 FALSE（默认值为 FALSE） + SUBCLONES：如果对克隆结构感兴趣，则为 TRUE，否则仅为恶性和非恶性细胞分类感兴趣，则为 FALSE（可选 - 默认值为 TRUE） + beta_vega：指定用于分割的 beta 参数，beta 值越高分割越粗粒化。

（可选 - 默认值为 0.5） + ClonalCN：从所有肿瘤细胞中获取克隆 CN 分析结果（可选） + plotTree：绘制系统发育树（可选 - 默认值为 FALSE） + AdditionalGeneSets：正常细胞类型附加签名的列表（可选） + SCEVANsignatures：如果只想使用 AdditionalGeneSets 中指定的签名，则为 FALSE（可选 - 默认值为 TRUE） + organism : 要进行分析的生物体（可选 - “mouse” 或 “human” - 默认为 “human”） 4.1 加载数据 # 示例：results <- pipelineCNA(count_mtx) # 加载数据 load ( "/data/09_scRNA-seq_cnv/data/MGH106_data.RData" ) 4.2 执行分析 single_result <- SCEVAN :: pipelineCNA ( count_mtx, # 单样本每个基因在每个细胞的表达矩阵 sample = "MGH106" , # 保存结果的样本名称 par_cores = 4 , # 分析所用核心数 SUBCLONES = TRUE , # 如果对克隆结构感兴趣，则为 TRUE，否则仅为恶性和非恶性细胞分类感兴趣，则为 FALSE（可选 - 默认值为 TRUE） plotTree = TRUE ) # 绘制系统发育树（可选 - 默认值为 FALSE） # 如果上一步出现报错提示：Error in f(...) : Graphics API version mismatch # 可以尝试在Rstudio 运行 install.packages("ragg") # 但是这时候可能还是会因为一些可能的权限原因提示无法安装 # 可以去linux端，执行sudo R ,启用管理员权限后再执行install.packages("ragg")，将ragg包的版本更新至‘1.4.0’ # 默认保存到当前工作目录下面的output文件夹，可以挪动结果目录到指定位置 unlink ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/" , recursive = TRUE ) dir.create ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample" , recursive = TRUE ) file.rename ( "output/" , "result_SCEVAN/single_sample/" ) 4.3 结果查看 # 返回结果包含每个细胞的分类（肿瘤/正常），它是否被用作确信的正常细胞，以及它属于哪个克隆亚群 head (single_result) # class confidentNormal subclone # MGH106-P7-E10 tumor <NA> 1 # MGH106-P2-B01 tumor <NA> 1 # MGH106-P2-B12 tumor <NA> 1 # MGH106-P2-H07 tumor <NA> 2 # MGH106-P2-E09 tumor <NA> 1 # MGH106-P2-E04 tumor <NA> 1 dim (single_result) # [1] 200 3 # 正常和肿瘤细胞数目 table (single_result $ class) # # normal tumor # 64 136 # 被用作确信的正常的正常细胞数目 table (single_result $ confidentNormal) # # yes # 30 # 肿瘤细胞每个克隆亚型的数目 table (single_result $ subclone) # # 1 2 3 # 58 57 21 拷贝数变异矩阵的热图，包含非恶性和恶性细胞的分类。

文件路径为：/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106heatmap.png 肿瘤细胞克隆亚群拷贝数变异矩阵的热图 (heatmap_subclones.png)，文件路径为：/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106heatmap_subclones.png 从亚克隆档案推断出的克隆树，文件路径为：/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106CloneTree.png 每个亚群体中存在的变化的紧凑共识图，文件路径为：/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106consensus.png OncoPrint样图突出显示特定改变、亚克隆之间的共享改变或克隆改变。

这些片段分为五种预定义的拷贝数状态之一：缺失、丢失、中性、增益或扩增。文件路径为：/data/09_scRNA-seq_cnv/img/MGH106OncoHeat2.png 特定改变的 DE 分析 (DEchr*_subclones.png)。火山图是从属于特定改变的基因的差异表达分析中获得的。

# 展示出指定目录下所有文件名称包含了MGH106-DEchr字符的文件 files <- list.files ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/" , pattern = "MGH106-DEchr" ) files # [1] "MGH106-DEchr1-155687960-160286348_subclones.png" # [2] "MGH106-DEchr1-28736621-36305357_subclones.png" # [3] "MGH106-DEchr1-74198212-92841924_subclones.png" # [4] "MGH106-DEchr11-87037844-118925608_subclones.png" # [5] "MGH106-DEchr11-95768942-117204782_subclones.png" # [6] "MGH106-DEchr12-109116595-120534434_subclones.png" # [7] "MGH106-DEchr12-280129-6961316_subclones.png" # [8] "MGH106-DEchr12-28133249-75511636_subclones.png" # [9] "MGH106-DEchr12-30709552-76853696_subclones.png" # [10] "MGH106-DEchr13-90060247-102876001_subclones.png" # [11] "MGH106-DEchr15-45587123-55508234_subclones.png" # [12] "MGH106-DEchr15-81159575-99135593_subclones.png" # [13] "MGH106-DEchr16-46955362-58283836_subclones.png" # [14] "MGH106-DEchr16-51135975-65126112_subclones.png" # [15] "MGH106-DEchr19-47713343-54102692_subclones.png" # [16] "MGH106-DEchr2-217730-27988087_subclones.png" # [17] "MGH106-DEchr2-46580937-74473322_subclones.png" # [18] "MGH106-DEchr21-17788967-46665124_subclones.png" # [19] "MGH106-DEchr21-44328944-46665124_subclones.png" # [20] "MGH106-DEchr22-29788608-49827512_subclones.png" # [21] "MGH106-DEchr3-113532296-139357223_subclones.png" # [22] "MGH106-DEchr3-183487531-198043720_subclones.png" # [23] "MGH106-DEchr4-108650584-127840733_subclones.png" # [24] "MGH106-DEchr4-153152342-183040119_subclones.png" # [25] "MGH106-DEchr4-68646630-77819615_subclones.png" # [26] "MGH106-DEchr8-73290267-93741017_subclones.png" # [27] "MGH106-DEchr9-88388419-98015943_subclones.png" # 提取第一个文件路径 first_file <- paste0 ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/" ,files[ 1 ]) 例如我这里展示第一个文件 亚克隆的富集分析 (pathway Analysis_subclones*.png)，利用GSEA获得的每个亚克隆的REACTOME途径活性。

# 展示出指定目录下所有文件名称包含了MGH106pathwayAnalysis字符的文件 files <- list.files ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/" , pattern = "MGH106pathwayAnalysis" ) files # [1] "MGH106pathwayAnalysis_subclones1.png" # [2] "MGH106pathwayAnalysis_subclones2.png" # [3] "MGH106pathwayAnalysis_subclones3.png" # 提取第一个文件路径 first_file <- paste0 ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/" ,files[ 1 ]) 这里同样以第一个文件为例 5、多样本分析 从胶质母细胞公共数据集 (GSE131928) 加载了三个样本 MGH102、MGH104 和 MGH105的scRNA-seq数据。

单个调用（multiSampleComparisonClonalCN）允许比较多个样品的克隆谱。

+ listCountMtx：要分析的样品的原始计数矩阵的命名列表 + listNormCells：要分析的每个样品的正常细胞列表 + analysisName：分析的名称（可选） + organism：可选“鼠”或“人”，默认为“人” + par_cores：内核数（默认为20） 5.1 加载数据 # 示例：multiSampleComparisonClonalCN(listCountMtx) # load(url("https://www.dropbox.com/s/esqvnltucdqajg1/listCountMtx.RData?raw=1")) # 此步骤可能出现读取失败的情况，因此我已经下载到了本地，可以直接从本地加载 load ( "/data/09_scRNA-seq_cnv/data/listCountMtx.RData" ) 5.2 执行分析 # 运行多样本分析步骤,约10min multi_result <- SCEVAN :: multiSampleComparisonClonalCN ( listCountMtx, # 多样本表达矩阵list analysisName = "all" , # 分析的名称 organism = "human" , # 分析物种 par_cores = 4 ) # 核心数 # [1] " raw data - genes: 23686 cells: 355" # [1] "1) Filter: cells > 200 genes" # [1] "2) Filter: genes > 10% of cells" # [1] "12813 genes past filtering" # [1] "3) Annotations gene coordinates" # [1] "found 0 confident non malignant cells" # [1] "11490 genes annotated" # [1] "4) Filter: genes involved in the cell cycle" # [1] "10928 genes past filtering " # [1] "5) Filter: cells > 5genes per chromosome " # [1] "6) Log Freeman Turkey transformation" # [1] "A total of 355 cells, 10928 genes after preprocessing" # [1] "7) Measuring baselines (pure tumor - synthetic normal cells)" # [1] "8) Smoothing data" # [1] "10) Adjust baseline" # [1] "11) plot heatmap" # [1] "found 355 tumor cells" # [1] "time classify tumor cells: 3.66527685721715" # [1] " raw data - genes: 23686 cells: 634" # [1] "1) Filter: cells > 200 genes" # [1] "2) Filter: genes > 10% of cells" # [1] "12189 genes past filtering" # [1] "3) Annotations gene coordinates" # [1] "found 30 confident non malignant cells" # [1] "10986 genes annotated" # [1] "4) Filter: genes involved in the cell cycle" # [1] "10481 genes past filtering " # [1] "5) Filter: cells > 5genes per chromosome " # [1] "6) Log Freeman Turkey transformation" # [1] "A total of 634 cells, 10481 genes after preprocessing" # [1] "7) Measuring baselines (confident normal cells)" # [1] "8) Smoothing data" # [1] "9) Segmentation (VegaMC)" # [1] "10) Adjust baseline" # [1] "11) plot heatmap" # [1] "found 427 tumor cells" # [1] "time classify tumor cells: 5.70456975301107" # [1] " raw data - genes: 23686 cells: 312" # [1] "1) Filter: cells > 200 genes" # [1] "2) Filter: genes > 10% of cells" # [1] "12359 genes past filtering" # [1] "3) Annotations gene coordinates" # [1] "found 30 confident non malignant cells" # [1] "11153 genes annotated" # [1] "4) Filter: genes involved in the cell cycle" # [1] "10641 genes past filtering " # [1] "5) Filter: cells > 5genes per chromosome " # [1] "6) Log Freeman Turkey transformation" # [1] "A total of 312 cells, 10641 genes after preprocessing" # [1] "7) Measuring baselines (confident normal cells)" # [1] "8) Smoothing data" # [1] "9) Segmentation (VegaMC)" # [1] "10) Adjust baseline" # [1] "11) plot heatmap" # [1] "found 216 tumor cells" # [1] "time classify tumor cells: 2.71259829600652" # [1] "Phylogenic tree cannot be plotted (ggtree requires dplyr version 1.0.5) try downgrading dprly to 1.0.5 or upgrading ggtree" # 同样如果上一步出现报错提示：Error in f(...) : Graphics API version mismatch # 可以尝试在Rstudio 运行 install.packages("ragg") # 但是这时候可能还是会因为一些可能的权限原因提示无法安装 # 可以去linux端，执行sudo R ,启用管理员权限后再执行install.packages("ragg")，将ragg包的版本更新至‘1.4.0’ # 默认保存到当前工作目录下面的output文件夹，可以挪动结果目录到指定位置 unlink ( "/data/09_scRNA-seq_cnv/result_SCEVAN/multi_sample" , recursive = TRUE ) dir.create ( "/data/09_scRNA-seq_cnv/result_SCEVAN/multi_sample" , recursive = TRUE ) file.rename ( "output" , "result_SCEVAN/multi_sample" ) # [1] TRUE 5.3 结果查看 class (single_result) # [1] "data.frame" # 与单样本不同的是，多样本分析生成的结果对象是一个list,总共包含了两个部分 class (multi_result) # [1] "list" # 包含了三个样本 names (multi_result[[ 1 ]]) # [1] "MGH104" "MGH105" "MGH102" # 每个样本包含了鉴定的正常/肿瘤细胞信息，以及当前样本细胞它是否被用作确信的正常细胞 head (multi_result[[ 1 ]][[ 1 ]]) # class confidentNormal # MGH104-P5-A02 tumor <NA> # MGH104-P5-A04 tumor <NA> # MGH104-P5-A05 tumor <NA> # MGH104-P5-A06 tumor <NA> # MGH104-P5-A07 tumor <NA> # MGH104-P5-A08 tumor <NA> # 包含了三个样本以及他们共享的位点 names (multi_result[[ 2 ]]) # [1] "all_MGH104" "all_MGH105" "all_MGH102" # [4] "all_shared" "all_shareSubshared" # 包含了染色体编号，染色体起点和终点以及拷贝数变异数值 head (multi_result[[ 2 ]][[ 1 ]]) # Chr Start End Alteration segm.mean # 5 3 47228026 49535618 2 0.174473 # 7 3 61561569 111665750 -2 -0.151862 # 11 4 152618632 168928457 -2 -0.225926 # 24 11 62832342 66958376 2 0.151061 # 27 12 53251251 57095476 2 0.262316 # 31 14 66507407 75282230 2 0.155383 # 所有样本都共有的变异位点 head (multi_result[[ 2 ]][[ 4 ]]) # Chr Start End Alteration segm.mean # 13 5 92151 147782848 2 0.040585 # 16 7 192969 159144957 2 0.088789 # 21 10 134465 133421307 -2 -0.114114 # 43 22 29767369 50801309 -2 -0.159207 # 部分样本共有的变异位点 head (multi_result[[ 2 ]][[ 5 ]]) # Chr Start End Alteration segm.mean sh_sub Mean # all_share2.6 2 1646536241817413 2 0.069152 2-3 0 # all_share2.11 5 148825245 181261139 -1 -0.098875 2-3 0 # all_share1.14 6 291630 28093664 -2 -0.099298 1-3 0 # all_share1.19 9 89318506 122913330 -2 -0.219364 1-3 0 # all_share2.23 11 78068304 134411918 1 0.098385 2-3 0 # all_share2.40 19 43646095 55442863 -2 -0.187000 2-3 0 三个样本拷贝数变异对比热图，文件路径为：/data/09_scRNA-seq_cnv/result_SCEVAN/multi_sample/all_compareClonalCN.png 三个样本之间的共享改变或不同克隆改变。

这些片段分为五种预定义的拷贝数状态之一：缺失、丢失、中性、增益或扩增。

文件路径为：/data/09_scRNA-seq_cnv/img/allOncoHeat2.png 6、联合Seurat对象 6.1 创建Seurat对象 library (Seurat) seurObj <- Seurat :: CreateSeuratObject (count_mtx, meta.data = single_result) # 或者可以使用下面方法增加metadata seurObj <- Seurat :: AddMetaData (seurObj, metadata = single_result) seurObj # An object of class Seurat # 23686 features across 200 samples within 1 assay # Active assay: RNA (23686 features, 0 variable features) # 2 layers present: counts, data 6.2 将拷贝数变异结果加入Seurat对象 # 加载前面的分析结果 # 基因注释结果 load ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106_count_mtx_annot.RData" ) # 拷贝数变异矩阵 load ( "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106_CNAmtx.RData" ) # 分别对应了基因名称、染色体编号、染色体对应起点和重点和symbol名称 count_mtx_annot[ 1 : 5 , 1 : 5 ] # seqnames start end gene_id gene_name # ENSG00000227232 1 14404 29570 ENSG00000227232 WASH7P # ENSG00000225880 1 826206 827522 ENSG00000225880 LINC00115 # ENSG00000230368 1 868071 876903 ENSG00000230368 FAM41C # ENSG00000188976 1 944204 959309 ENSG00000188976 NOC2L # ENSG00000187608 1 1001138 1014541 ENSG00000187608 ISG15 # 查看数据维度 dim (count_mtx_annot) # [1] 10566 5 # 列分别对应了细胞名称 CNA_mtx_relat[ 1 : 5 , 1 : 5 ] # MGH106-P7-E10 MGH106-P2-B01 MGH106-P2-B12 MGH106-P2-H07 MGH106-P2-E09 # [1,] 0.197378 -0.09338299 0.01316769 -0.07037455 -0.08075804 # [2,] 0.197378 -0.09338299 0.01316769 -0.07037455 -0.08075804 # [3,] 0.197378 -0.09338299 0.01316769 -0.07037455 -0.08075804 # [4,] 0.197378 -0.09338299 0.01316769 -0.07037455 -0.08075804 # [5,] 0.197378 -0.09338299 0.01316769 -0.07037455 -0.08075804 # 查看数据维度 dim (CNA_mtx_relat) # [1] 10566 200 # 从CNA_mtx_relat中筛选出染色体3上特定区域的数据，计算每列的均值，然后按照[email protected]的行名对结果进行重新排序和命名，最后将其转换为数据框 chr3 <- apply (CNA_mtx_relat[count_mtx_annot $ seqnames == 3 & count_mtx_annot $ start >= 158644278 & count_mtx_annot $ end <= 194498364 ,], 2 , mean) head (chr3) # MGH106-P7-E10 MGH106-P2-B01 MGH106-P2-B12 MGH106-P2-H07 MGH106-P2-E09 # -0.073484354 0.143933614 0.006989981 -0.029519201 -0.074089242 # MGH106-P2-E04 # 0.022922624 chr3 <- chr3[ rownames (seurObj @ meta.data)] names (chr3) <- rownames (seurObj @ meta.data) chr3 <- as.data.frame (chr3) head (chr3) # chr3 # MGH106-P7-E10 -0.073484354 # MGH106-P2-B01 0.143933614 # MGH106-P2-B12 0.006989981 # MGH106-P2-H07 -0.029519201 # MGH106-P2-E09 -0.074089242 # MGH106-P2-E04 0.022922624 # 将chr3结果加入Seurat对象中 seurObj <- AddMetaData (seurObj, metadata = chr3) 6.3 降维分群 seurObj <- NormalizeData (seurObj , normalization.method = "LogNormalize" , scale.factor = 10000 ) # 寻找高变基因 seurObj <- FindVariableFeatures (seurObj , selection.method = "vst" , nfeatures = 2000 ) # 归一化数据 seurObj <- ScaleData (seurObj , features = VariableFeatures (seurObj)) # 计算线粒体含量 seurObj [[ "percent.mt" ]] <- PercentageFeatureSet (seurObj , pattern = "^MT-" ) seurObj <- ScaleData (seurObj , vars.to.regress = "percent.mt" ) [["RNA"]]@scale.data seurObj <- RunPCA (seurObj , features = VariableFeatures ( object = seurObj )) ElbowPlot (seurObj ) # 寻找邻近值 seurObj <- FindNeighbors (seurObj , dims = 1 : 15 ) seurObj <- RunUMAP (seurObj , dims = 1 : 15 ) 6.4 可视化 # 展示正常/肿瘤的降维图 DimPlot (seurObj , reduction = "umap" , group.by = "class" , label = TRUE ) # 展示chr3染色体选取区域拷贝数变异VlnPlot # 将[email protected]$subclone的NA对应的正常细胞名称调整为normal seurObj @ meta.data $ subclone[ is.na (seurObj @ meta.data $ subclone)] <- "normal" seurObj @ meta.data $ subclone[seurObj @ meta.data $ subclone == 1 ] <- "clone1" seurObj @ meta.data $ subclone[seurObj @ meta.data $ subclone == 2 ] <- "clone2" seurObj @ meta.data $ subclone[seurObj @ meta.data $ subclone == 3 ] <- "clone3" table (seurObj @ meta.data $ subclone) # # clone1 clone2 clone3 normal # 58 57 21 64 VlnPlot (seurObj, "chr3" , group.by = "subclone" ) # 热图 plotCNA_withAnnotCells ( SampleName = "MGH106" , # 分析时用的样本名称 outputPATH= "/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/" , # 输出结果目录，和前面保持一致 metadata = seurObj @ meta.data, # 包含分类信息的列表 COLUMNS_TO_PLOT = c ( "class" , "subclone" )) # 需要展示分类变量名称 # [1] "class" # [1] "subclone" # png # 2 热图的文件路径为：/data/09_scRNA-seq_cnv/result_SCEVAN/single_sample/MGH106heatmap_with_annotation.png Sys.time # [1] "2025-08-29 18:38:18 CST" save.image ( "result_SCEVAN/SCEVAN.rdata" )

七、numbat 本文介绍了一种创新的计算方法Numbat，通过整合基于群体的单倍型信息、等位基因信号和表达信号，显著提升了CNV检测的准确性。Numbat利用肿瘤亚克隆间的进化关系，迭代推断单细胞拷贝数谱及其克隆系统发育树，成功应用于多发性骨髓瘤、乳腺癌和甲状腺癌等多种肿瘤类型。本文旨在阐述Numbat的工作原理及其在揭示肿瘤微环境和治疗耐药机制中的应用潜力。

1、简介 DOI: 10.1038/s41587-022-01468-y 链接：https://pubmed.ncbi.nlm.nih.gov/36163550/ 1.1 利用单倍型信息实现灵敏的CNV检测 Numbat算法的核心创新是：整合单倍型（haplotype）信息来显著提升从稀疏的scRNA-seq数据中检测拷贝数变异（CNV）的灵敏度。

作者首先验证了尽管scRNA-seq数据中的SNP覆盖度很低，但利用群体水平的单倍型推断（population-based haplotype phasing）工具（如Eagle2）依然能够有效地推断出等位基因在染色体上的连锁关系，甚至能够连接不同基因座上的SNP。这一跨基因的定相能力至关重要。

因为它使得算法能够区分真实的、影响一整段染色体的CNV信号，与单个基因随机发生的等位基因特异性表达（Allele-Specific Expression,ASE）所带来的噪声。基于这一原理，Numbat构建了一个单倍型感知的隐马尔可夫模型（haplotype-aware HMM），该模型通过分析定相后单倍型频率的定向偏移，而非仅仅观察等位基因频率的方差，从而获得了更强的统计效力（图1a）。

在对“伪群体”肿瘤-正常细胞混合样本的基准测试中，该模型成功识别出了在低肿瘤细胞比例下传统方法无法检出的微弱CNV信号（图1b,c），并且在单细胞层面实现了更准确的等位基因分配，使得单个肿瘤细胞的等位基因不平衡信号更加清晰可见（图1f,g）。

1.2 从转录组推断等位基因特异性拷贝数 为了获得更稳健、更全面的拷贝数变异（CNV）图谱，Numbat不仅利用单倍型信息，还创新性地将等位基因不平衡（allelic imbalance）和基因表达丰度（expression magnitude）这两种正交的信号整合到一个统一的分析框架中。

其核心是一个基于生成式统计框架的联合隐马尔可夫模型（joint HMM），该模型能够同时对两种信号进行建模，从而推断出等位基因特异性的拷贝数状态（例如，父源染色体扩增、母源染色体缺失等），而不仅仅是总拷贝数的变化。

此外，为解决其他工具常因肿瘤整体倍性变化而导致基线判断错误的问题，Numbat采用了一种巧妙的两步法来识别真实的二倍体基线（diploid baseline），即先通过等位基因信息找到平衡区域，再从中选择表达水平最低的区域作为参照。

在与拥有“金标准”全基因组测序（WGS）数据的五个多发性骨髓瘤样本的比较中，Numbat推断出的CNV图谱与真实的DNA图谱高度一致，其准确率（精确率99.2%，召回率95.4%）显著优于HoneyBADGER、InferCNV和CopyKAT等其他主流方法（图2d）。特别是，Numbat能够准确识别出拷贝数中性杂合性丢失（copy-neutral LoH）等仅依靠表达量无法发现的事件（图2c），展现了其整合多信息的强大能力。

1.3 推断肿瘤克隆架构与演化历史 Numbat算法解决了肿瘤异质性分析中的一个核心难题：拷贝数变异（CNV）的检测与肿瘤克隆结构的推断是相互依赖的。在解决这两个问题的过程中，Numbat采用了一种创新的迭代优化程序（alternating optimization procedure），以联合推断单细胞CNV图谱和相关的亚克隆系统发育树（subclonal phylogeny）。这个过程始于一个基于表达谱聚类的初始演化树（图3a）。

在每一次迭代中，算法首先会根据当前的演化树结构，将细胞聚合成更精确的、谱系特异性的“伪群体（pseudobulks）”，并在此基础上运行其联合HMM模型以识别共享的CNV事件（图3b）。随后，算法会利用这些新发现的CNV事件，反过来为每一个单细胞计算其携带该事件的后验概率。

这个包含不确定性信息的概率矩阵将被输入到一个最大似然完美系统发育（maximum-likelihood perfect phylogeny）的框架中（ScisTree），以构建一个更准确、更精细的演化树。这个新生成的演化树将作为下一次迭代的起点，整个过程循环往复，直至克隆结构和CNV图谱收敛稳定。通过这种相互优化的迭代策略，Numbat能够同时精准地解析出肿瘤的克隆架构和每个细胞的CNV状态（图3c）。

1.4 肿瘤与正常细胞的可靠分类 作者在一个包含18个肿瘤样本（涵盖三阴性乳腺癌、甲状腺癌和多发性骨髓瘤）的数据集上，将Numbat的分类准确性与CopyKAT进行了比较。在乳腺癌和甲状腺癌这两个实体瘤样本中，Numbat的准确率极高（平均约98.5%），与CopyKAT的表现相当。

然而，在多发性骨髓瘤样本中，两者表现出了显著差异：Numbat依然保持了稳定的高准确率（98.7%），而CopyKAT在八个样本中的五个都出现了细胞簇的错误分类，导致其平均准确率骤降至74.7%。作者推测，这可能是因为多发性骨髓瘤样本的单细胞测序深度较低，且其染色体畸变程度不如实体瘤那般剧烈。

这一结果凸显了Numbat的稳健性，因为它通过整合基因表达和等位基因这两种正交的证据来源，增强了信号的可靠性，从而在更具挑战性的数据中也能做出准确的判断。1.5 单倍型感知CNV分析揭示亚克隆复杂性 作者向读者展示了Numbat算法最核心的优势：利用其独特的单倍型感知（Haplotype-aware）能力，揭示了传统scRNA-seq分析方法无法识别的、复杂的肿瘤亚克隆结构。

作者首先通过与scDNA-seq数据的直接比较，证实了Numbat能够仅从scRNA-seq数据中就准确地重构出与DNA层面一致的亚克隆结构（Extended Data Fig.9）。在对一个三阴性乳腺癌样本（TNBC1）的分析中，Numbat揭示了一个复杂的、具有两个主要演化分支的亚克隆谱系。

重要的是它识别出了区分这些亚克隆的、极为精细的等位基因特异性事件，包括亚克隆级别的拷贝数中性杂合性丢失（copy-neutral LoH, CNLoH），以及此前从未在scRNA-seq数据中被报道过的镜像的单倍型特异性扩增（mirrored haplotype-specific amplifications）——即两个不同的亚克隆虽然都扩增了同一条染色体，但它们扩增的分别是来自不同亲本（父源或母源）的染色体拷贝（图4a-e）。

在另一个甲状腺癌样本（ATC1）中，Numbat同样发现了具有相互排斥畸变（reciprocal aberrations）的两个亚克隆，一个克隆扩增7号染色体并丢失17号，另一个则反之（图4f-i）。这些复杂的、仅在等位基因层面可见的演化事件，凸显了整合单倍型信息对于深入理解肿瘤异质性的巨大威力。

1.6 遗传异质性与转录异质性之间的相互作用 研究团队展示了Numbat的终极应用：通过整合分析，深入探究肿瘤的遗传异质性（genetic heterogeneity）与转录异质性（transcriptional heterogeneity）之间的复杂相互作用，并追踪其在治疗过程中的演化。作者以一个多发性骨髓瘤（multiple myeloma）患者的连续样本（包含初诊、缓解期和两次复发期）为案例进行了深入分析。

Numbat首先根据拷贝数变异（CNV）识别出了三个遗传亚克隆（three genetic subclones），并清晰地描绘了它们在治疗压力下的演化路径：一个祖先克隆（g1）在初次治疗后幸存下来，并在复发过程中演化出一个新的、更具侵袭性的克隆（g3），最终在第二次复发时成为优势克隆（图5c）。

通过将遗传克隆信息与细胞的转录状态相结合，团队成员发现，在初诊样本中，这些遗传克隆分布在两种不同的转录状态中，并据此推断出肿瘤演化的先后顺序：一次大规模的转录程序转变先于一次亚克隆CNV事件的获得（图5d）。通过比较不同遗传亚克隆的基因表达谱，Numbat成功地将特定的CNV事件与其在反式（trans）作用下调控的关键通路联系起来。

例如，最终耐药克隆（g3）所获得的16q染色体缺失，与干扰素γ响应通路（interferon-gamma response pathway）的下调显著相关，这是一种已知的肿瘤免疫逃逸机制（图5i）。这一系列分析充分展示了Numbat的强大能力，它能够将肿瘤的基因型、转录表型与临床过程（如耐药）联系起来，从scRNA-seq数据中挖掘出深刻的生物学机制。

1.7 小结 Numbat通过整合单倍型信息和scRNA-seq数据，显著提升了肿瘤克隆结构的解析能力，成功区分恶性与非恶性细胞，并揭示了与肿瘤进展和治疗耐药相关的遗传和转录亚群。这项研究在22个肿瘤样本中的应用验证了其在多种癌症类型中的鲁棒性，尤其是在识别单倍型特异性CNV和克隆复杂性方面展现出独特优势。不过Numbat在处理复杂基因组变异（如全基因组倍增）时仍需手动校正，提示未来需改进基线倍性估计方法。

展望未来，结合多组学数据（如表观遗传学信息）将进一步增强Numbat的功能，为解析基因组不稳定性对肿瘤细胞状态的影响提供更全面视角。Numbat的开源特性及其广泛适用性，使其有望成为肿瘤异质性研究和精准医疗的重要工具。

2、安装R包 install.packages ( 'memuse' , dependencies = TRUE ) install.packages ( 'vcfR' , dependencies = TRUE ) install.packages ( 'numbat' , dependencies = TRUE ) install.packages ( 'optparse' , dependencies = TRUE ) 3、加载R包 # 清空环境变量 rm ( list = ls ) gc # used (Mb) gc trigger (Mb) max used (Mb) # Ncells 11072010 591.4 25853280 1380.8 25853280 1380.8 # Vcells 21354993 163.0 142928480 1090.5 1006972481 7682.6 # 展示系统时间 Sys.time # [1] "2025-08-29 18:38:30 CST" # 设置工作目录 setwd ( "/data/09_scRNA-seq_cnv/" ) library (numbat) library (Seurat) # 如果不清楚自己的numbat包的安装路径，可以使用下面命令查找 find.package ( "numbat" ) # [1] "/usr/local/lib/R/site-library/numbat" 4、数据下载 numbat需要一个包含了SNP位点的VCF文件和1000G Reference Panel。

这里使用来自千人基因组计划（1000 Genomes Project）第三阶段（Phase 3）的常见单核苷酸多态性（SNP）数据，覆盖染色体1号到X号（不包括Y染色体和线粒体基因组）。

以及1000G Reference Panel（1000 Genomes Reference Panel）这是用于 基因型填补（genotype imputation） 的 参考数据集，包括了来自1000 Genomes Project的高质量基因型数据，包括不同人群的数百万个SNP。第三个文为人类遗传图谱文件。注意这里基因组包含了两个版本，根据自己的实际分析参考基因组进行选择。

# hg38 http : // pklab.med.harvard.edu / teng / data / 1000 G_hg38.zip https : // sourceforge.net / projects / cellsnp / files / SNPlist / genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz https : // storage.googleapis.com / broad - alkesgroup - public / Eagle / downloads / tables / genetic_map_hg38_withX.txt.gz # hg19 http : // pklab.med.harvard.edu / teng / data / 1000 G_hg19.zip https : // sourceforge.net / projects / cellsnp / files / SNPlist / genome1K.phase3.SNP_AF5e2.chr1toX.hg19.vcf.gz https : // storage.googleapis.com / broad - alkesgroup - public / Eagle / downloads / tables / genetic_map_hg19_withX.txt.gz 5、等位基因数据准备 5.1 运行参数解释：label : 个体标签，每次运行一个 samples : 样本名称，若多个则以逗号分隔 bams : BAM 文件，每个样本一个，若多个则以逗号分隔 barcodes : 细胞barcode文件，每个样本一个，若多个则以逗号分隔 gmap : 遗传图谱路径（例如，Eagle_v2.4.1/tables/genetic_map_hg38_withX.txt.gz） eagle : Eagle2 二进制文件路径 snpvcf : 包含了SNP信息的vcf文件 paneldir : reference panel文件夹目录 (bcf文件) outdir : 输出文件路径 ncores : 运行核心数 smartseq : 以SMART-seq模式运行 5.2 示例脚本 # 修改文件权限，注意这里需要切换至linux终端运行 sudo chmod 777 / usr / local / lib / R / site - library / numbat / bin / pileup_and_phase.R # 执行分析，注意这里同样需要在linux终端运行 Rscript / numbat / inst / bin / pileup_and_phase.R \ -- label {sample} \ -- samples {sample} \ -- bams / mnt / mydata / {sample}.bam \ -- barcodes / mnt / mydata / {sample}_barcodes.tsv \ -- outdir / mnt / mydata / {sample} \ -- gmap / Eagle_v2. 4.1 / tables / genetic_map_hg38_withX.txt.gz \ -- snpvcf / data / genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf \ -- paneldir / data / 1000 G_hg38 \ -- ncores 5.3 可能出现的报错 5.3.1 optparse问题 # 如果报错:Error in library(optparse) : there is no package called ‘optparse’ # find.package("optparse")又确认已经安装完成，执行下面脚本，查看你的当前调用R包路径位置 Rscript - e '.libPaths' # 此时你可能位于某个conda环境，可以执行 conda deactivate 5.3.2 cellsnp-lite问题 cellsnp-lite 是一个专为单细胞测序数据设计的高性能SNP calling工具，主要功能包括：SNP pileup / 基因型计数（最主要功能）：从BAM文件中提取参考碱基和变异碱基的计数信息（allele counts），适用于单细胞数据 支持多种模式 mode 1a: 使用已知 SNP（VCF），在基于每个细胞分析 mode 1b: 使用已知 SNP（VCF），对所有细胞整体（不分 barcode）分析 mode 2a: droplet-based的单细胞数据，没有给定SNPs mode 2b: well-based单细胞数据或者bulk数据，没有给定SNPs 支持并行计算（多线程） 使用 -p 参数可以启用多线程加速处理大数据集 支持多种标签 自动识别或手动指定 cell barcode (–cellTAG) 和 UMI (–UMItag)，常见设置为：–cellTAG CB 和 –UMItag UB # 如果报错：/data/09_scRNA-seq_cnv/data/result_numbat/pbmc3k/run_pileup.sh: 1: cellsnp-lite: Permission denied，说明缺少了cellsnp-lite这个软件，接下来需要先安装这个软件，执行命令：conda create - n cellsnp - lite - c bioconda cellsnp - lite # 激活环境 conda activate cellsnp - lite # 测试是否安装成功 cellsnp - lite -- help 5.3.3 eagle问题 # 如果报错：/data/09_scRNA-seq_cnv/result_numbat/pbmc3k/run_phasing.sh: 22: eagle: Permission denied # 激活环境 conda activate cellsnp - lite conda install - c bioconda eagle -- help # 但是，还是有问题，这是因为当前我的conda只有V1版本的eagle，但是其实我需要V2版本的 cd / data / 09 _scRNA - seq_cnv / data wget https : // data.broadinstitute.org / alkesgroup / Eagle / downloads / Eagle_v2. 4 . 1. tar.gz # 解压缩 tar - xzvf Eagle_v2. 4 . 1. tar.gz cd Eagle_v2. 4.1 # 检测是否正常使用 / data / 09 _scRNA - seq_cnv / data / Eagle_v2. 4.1 / eagle -- help # 为了便于后续继续使用当前最新版本，将软件地址添加到环境变量中去 export PATH = / data / 09 _scRNA - seq_cnv / data / Eagle_v2. 4.1 : $ PATH source ~ / .bashrc conda activate cellsnp - lite eagle -- help 5.4 执行分析 # 修改文件权限，注意这里需要切换至linux终端运行 sudo chmod 777 /usr/local/lib/R/site-library/numbat/bin/pileup_and_phase.R # 执行分析，注意这里同样需要在linux终端运行 Rscript /usr/local/lib/R/site-library/numbat/bin/pileup_and_phase.R \ --label pbmc3k \ --samples pbmc3k \ --bams /data/09_scRNA-seq_cnv/data/pbmc3k_possorted_genome_bam.bam \ --barcodes /data/09_scRNA-seq_cnv/data/hg19/barcodes.tsv \ --outdir /data/09_scRNA-seq_cnv/result_numbat/pbmc3k_new \ --gmap /data/09_scRNA-seq_cnv/data/genetic_map_hg19_withX.txt.gz \ --snpvcf /data/09_scRNA-seq_cnv/data/genome1K.phase3.SNP_AF5e2.chr1toX.hg19.vcf \ --paneldir /data/09_scRNA-seq_cnv/data/1000G_hg19 \ --ncores 4 6、运行numbat 6.1 运行参数解释 CNV 检测 t：HMM 中使用的转移概率。

较低的t更适合于拷贝数景观更复杂的肿瘤 (从中可以预期更多的断点), 有时更适合检测亚克隆事件。较高的t对于控制CNV调用的假阳性率更有效。gamma：等位基因计数的过度离散 (默认值：20)。对于10X的数据，建议使用20。非UMI协议 (例如SMART-Seq) 通常会产生更嘈杂的等位基因数据，建议使用较小的gamma值（例如5）。min_cells：运行pseudobulk HMM的最小单元数。

如果对于一个数据集来说，每个单元的等位基因覆盖率非常稀疏，可以考虑将这个阈值设置得更高。multi_allelic：是否启用多等位基因CNV 筛选 min_LLR：最小对数似然比阈值用于过滤 CNV (默认值：5)。为了确保系统发育推理的质量，我们只使用可靠的CNV来重建系统发育。max_entropy：在系统发育学构建之前用于过滤CNV的另一个标准 (默认值：0.5)。二进制后验的熵量化了单个细胞中事件的不确定性。

该值可以是0到1之间的任意值，其中1是最不严格的。系统发育 n_cut：系统发育中用于定义亚克隆的切割次数。请注意，n次切割应导致n+1个克隆 (顶级正常二倍体克隆始终包括在内)。或者，可以指定max_cost或tau来控制亚克隆的数量。max_cost：内部分支崩溃的可能性阈值。较高的max_cost通常会导致较少的克隆。tau：简化突变历史的严格性。基本上将max_cost设置为tau乘以细胞数。

迭代优化 init_k：用于hclust初始化的初始子簇数。默认情况下，Numbat使用平滑表达式值的层次聚类(hclust)来近似初始系统发育谱系。这将将初始树切割成k个簇。更多的簇意味着在亚克隆CNV检测的初始阶段有更高的分辨率。默认情况下，我们将init_k设置为3。max_iter：最大迭代次数。Numbat 迭代优化系统发育和拷贝数估计。在实践中，我们发现 2 次迭代后的结果通常是稳定的。

check_convergence：如果结果已收敛，请停止迭代 (基于共识 CNV 段落)。并行 ncores：用于单细胞CNV测试的核心数量 ncores_nni：用于系统发育推理的核心数量 其他运行模式 call_clonal_loh：是否调用克隆 LOH 区域。建议用于没有匹配正常细胞的高纯度样品。有关详细信息，请参阅 检测克隆LOH。

segs_consensus_fix：一个预先确定的共识CNV数据框架 (例如，来自批量WGS/WES分析), 在分析过程中保持固定。有关详细信息，请参阅 使用现有的CNV调用。6.2 运行numbat 示例数据来源于 Gao等人的ATC2样本 , 基因表达计数矩阵和等位基因数据已经准备好，可以直接读取。

# 在线读取可能会失败 # count_mat_ATC2 = readRDS(url('http://pklab.med.harvard.edu/teng/data/count_mat_ATC2.rds')) # df_allele_ATC2 = readRDS(url('http://pklab.med.harvard.edu/teng/data/df_allele_ATC2.rds')) # 直接从本地读取文件 count_mat_ATC2 <- readRDS ( "data/count_mat_ATC2.rds" ) df_allele_ATC2 <- readRDS ( "data/df_allele_ATC2.rds" ) # 运行 out = run_numbat ( count_mat_ATC2, # 基因 x 细胞的UMI矩阵 ref_hca, # 参考表达谱，一个基因x细胞类型归一化表达水平矩阵 df_allele_ATC2, # 由pileup_and_phase脚本生成的等位基因数据 genome = "hg38" , # 参考基因组 t = 1e-5 , # HMM 中使用的转移概率 ncores = 4 , # 分析核心数 plot = TRUE , # 是否绘制图形 out_dir = './result_numbat/' ) 7、结果解读 7.1 拷贝数景观和单细胞系统发育 在这个可视化中，单细胞系统发育(左)与单细胞CNV调用的热图(右)并列。

CNV调用按照改变类型(AMP：扩增，BAMP：平衡扩增，DEL：缺失，CNLoH：拷贝中性杂合性缺失)进行着色。中间的颜色条区分了不同的遗传群体(基因型)。虚线蓝线将预测的肿瘤与正常细胞分开。这告诉我们，数据集主要由三个细胞群体、一个正常群体(灰色)和两个肿瘤亚克隆(绿色和紫色)组成。

library (ggplot2) library (numbat) library (dplyr) library (glue) library (data.table) library (ggtree) library (stringr) library (tidygraph) library (patchwork) nb = Numbat $ new ( out_dir = './result_numbat/' ) # Numbat概率性地评估单个细胞中CNV的存在/不存在。

单个细胞级别的CNV后验值存储在nb$joint_post中。其中包含特定CNV片段 (seg)、其改变类型(cnv_state)、CNV的联合后验概率(p_cnv)、基于表达的后验(p_cnv_x)和基于等位基因的后验 (p_cnv_y) 的细胞级信息。# p_cnv_x：expression-based posterior：是基于表达数据（scRNA-seq）计算出的后验概率，表示该基因组片段属于某种CNV状态的置信度。

数据基因表达矩阵（例如表达下降暗示缺失，表达升高可能提示扩增） # p_cnv_y：allele-based posterior是基于等位基因信息（SNP pileup）计算出的后验概率，通过等位基因表达失衡推断 CNV。

数据由cellsnp-lite生成的 allele counts head (nb $ joint_post) %>% select (cell, CHROM, seg, cnv_state, p_cnv, p_cnv_x, p_cnv_y) # cell CHROM seg cnv_state p_cnv p_cnv_x # <char> <fctr> <char> <char> <num> <num> # 1: ATC2_AACCACAGTCCATCTC-1 1 1a del 1.000000e+00 1.000000e+00 # 2: ATC2_AACCACAGTCCATCTC-1 2 2a del 1.000000e+00 2.580679e-15 # 3: ATC2_AACCACAGTCCATCTC-1 3 3a del 1.000000e+00 9.968479e-01 # 4: ATC2_AACCACAGTCCATCTC-1 4 4a del 1.000000e+00 9.999929e-01 # 5: ATC2_AACCACAGTCCATCTC-1 5 5b del 4.998742e-01 4.998742e-01 # 6: ATC2_AACCACAGTCCATCTC-1 5 5c del 2.332953e-27 7.654833e-37 # p_cnv_y # <num> # 1: 1.000000e+00 # 2: 1.000000e+00 # 3: 1.000000e+00 # 4: 1.000000e+00 # 5: 5.000000e-01 # 6: 5.188836e-24 mypal = c ( '1' = ' , '2' = " , '3' = " , '4' = " , '5' = " ) nb $ plot_phylo_heatmap ( clone_bar = TRUE , p_min = 0.9 , pal_clone = mypal ) 7.2 在系统发育上完善亚克隆 在cutree中可以指定n_cut或max_cost。

注意，n次切割应该产生n+1个克隆（顶级正常二倍体克隆总是包含在内）。第一次切割应分离出肿瘤细胞作为单个克隆，第二次切割产生两个肿瘤亚克隆，依此类推。或者，可以指定max_cost，这是减少系统发育的最大可能性成本阈值。

plots = lapply ( 1 : 4 , function (k) { nb $ cutree ( n_cut = k) nb $ plot_phylo_heatmap + ggtitle ( paste0 ( 'n_cut=' , k)) } ) wrap_plots (plots) 7.3 共识拷贝数区段 指在多个细胞、样本或分析方法中共同识别出的一致性拷贝数变化区域 nb $ plot_consensus 7.4 pseudobulks水平拷贝数变异 我们还可以在数据更丰富的pseudobulks中可视化这些CNV事件，通过克隆对单元进行聚合。

nb $ bulk_clones %>% filter (n_cells > 50 ) %>% # 只保留那些 细胞数大于50的克隆群体，避免绘制掉小而可能不可靠的克隆 plot_bulks ( min_LLR = 10 , # 只显示那些具有较强统计支持（对数似然比 > 10）的CNV legend = TRUE ) 7.5 进化树 nb $ plot_sc_tree ( label_size = 3 , branch_width = 0.5 , tip_length = 0.5 , pal_clone = mypal, tip = TRUE ) Sys.time # [1] "2025-08-29 18:55:39 CST" save.image ( "result_numbat/numbat.rdata" ) 环境版本信息 sessionInfo # R version 4.4.2 (2024-10-31) # Platform: x86_64-pc-linux-gnu # Running under: Ubuntu 20.04.6 LTS # # Matrix products: default # BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 # LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0 # # locale: # [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C # [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 # [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 # [7] LC_PAPER=en_US.UTF-8 LC_NAME=C # [9] LC_ADDRESS=C LC_TELEPHONE=C # [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C # # time zone: Asia/Shanghai # tzcode source: system (glibc) # # attached base packages: # [1] stats4 parallel grid stats graphics grDevices utils # [8] datasets methods base # # other attached packages: # [1] patchwork_1.3.0 tidygraph_1.2.2 # [3] data.table_1.16.4 glue_1.8.0 # [5] numbat_1.4.2 Matrix_1.7-1 # [7] fgsea_1.24.0 ggrepel_0.9.6 # [9] ggtree_3.6.2 ape_5.8-1 # [11] tidytree_0.4.1 Rtsne_0.17 # [13] igraph_2.0.1.1 scran_1.26.0 # [15] scuttle_1.8.0 SingleCellExperiment_1.20.0 # [17] SummarizedExperiment_1.28.0 Biobase_2.66.0 # [19] GenomicRanges_1.50.1 GenomeInfoDb_1.42.3 # [21] IRanges_2.32.0 S4Vectors_0.36.0 # [23] BiocGenerics_0.52.0 MatrixGenerics_1.18.1 # [25] matrixStats_1.5.0 yaGST_2017.08.25 # [27] doParallel_1.0.17 iterators_1.0.14 # [29] foreach_1.5.2 SCEVAN_1.0.3 # [31] copykat_1.1.0 devtools_2.4.5 # [33] usethis_2.1.6 ggsignif_0.6.4 # [35] AnnoProbe_0.1.7 RColorBrewer_1.1-3 # [37] circlize_0.4.16 ComplexHeatmap_2.14.0 # [39] forcats_1.0.0 stringr_1.5.1 # [41] dplyr_1.1.4 purrr_1.0.2 # [43] readr_2.1.5 tidyr_1.3.1 # [45] tibble_3.2.1 ggplot2_3.5.1 # [47] tidyverse_1.3.2 Seurat_5.0.3 # [49] SeuratObject_5.0.1 sp_2.1-4 # [51] infercnv_1.14.2 # # loaded via a namespace (and not attached): # [1] ica_1.0-3 plotly_4.10.1 # [3] zlibbioc_1.52.0 tidyselect_1.2.1 # [5] rvest_1.0.3 clue_0.3-66 # [7] lattice_0.22-6 rjson_0.2.23 # [9] urlchecker_1.0.1 png_0.1-8 # [11] cli_3.6.3 ggplotify_0.1.0 # [13] usdata_0.3.1 goftest_1.2-3 # [15] gargle_1.2.1 textshaping_1.0.1 # [17] kernlab_0.9-33 bluster_1.8.0 # [19] BiocNeighbors_2.0.1 uwot_0.1.14 # [21] dendextend_1.19.0 mime_0.12 # [23] evaluate_1.0.3 leiden_0.4.3 # [25] coin_1.4-2 openintro_2.4.0 # [27] stringi_1.8.4 backports_1.5.0 # [29] rjags_4-16 parallelDist_0.2.6 # [31] lubridate_1.9.0 httpuv_1.6.6 # [33] magrittr_2.0.3 splines_4.4.2 # [35] ggraph_2.1.0 logger_0.4.0 # [37] DT_0.26 sctransform_0.4.1 # [39] ggbeeswarm_0.7.2 sessioninfo_1.2.2 # [41] DBI_1.2.3 jquerylib_0.1.4 # [43] withr_3.0.2 systemfonts_1.1.0 # [45] lmtest_0.9-40 formatR_1.14 # [47] htmlwidgets_1.5.4 fs_1.6.5 # [49] segmented_2.1-3 cherryblossom_0.1.0 # [51] labeling_0.4.3 cellranger_1.1.0 # [53] mixtools_1.2.0 reticulate_1.26 # [55] zoo_1.8-12 XVector_0.38.0 # [57] knitr_1.49 hahmmr_1.0.0 # [59] UCSC.utils_1.2.0 airports_0.1.0 # [61] RhpcBLASctl_0.23-42 timechange_0.3.0 # [63] caTools_1.18.3 R.oo_1.27.0 # [65] quantreg_5.94 RSpectra_0.16-2 # [67] irlba_2.3.5.1 ggrastr_1.0.1 # [69] fastDummies_1.7.4 gridGraphics_0.5-1 # [71] ellipsis_0.3.2 lazyeval_0.2.2 # [73] yaml_2.3.10 phyclust_0.1-34 # [75] survival_3.4-0 scattermore_1.2 # [77] crayon_1.5.3 RcppAnnoy_0.0.22 # [79] progressr_0.15.1 tweenr_2.0.3 # [81] scistreer_1.2.0 later_1.4.1 # [83] ggridges_0.5.6 codetools_0.2-20 # [85] GlobalOptions_0.1.2 profvis_0.3.7 # [87] shape_1.4.6.1 limma_3.62.2 # [89] pkgconfig_2.0.3 xml2_1.3.6 # [91] ggpubr_0.5.0 aplot_0.1.8 # [93] spatstat.sparse_3.1-0 viridisLite_0.4.2 # [95] xtable_1.8-4 car_3.1-1 # [97] fastcluster_1.2.6 plyr_1.8.9 # [99] httr_1.4.4 tools_4.4.2 # [101] globals_0.16.3 pkgbuild_1.3.1 # [103] beeswarm_0.4.0 broom_1.0.7 # [105] nlme_3.1-160 futile.logger_1.4.3 # [107] lambda.r_1.2.4 dbplyr_2.5.0 # [109] MatrixModels_0.5-3 digest_0.6.37 # [111] farver_2.1.2 tzdb_0.4.0 # [113] reshape2_1.4.4 yulab.utils_0.1.9 # [115] viridis_0.6.5 cachem_1.1.0 # [117] polyclip_1.10-7 generics_0.1.3 # [119] mvtnorm_1.3-3 googledrive_2.0.0 # [121] presto_1.0.0 parallelly_1.41.0 # [123] pkgload_1.4.0 statmod_1.5.0 # [125] RcppHNSW_0.6.0 ragg_1.4.0 # [127] ScaledMatrix_1.6.0 carData_3.0-5 # [129] pbapply_1.7-2 mcmc_0.9-8 # [131] MCMCpack_1.7-1 spam_2.11-0 # [133] dqrng_0.4.1 graphlayouts_0.8.3 # [135] gtools_3.9.5 readxl_1.4.3 # [137] gridExtra_2.3 shiny_1.10.0 # [139] GenomeInfoDbData_1.2.13 R.utils_2.12.3 # [141] RCurl_1.98-1.16 memoise_2.0.1 # [143] rmarkdown_2.29 pheatmap_1.0.12 # [145] scales_1.3.0 R.methodsS3_1.8.2 # [147] googlesheets4_1.0.1 future_1.34.0 # [149] RANN_2.6.2 Cairo_1.6-2 # [151] spatstat.data_3.1-4 rstudioapi_0.17.1 # [153] cluster_2.1.8 spatstat.utils_3.1-2 # [155] hms_1.1.3 fitdistrplus_1.2-2 # [157] munsell_0.5.1 cowplot_1.1.3 # [159] colorspace_2.1-1 rlang_1.1.4 # [161] quadprog_1.5-8 DelayedMatrixStats_1.20.0 # [163] sparseMatrixStats_1.10.0 dotCall64_1.2 # [165] ggforce_0.4.2 xfun_0.50 # [167] coda_0.19-4.1 TH.data_1.1-2 # [169] modelr_0.1.10 remotes_2.5.0 # [171] modeltools_0.2-23 abind_1.4-8 # [173] libcoin_1.0-10 treeio_1.22.0 # [175] ggsci_3.2.0 futile.options_1.0.1 # [177] bitops_1.0-9 ps_1.8.1 # [179] promises_1.3.2 sandwich_3.1-1 # [181] reprex_2.0.2 DelayedArray_0.24.0 # [183] compiler_4.4.2 prettyunits_1.2.0 # [185] beachmat_2.14.0 SparseM_1.84-2 # [187] listenv_0.9.1 Rcpp_1.0.14 # [189] edgeR_4.4.1 BiocSingular_1.14.0 # [191] tensor_1.5 MASS_7.3-64 # [193] BiocParallel_1.40.0 spatstat.random_3.0-1 # [195] R6_2.5.1 fastmap_1.2.0 # [197] multcomp_1.4-20 fastmatch_1.1-6 # [199] rstatix_0.7.1 vipor_0.4.7 # [201] ROCR_1.0-11 rsvd_1.0.5 # [203] gtable_0.3.6 phangorn_2.10.0 # [205] KernSmooth_2.23-26 miniUI_0.1.1.1 # [207] deldir_2.0-4 htmltools_0.5.8.1 # [209] RcppParallel_5.1.9 spatstat.explore_3.0-5 # [211] lifecycle_1.0.4 processx_3.8.5 # [213] pryr_0.1.6 callr_3.7.6 # [215] sass_0.4.9 vctrs_0.6.5 # [217] spatstat.geom_3.0-3 haven_2.5.4 # [219] ggfun_0.0.8 future.apply_1.11.3 # [221] bslib_0.9.0 pillar_1.10.1 # [223] gplots_3.1.3 magick_2.8.5 # [225] metapod_1.14.0 locfit_1.5-9.10 # [227] jsonlite_1.8.9 argparse_2.2.5 # [229] GetoptLong_1.0.5 参考：https://github.com/broadinstitute/infercnv https://github.com/harbourlab/UPhyloplot2 https://pmc.ncbi.nlm.nih.gov/articles/PMC7210975/ https://pmc.ncbi.nlm.nih.gov/articles/PMC8122019/ https://github.com/AntonioDeFalco/SCEVAN https://kharchenkolab.github.io/numbat/articles/numbat.html