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" )