本文介绍了一种创新的计算方法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" )