单细胞转录组鉴定肿瘤细胞：CopyKAT和InferCNV
如何区分肿瘤微环境中的正常细胞和恶性细胞，相信是每个肿瘤单细胞研究者非常关注的问题。
这儿我们介绍两个单细胞CNV分析的软件，一个是2021年1月发在nature biotechnology上的CopyKAT，另一个是之前很多单细胞paple用到R包InferCNV。
CopyKAT (Copynumber Karyotyping of Tumors) 是一种计算工具，使用综合贝叶斯方法以 5MB 分辨率识别单细胞中的全基因组非整倍性，以使用高通量 sc-RNAseq 数据将肿瘤细胞与正常细胞和肿瘤亚克隆分离。从 RNAseq 数据计算 DNA 拷贝数事件的基本逻辑是，许多相邻基因的基因表达水平可以提供深度信息来推断该区域的基因组拷贝数。 CopyKAT 估计的拷贝数图谱可以与通过全基因组 DNA 测序获得的实际 DNA 拷贝数达到高度一致 (80%)。预测肿瘤/正常细胞状态的基本原理是非整倍体在人类癌症中很常见 (90%)。具有广泛全基因组拷贝数畸变（非整倍体）的细胞被认为是肿瘤细胞，而基质正常细胞和免疫细胞通常具有 2N 二倍体或接近二倍体的拷贝数分布。
https://github.com/navinlabcode/copykat
InferCNV 用于探索肿瘤单细胞 RNA-Seq 数据，以识别体细胞大规模染色体拷贝数改变的证据，例如整个染色体或大片段染色体的获得或缺失。这是通过与一组参考“正常”细胞相比，探索跨肿瘤基因组位置的基因表达强度来完成的。生成的热图显示了每条染色体的相对表达强度，并且与正常细胞相比，肿瘤基因组的哪些区域过多或较少，这通常变得很明显。
https://bioconductor.org/packages/devel/bioc/vignettes/infercnv/inst/doc/inferCNV.html
copyKAT相较于inferCNV是不需要正常细胞的，可以自动寻找diploid细胞作为正常细胞。不过原理都是通过单细胞转录组数据来推断细胞的染色体倍数，进而推断是正常细胞还是肿瘤细胞。
相关包安装
##安装copykat
library(devtools)
install_github("navinlabcode/copykat")
library("copykat") #安装成功可以简单测试这个包
copykat.test <- copykat(rawmat=exp.rawdata, sam.name="test")
##安装InferCNV
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("infercnv")
#如果没有安装Seurat的，最好装一个Seurat
# conda install -c bioconda r-seurat
##我安装的时候一个依赖包一直报错，用conda重新创建了一个环境，才安装成功
在这里，我们还是使用我们前面经常使用的pbmc3k数据集先做测试，看看两个包的运算性能和运算结果
pbmc数据集相关下载，seurat聚类都可参照前面的简书：https://www.jianshu.com/p/adda4536b2cb
copykat 使用
library(Seurat)
library(copykat)
library(tidyverse)
rm(list = ls())
setwd("/home/wucheng/jianshu/function/data")
pbmc <-readRDS("pbmc.rds") ##这儿我们直接导入Seurat标准化，聚类的pbmc数据
current.cluster.ids <- c(0:8)
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
pbmc@meta.data$celltype <- plyr::mapvalues(x = pbmc@meta.data[,"seurat_clusters"], from = current.cluster.ids, to = new.cluster.ids)
head(pbmc@meta.data)
scRNA <- pbmc
counts <- as.matrix(scRNA@assays$RNA@counts)
setwd("/home/wucheng/jianshu/pbmc/copykat")
cnv <- copykat(rawmat=counts, ngene.chr=5, sam.name="SCLC", n.cores=8)
# ngene.chr参数是过滤细胞的一个标准，它要求被用来做CNV预测的细胞，一个染色体上至少有5个基因。
# sam.name定义样本名称 (sample name)，会给出来的文件加前缀
saveRDS(cnv, "cnv.rds")
# SCLC_copykat_prediction.txt是上一步生成的结果
mallignant <- read.delim("SCLC_copykat_prediction.txt", row.names = 1)
# 把细胞的良恶性信息加入metadata
scRNA <- AddMetaData(scRNA, metadata = mallignant)
p1 <- DimPlot(scRNA, group.by = "celltype", label = T)
p2 <- DimPlot(scRNA, group.by = "copykat.pred") + scale_color_manual(values = c("red", "gray"))
pc <- p1 + p2
ggsave("pred_mallignant.pdf", pc, width = 12, height = 5)
pbmc_copykat
pbmc_copykat
可以看到部分细胞被预测出拷贝数畸变（非整倍体），但是实际同为正常细胞变化其实是不明显的。
InferCNV使用
获取InferCNV所需文件，主要包含一个基因表达矩阵，meta矩阵和一个基因位置矩阵
library(infercnv)
library(Seurat)
pbmc <-readRDS("/home/wucheng/jianshu/function/data/pbmc.rds")
current.cluster.ids <- c(0:8)
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
pbmc@meta.data$celltype <- plyr::mapvalues(x = pbmc@meta.data[,"seurat_clusters"], from = current.cluster.ids, to = new.cluster.ids)
head(pbmc@meta.data)
#第二个文件:单细胞RNA-seq表达量的原始矩阵
dfcount = as.data.frame(pbmc@assays$RNA@counts)
# 第二个文件:注释文件，记录肿瘤和正常细胞
groupinfo= data.frame(cellId = colnames(dfcount),cellType= pbmc@meta.data$celltype )
# 第三文件基因注释文件
library(AnnoProbe)
geneInfor=annoGene(rownames(dfcount),"SYMBOL",'human')
geneInfor=geneInfor[with(geneInfor, order(chr, start)),c(1,4:6)]
geneInfor=geneInfor[!duplicated(geneInfor[,1]),]
## 这里可以去除性染色体
# 也可以把染色体排序方式改变
dfcount =dfcount [rownames(dfcount ) %in% geneInfor[,1],]
dfcount =dfcount [match( geneInfor[,1], rownames(dfcount) ),]
# 输出
setwd("/home/wucheng/jianshu/pbmc/InferCNV")
write.table(dfcount ,file ='expFile.txt',sep = '\t',quote = F)
write.table(groupinfo,file = 'metaFiles.txt',sep = '\t',quote = F,col.names = F,row.names = F)
write.table(geneInfor,file ='geneFile.txt',sep = '\t',quote = F,col.names = F,row.names = F)
运行InferCNV，获取CNV表达热图
##
library(Seurat)
library(ggplot2)
library(infercnv)
infercnv_obj = CreateInfercnvObject(raw_counts_matrix="/home/wucheng/jianshu/pbmc/InferCNV/expFile.txt",
annotations_file="/home/wucheng/jianshu/pbmc/InferCNV/metaFiles.txt",
delim="\t",
gene_order_file= "/home/wucheng/jianshu/pbmc/InferCNV/geneFile.txt",
ref_group_names=NULL) # 如果有正常细胞的话，把正常细胞的分组填进去
future::plan("multiprocess",workers=15)# 多核并行处理
infercnv_obj = infercnv::run(infercnv_obj,
cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
out_dir='infercnv_out/',
cluster_by_groups=TRUE, # 是否根据细胞注释文件的分组对肿瘤细胞进行
denoise=TRUE, #去噪
HMM=TRUE, #是否基于HMM预测CNV,选择F会加快运算速度
output_format = "pdf")
infercnv_heatmap
我们可以发现正常细胞之间其实同样存在拷贝数变化。InferCNV不会直接判断哪些细胞为另一个状态（恶性）的细胞。
几千个细胞两个方法运行的都比较快，可以调整运行的线程数。
思考：这儿是否说明这些算法都不够精准，如果把拷贝数变化的细胞作为肿瘤细胞，那正常细胞一部分也能预测为肿瘤细胞，这不符合实际？？？
可能算法适用肿瘤数据，来预测肿瘤细胞，我们使用一个肿瘤数据试试。
我们从copykat提供的数据集中，下载了DCIS表达谱GSM4476485（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4476485）
setwd("/home/wucheng/jianshu/pbmc/TEST")
counts <- read.table("/home/wucheng/jianshu/pbmc/GSM4476485_combined_UMIcount_CellTypes_DCIS1.txt",header=T,sep="\t",check.names=F)
cnv <- copykat(rawmat=counts, ngene.chr=5, sam.name="SCLC", n.cores=8)
# ngene.chr参数是过滤细胞的一个标准，它要求被用来做CNV预测的细胞，一个染色体上至少有5个基因。
# sam.name定义样本名称 (sample name)，会给出来的文件加前缀
saveRDS(cnv, "cnv.rds")
##前面预测了每个细胞的拷贝数
pred.test <- data.frame(cnv$prediction)
CNA.test <- data.frame(cnv$CNAmat)
head(CNA.test[ , 1:5])
head(pred.test)
##对单细胞样本进行普通聚类分析
scRNA = CreateSeuratObject(counts[c(-1,-2),])
scRNA = NormalizeData(scRNA)
scRNA = FindVariableFeatures(scRNA)
scRNA = ScaleData(scRNA)
scRNA = RunPCA(scRNA,npcs = 50)
scRNA = RunTSNE(scRNA, reduction = "pca", dims = 1:30)
scRNA = FindNeighbors(scRNA,reduction = "pca",dims = 1:30)
scRNA = FindClusters(scRNA,resolution =0.8)
##copykat预测的细胞状态添加到单细胞meta表中
scRNA@meta.data$copykat.pred <- pred.test$copykat.pred
tumor.cells <- pred.test$cell.names[which(pred.test$copykat.pred=="aneuploid")] #这儿可以把恶性细胞聚类，识别亚克隆，这也是copykat提到的可以识别恶性细胞亚克隆功能
tumor.mat <- CNA.test[, which(colnames(CNA.test) %in% tumor.cells)]
hcc <- hclust(parallelDist::parDist(t(tumor.mat),threads =4, method = "euclidean"), method = "ward.D2")
hc.umap <- cutree(hcc,2) ##这儿把恶性细胞分为了两个簇
scRNA@meta.data$copykat.tumor.pred <- rep("normal", nrow(scRNA@meta.data))
scRNA@meta.data$copykat.tumor.pred[rownames(scRNA@meta.data) %in% names(hc.umap[hc.umap==1])] <- "tumor cluster 1"
scRNA@meta.data$copykat.tumor.pred[rownames(scRNA@meta.data) %in% names(hc.umap[hc.umap==2])] <- ""##把亚克隆信息也添加到单细胞meta表中
p1 <- DimPlot(scRNA, label = T)
p2 <- DimPlot(scRNA, group.by = "copykat.pred")
p3 <- DimPlot(scRNA, group.by = "copykat.tumor.pred")
pc <-p1 + p2 + p3
ggsave("pred_celltype.pdf", pc, width = 15, height = 5)
这儿把恶性细胞重新分为不同亚克隆tumor cluster 1，tumor cluster 2，后面实际是可以比较不同亚克隆类别间的基因功能差异的
copykat
copykat_heatmap
可以看到分的几个细胞簇中，恶性细胞和非恶性细胞还是分的很开的，恶性细胞只出现某些簇中。热图也显示了恶性细胞相比非恶性细胞一些区域有拷贝数的扩增或者缺失。说明这个方法在它们数据集中性能还是不错的。
InferCNV运行和之前类似，这儿不做重复了。
有趣的是，把copykat用到了我们测序的一套数据集中，具体数据集我就不说了，这儿包含了正常和肿瘤样本。
image.png
image.png
N为正常样本，S，T为肿瘤样本，同一个病人。
如果将正常癌旁和肿瘤一起分析运行copykat，我们发现恶性主要集中在上皮，不过正常样本中也会存在一部分细胞预测为恶性细胞，这其实是不符合实际的
index1 <-which(scRNA@meta.data[,"copykat.pred"]=="aneuploid")
index2 <-which(scRNA@meta.data[,"copykat.pred"]=="diploid")
table(scRNA@meta.data[index1,"orig.ident"],scRNA@meta.data[index1,"cell"]) ##恶性
B cell Endothelial Epithelial Fibroblast Myeliod cell T cell
N 4 2 678 3 0 11
S 10 3 2086 3 7 6
T 3 6 605 7 7 4
table(scRNA@meta.data[index2,"orig.ident"],scRNA@meta.data[index2,"cell"]) ##非恶性
B cell Endothelial Epithelial Fibroblast Myeliod cell T cell
N 788 121 706 72 115 281
S 885 55 1107 90 352 529
T 594 987 728 307 473 521
或许设定对照（已知的正常细胞），或者干脆只运行肿瘤样本，来预测肿瘤样本中存在的恶性细胞，是一个更好的选择。
注意：儿童肿瘤和血液肿瘤中基本没有拷贝数改变，不太适合copyKAT或inferCNV寻找恶性细胞。
想知道单个细胞为正常或者恶性，copykat会更简单容易很多；如果利用InferCNV预测，一方面最好有正常细胞对照，然后后面需要做一个对每个细胞进行拷贝数打分，判断恶性程度，会复杂一些。
希望大家点赞关注，谢谢！！！！！！！！！！！！
有不清楚的可以一起探讨，谢谢！！！！
