单细胞分析十八般武艺13：CopyKAT
单细胞测序技术的发展日新月异，新的分析工具也层出不穷。每个工具都有它的优势与不足，在没有权威工具和流程的单细胞生信江湖里，多掌握几种分析方法和工具，探索数据时常常会有意想不到的惊喜。
往期专题
单细胞转录组基础分析专题
CopyKAT简介
CopyKAT是一个核心算法原理与inferCNV基本相同，但是作为后来者依然可以发表在Nat Biotechnol上的CNV分析工具。我们先来看看它的原理：
a：copyKAT首先会对输入的表达矩阵按基因在基因组的顺序排序；然后进行数据均一化处理，所用算法是多项式动态线性模型（polynomial dynamic linear modeling，DLM），感觉类似sctransform的方法。
b：接下来是基于欧式距离或相关性对细胞聚类，然后使用高斯混合模型（Gaussian mixture model，GMM）估计每个cluster的方差，方差最小的cluster用作后续分析的二倍体参考数据（即正常细胞）。如果数据中正常细胞过少，或肿瘤细胞的染色体特征接近二倍体（即CNA较少），copyKAT会增加误判的概率，这种情况下建议使用 ‘GMM definition’ 模式逐个识别二倍体的正常细胞。
c：copyKAT通过基因表达量计算CNA时，每25个基因合并成一个检测窗口，因此最终的CNA分辨率是5M（基因组上基因之间平均间距20kb*25基因/窗口）。copyKAT会计算相邻检测窗口的基因表达量均值之间的显著性，存在显著性差异的窗口之间识别为染色体断点（chromosome breakpoints）。
d：使用CNA数据进行层次聚类或GMM（适用于正常细胞与肿瘤细胞CNA差异不大的数据）区分非整倍体肿瘤细胞和二倍体正常细胞。
e：最后通过对肿瘤细胞进一步聚类识别肿瘤亚克隆。
使用copyKAT几个月以来，我的感觉是它比inferCNV更快、更准，使用起来也更方便，只需输入单细胞测序的原始UMI计数矩阵就可以了。
安装与使用
library(devtools)install_github("navinlabcode/copykat")# 使用包自带的数据测试copykat.test <- copykat(rawmat=exp.rawdata, sam.name="test")
## 运行步骤如下：#[1] "running copykat v1.0.5 updated 07/15/2021"#[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 3.175674 mins参数说明
copykat.test <- copykat( rawmat=exp.rawdata, # 原始表达矩阵，TPM数据也能使用 id.type="S", # 基因名称的类型，S代表symbol ngene.chr=5, # 过滤细胞的标准，要求每条染色体上最少有5个基因，实测会过滤不少细胞 win.size=25, # 检测窗口的基因数，默认25个基因，可以在15-150的范围内选择 KS.cut=0.15, # 差异检验的阈值，越大灵敏度越低 sam.name="test", # 样本名称，copykat运行的结果文件会以此为前缀 distance="euclidean", # 聚类的距离，可选c("euclidean", "pearson", "spearman") norm.cell.names="", # 正常细胞的名称向量 n.cores=16) # 线程数量数据实测
library(copykat)download.file("https://cf.10xgenomics.com/samples/cell-exp/6.0.0/Brain_Tumor_3p/Brain_Tumor_3p_filtered_feature_bc_matrix.h5", destfile = "matrix.h5")mat <- Read10X_h5("matrix.h5")sco <- CreateSeuratObject(mat, project = "Glioblastoma")Glioblastoma <- copykat(rawmat = sco@assays$RNA@counts, sam.name = "Glioblastoma", n.cores = 18)大约20分钟就得到了分析结果
常见问题
CopyKAT会直接生成一个肿瘤细胞判断结果的文件copykat_prediction.txt，用过的朋友很容易发现这个结果中的细胞数目与我们输入的细胞数目不同。
dim(sco)#36601 1615 #输入copykat的细胞数为1615个pred <- data.table::fread("Glioblastoma_copykat_prediction.txt")nrow(pred)# 1299 #有预测结果的细胞只有1299个table(pred$copykat.pred)#aneuploid diploid # 512 787这是因为copykat会自动过滤低质量细胞，它的默认判定标准是每条染色体上最少有5个基因表达的细胞予以保留。如果觉得这个标准太严了，可以通过ngene.chr参数调整，例如ngene.chr=3就可以少过滤一些细胞。
肿瘤细胞分亚型
### 肿瘤细胞分亚型cnv.scores <- data.table::fread("Glioblastoma_copykat_CNA_results.txt", data.table = F)colnames(cnv.scores) <- gsub(".", "-", colnames(cnv.scores), fixed = T)pred.res <- data.table::fread("Glioblastoma_copykat_prediction.txt")tumor.cells <- pred.res$cell.names[which(pred.res$copykat.pred=="aneuploid")]tumor.mat <- cnv.scores[, tumor.cells]hcc <- hclust(parallelDist::parDist(t(tumor.mat),threads =10, 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)
my_palette <- colorRampPalette(rev(RColorBrewer::brewer.pal(n = 3, name = "RdBu")))(n = 999)chr <- as.numeric(cnv.scores$chrom) %% 2+1rbPal1 <- colorRampPalette(c('black','grey'))CHR <- rbPal1(2)[as.numeric(chr)]chr1 <- cbind(CHR,CHR)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))
png("Glioblastoma_subcluster.png", width = 1500, height = 900)heatmap.3(t(tumor.mat), dendrogram="r", distfun = function(x) parallelDist::parDist(x, threads =10, 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=rbPal6, cex=0.9, bty='n')dev.off()这个数据没有明显的CNV亚型
如果您阅读此文有所疑惑，或有不同见解，亦或其他问题，欢迎添加我的微信探讨。我工作的公司2016年开始单细胞测序服务，至今已完成近万例样本的单细胞测序，服务质量经过了10X Genomics公司的权威认证。欢迎大家联系Kinesin洽谈单细胞测序及数据分析业务！