本文共分3部分。
首先介绍一下copyKAT原理和使用方法。
然后提出一个关于细胞过滤的问题。
最后提出一个解决问题的方法供参考。
1.copyKAT及使用简介
copyKAT是一种推断肿瘤细胞CNV的方法,和inferCNV类似。文章于2021年1月发在nature biotechnology,链接见文末。也是大名鼎鼎的期刊。其原理同样是通过单细胞转录组数据来推断细胞的染色体倍数,进而推断是正常细胞(diploid)还是肿瘤细胞(aneuploid)。
1.1原理及工作流程图
image
(a)copyKAT首先对基因表达量做标准化并稳定其方差。
(b)它可以自动寻找diploid cells作为正常细胞。
(c)对每个非正常细胞,利用MCMC寻找其CNV的断点(breakpoints)并得到segments。
(d)正常细胞和肿瘤细胞由于其基因表达量分布的不同可以被分开。
(e)肿瘤细胞通常还可以继续聚类得到其亚群。
1.2正常运行示例
我们首先导入示例数据运行一下,用文章里的一个TNBC1数据(下载链接:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4476486)。运行非常简单,缺点是慢,且windows只能用n.cores=1来运行,建议使用服务器来做。我用16G内存电脑跑了一次5000个细胞花费了8h。4000个细胞5h。
exp.rawdata <- read.table("./copykat/data/GSM4476486_filtered_UMIcount_TNBC1.txt", header=T, sep='\t', check.names = F)copykat.test <- copykat(rawmat=exp.rawdata, id.type="S", cell.line="no", ngene.chr=5, win.size=25, KS.cut=0.15, sam.name="TNBC1", distance="euclidean", n.cores=4)pred.test <- data.frame(copykat.test$prediction)CNA.test <- data.frame(copykat.test$CNAmat)
看一下结果,主要有两个,
1)预测的结果正常细胞(diploid)还是肿瘤细胞(aneuploid),
2) 每个CNV segment在每个细胞的表达量。
这里看出和inferCNV的不同来了,是基于genomic coordinate而不是gene level的表达量。。
head(pred.test) cell.names copykat.pred#AAACCTGCACCTTGTC AAACCTGCACCTTGTC aneuploid#AAACGGGAGTCCTCCT AAACGGGAGTCCTCCT diploid#AAACGGGTCCAGAGGA AAACGGGTCCAGAGGA aneuploid#AAAGATGCAGTTTACG AAAGATGCAGTTTACG aneuploid#AAAGCAACAGGAATGC AAAGCAACAGGAATGC aneuploid#AAAGCAATCGGAATCT AAAGCAATCGGAATCT aneuploidhead(CNA.test[,1:5]) chrom chrompos abspos AAACCTGCACCTTGTC AAACGGGAGTCCTCCT#1 1 1042457 1042457 -0.03206638 0.03170166#2 1 1265484 1265484 -0.03206638 0.03170166#3 1 1519859 1519859 -0.03206638 0.03170166#4 1 1826619 1826619 -0.03206638 0.03170166#5 1 2058465 2058465 -0.03206638 0.03170166#6 1 2280372 2280372 -0.03206638 0.03170166
画个热图看看。很明显正常细胞和肿瘤细胞分开了。
my_palette <- colorRampPalette(rev(RColorBrewer::brewer.pal(n = 3, name = "RdBu")))(n = 999)chr <- as.numeric(CNA.test$chrom) %% 2+1rbPal1 <- colorRampPalette(c('black','grey'))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.predpred <- 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"), 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")
![](//image109.360doc.cn/DownloadImg/2022/10/0820/253279507_2_20221008082100289.png)
image
再对肿瘤细胞再聚类并画热图,又能分成两群。
tumor.cells <- pred.test$cell.names[which(pred.test$copykat.pred=="aneuploid")]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)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')
image
最后把CNV的结果投射到单细胞聚类结果上看一看是否合理,Seurat标准流程走一遍,聚类结果和copyKAT分群结果投射到TSNE上。
standard10X = function(dat,nPCs=50,res=1.0,verbose=FALSE){ srat = CreateSeuratObject(dat) srat = NormalizeData(srat,verbose=verbose) srat = ScaleData(srat,verbose=verbose) srat = FindVariableFeatures(srat,verbose=verbose) srat = RunPCA(srat,verbose=verbose) srat = RunTSNE(srat,dims=seq(nPCs),verbose=verbose) srat = FindNeighbors(srat,dims=seq(nPCs),verbose=verbose) srat = FindClusters(srat,res=res,verbose=verbose) return(srat)}TNBC1 <- standard10X(exp.rawdata, nPCs=30, res=0.6)TNBC1@meta.data$copykat.pred <- pred.test$copykat.predTNBC1@meta.data$copykat.tumor.pred <- rep("normal", nrow(TNBC1@meta.data))TNBC1@meta.data$copykat.tumor.pred[rownames(TNBC1@meta.data) %in% names(hc.umap[hc.umap==1])] <- "tumor cluster 1"TNBC1@meta.data$copykat.tumor.pred[rownames(TNBC1@meta.data) %in% names(hc.umap[hc.umap==2])] <- "tumor cluster 2"p1 <- DimPlot(TNBC1, label = T)p2 <- DimPlot(TNBC1, group.by = "copykat.pred")p3 <- DimPlot(TNBC1, group.by = "copykat.tumor.pred")p1 + p2 + p3
image
从免疫细胞和肿瘤细胞的标记基因表达来看,copyKAT可以正确找出正常细胞和肿瘤细胞。
FeaturePlot(TNBC1,features=c("PTPRC", "EPCAM"), order = T)
image
最后作者提到一个需要注意的点,不是所有的肿瘤都存在CNV。儿童肿瘤和血液肿瘤中基本没有copy number event,所以是不适合用这些方法(copyKAT或inferCNV)来寻找肿瘤细胞的。
2. 错误示例
此处不放具体数据。遇到这个问题的人一看就懂,遇不到的就更不用看了。说明你的数据很鲁棒,棒的根本过滤不掉啊。
简单说就是我从seurat结果拿来已经过滤好的100个细胞,我希望这100个细胞都被推断出肿瘤细胞和非肿瘤细胞,
但是这个包会把我按照自己的标准选好的细胞再过滤一遍,导致不到100个细胞,然后最后就没办法和我们原来seurat画出来的tsne聚类图很好的overlap。
如下,我输入了4422个细胞,希望不要过滤全部推断CNV.
image.png
但是copyKAT之后只剩下4364个细胞被用来做后续推断,这不是我想要的结果。我试图改参数,发现所有参数改到最低0,也无法保证所有细胞被保留。于是尝试去看了源代码。。。发现代码中是有几处过滤的。见下一部分。
image.png
image.png
3.解决方案
我找到了源代码,将过滤细胞的代码行注释掉或删除即可完全保留细胞。
这里说明下(注释掉或删除)。建议注释,这是一个好习惯,可以保留代码修改记录,方便以后修改和理解。
有人可能还不知道如何查看源码,对着要看的函数按下F2即可。
过滤的地方有3处。如下,找到这三处注释掉,运行下函数就可以保留全部细胞了,不信你试试?!供大家参考。不足之处,欢迎讨论学习!
代码块1
# genes.raw <- apply(rawmat, 2, function(x) (sum(x > 0))) # if (sum(genes.raw > 200) == 0) # stop("none cells have more than 200 genes") # if (sum(genes.raw < 100) > 1) { # rawmat <- rawmat[, -which(genes.raw < 200)] # print(paste("filtered out ", sum(genes.raw <= 200), # " cells with less than 200 genes; remaining ", ncol(rawmat), # " cells", sep = "")) # }
代码块1
# if (length(toRev) > 0) { # anno.mat <- anno.mat[-toRev, ] # } # ToRemov2 <- NULL # for (i in 8:ncol(anno.mat)) { # cell <- cbind(anno.mat$chromosome_name, anno.mat[, i]) # cell <- cell[cell[, 2] != 0, ] # if (length(as.numeric(cell)) < 5) { # rm <- colnames(anno.mat)[i] # ToRemov2 <- c(ToRemov2, rm) # } # else if (length(rle(cell[, 1])$length) < 23 | min(rle(cell[, # 1])$length) < ngene.chr) { # rm <- colnames(anno.mat)[i] # ToRemov2 <- c(ToRemov2, rm) # } # i <- i + 1 # } # if (length(ToRemov2) == (ncol(anno.mat) - 7)) # stop("all cells are filtered") # if (length(ToRemov2) > 0) { # anno.mat <- anno.mat[, -which(colnames(anno.mat) %in% # ToRemov2)] # }
代码块2
# ToRemov3 <- NULL # for (i in 8:ncol(anno.mat2)) { # cell <- cbind(anno.mat2$chromosome_name, anno.mat2[, # i]) # cell <- cell[cell[, 2] != 0, ] # if (length(as.numeric(cell)) < 5) { # rm <- colnames(anno.mat2)[i] # ToRemov3 <- c(ToRemov3, rm) # } # else if (length(rle(cell[, 1])$length) < 23 | min(rle(cell[, # 1])$length) < ngene.chr) { # rm <- colnames(anno.mat2)[i] # ToRemov3 <- c(ToRemov3, rm) # } # i <- i + 1 # } # if (length(ToRemov3) == ncol(norm.mat.relat)) # stop("all cells are filtered") # if (length(ToRemov3) > 0) { # norm.mat.relat <- norm.mat.relat[, -which(colnames(norm.mat.relat) %in% # ToRemov3)] # }
References:
https://www.nature.com/articles/s41587-020-00795-2#data-availabilityhttps://github.com/navinlabcode/copykatcopyKAT推断单细胞转录组肿瘤细胞CNV