copyKAT推断单细胞转录组肿瘤细胞CNV（自动识别肿瘤normal和tumor）

hi，各位，2021年1月份发表了一篇文献Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes ,2021年1月发表于nature biotechnology，影响因子36分多，非常高，前几位的作者都是中国人，非常值得一看，文章中推出新的推断肿瘤细胞CNV的方法：copyKAT，也是通过单细胞转录组数据来推断细胞的染色体倍数，进而推断是正常细胞（diploid）还是肿瘤细胞（aneuploid）。它还可以进一步对肿瘤细胞进行聚类，找出不同的亚群。

我们来看一下原理图

首先对基因表达量做标准化并稳定其方差（a）。相较于inferCNV（见之前的分享：inferCNV），copyKAT可以自动寻找diploid cells作为正常细胞（b）。对每个非正常细胞，利用MCMC寻找其CNV的断点（breakpoints）并得到segments（c）。正常细胞和肿瘤细胞由于其基因表达量分布的不同可以被分开（d）。肿瘤细胞通常还可以继续聚类得到其亚群（e）。

接下来我们来运行一下这个代码，看看结果如何

加载模块并读取数据（以10X数据为例）

library(Seurat)
library(copykat)
raw <- Read10X(data.dir = data.path.to.cellranger.outs)
raw <- CreateSeuratObject(counts = raw, project = "copycat.test", min.cells = 0, min.features = 0)
exp.rawdata <- as.matrix(raw@assays$RNA@counts)

running copykat

已经准备好了输入矩阵，原始的UMI计数矩阵，已经准备运行copykat。 cellranger输出中的默认基因ID是基因符号，因此输入"Symbol”或“ S”。 为了滤除细胞，需要每个染色体至少5个基因来计算DNA拷贝数。 可以将其调低到ngene.chr = 1以保持尽可能多的细胞，但是使用至少5个基因来代表一个染色体并不是很严格。 为了滤除基因，可以调整参数以仅保留在细胞的LOW.DR至UP.DR部分中表达的基因。 我把默认值LOW.DR = 0.05，UP.DR = 0.2。 我可以调低这些值以在分析中保留更多基因。 我需要确保LOW.DR小于UP.DR。
要求copykat每个片段至少接受25个基因。 我可以尝试使用其他选项，每个bin的范围为15-150个基因。 KS.cut是分段参数，范围从0到1。增加KS.cut会降低灵敏度，即减少分段/断点。 通常它在0.05-0.15的范围内工作。 (还记得inferCNV每个bin多少个基因么？童鞋们)
One struggling and intesting observation is that none of one clustering method could fit all datasets. In this version, I add a distance parameters for clustering that include "euclidean" distannce and correlational distance, ie. 1-"pearson" and "spearman" similarity. In general, corretional distances tend to favor noisy data, while euclidean distance tends to favor data with larger CN segments.（确实是这样）。
I add an option to input known normal cell names as a vector object. Default is NULL.（我们这里就使用默认参数）。

我们先来看一下这个函数的参数

help(copykat)
rawmat: raw data matrix; genes in rows; cell names in columns.
id.type: gene id type: Symbol or Ensemble.
ngene.chr: minimal number of genes per chromosome for cell filtering.
LOW.DR: minimal population fractoins of genes for smoothing.
UP.DR: minimal population fractoins of genes for segmentation.
win.size: minimal window sizes for segmentation.
norm.cell.names: a vector of normal cell names.
KS.cut: segmentation parameters, input 0 to 1; larger looser criteria.
sam.name: sample name.
n.cores: number of cores for parallel computing.
###参数都非常容易理解

那我们运行一下

copykat.test <- copykat(rawmat=exp.rawdata, id.type="S", ngene.chr=5, win.size=25, KS.cut=0.1, sam.name="test", distance="euclidean", norm.cell.names="", n.cores=4)
###看看有什么（得到一个列表）
pred.test <- data.frame(copykat.test$prediction)
head(pred.test)

这部分是对细胞类型是否是恶性细胞的判定。

CNA.test <- data.frame(copykat.test$CNAmat)
head(CNA.test)

拷贝数矩阵包含了基因的染色体信息和每个细胞的变异情况

###聚类信息
head(copykat.test$hclustering$)

这部分包含了聚类的顺序，距离等信息，非常有用。

热图展示

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

看起来层次分明

再对肿瘤细胞再聚类并画热图，又能分成两群。

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

展现肿瘤内部的异质性。

最后把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.pred
TNBC1@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

最后作者提到一个需要注意的点，不是所有的肿瘤都存在CNV。儿童肿瘤和血液肿瘤中基本没有copy number event，所以是不适合用这些方法（copyKAT或inferCNV）来寻找肿瘤细胞的。

看看英文版原理

The statistical workflow of CopyKAT combines a Bayesian approach with hierarchical clustering to calculate genomic copy number profiles of single cells and define clonal substructure from high-throughput 3′ scRNA-seq data
原理跟inferCNV一样

表现基本一致，更优一点。