hi募舟，各位筑舅，2021年1月份發(fā)表了一篇文獻Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes ,2021年1月發(fā)表于nature biotechnology茵典，影響因子36分多，非常高邪狞，前幾位的作者都是中國人挖垛，非常值得一看，文章中推出新的推斷腫瘤細胞CNV的方法：copyKAT婉支，也是通過單細胞轉錄組數(shù)據(jù)來推斷細胞的染色體倍數(shù)鸯隅，進而推斷是正常細胞（diploid）還是腫瘤細胞（aneuploid）。它還可以進一步對腫瘤細胞進行聚類向挖，找出不同的亞群蝌以。

我們來看一下原理圖

24927444-5fbf08d4fc6ef6d0.png

接下來我們來運行一下這個代碼蒜危，看看結果如何

加載模塊并讀取數(shù)據(jù)（以10X數(shù)據(jù)為例）

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

已經(jīng)準備好了輸入矩陣袁串，原始的UMI計數(shù)矩陣箕别，已經(jīng)準備運行copykat。 cellranger輸出中的默認基因ID是基因符號，因此輸入"Symbol”或“ S”源譬。 為了濾除細胞，需要每個染色體至少5個基因來計算DNA拷貝數(shù)。 可以將其調(diào)低到ngene.chr = 1以保持盡可能多的細胞，但是使用至少5個基因來代表一個染色體并不是很嚴格捎迫。 為了濾除基因，可以調(diào)整參數(shù)以僅保留在細胞的LOW.DR至UP.DR部分中表達的基因表牢。 我把默認值LOW.DR = 0.05窄绒，UP.DR = 0.2。 我可以調(diào)低這些值以在分析中保留更多基因崔兴。 我需要確保LOW.DR小于UP.DR彰导。
要求copykat每個片段至少接受25個基因。 我可以嘗試使用其他選項敲茄，每個bin的范圍為15-150個基因位谋。 KS.cut是分段參數(shù)，范圍從0到1堰燎。增加KS.cut會降低靈敏度掏父，即減少分段/斷點。 通常它在0.05-0.15的范圍內(nèi)工作秆剪。 (還記得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.（我們這里就使用默認參數(shù)）仅讽。

我們先來看一下這個函數(shù)的參數(shù)

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.
###參數(shù)都非常容易理解

那我們運行一下

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)

圖片.png

這部分是對細胞類型是否是惡性細胞的判定陶缺。

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

圖片.png

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

圖片.png

這部分包含了聚類的順序，距離等信息洁灵，非常有用饱岸。

熱圖展示

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

圖片.png

看起來層次分明

再對腫瘤細胞再聚類并畫熱圖，又能分成兩群徽千。

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

圖片.png

展現(xiàn)腫瘤內(nèi)部的異質(zhì)性苫费。

最后把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

圖片.png

最后作者提到一個需要注意的點，不是所有的腫瘤都存在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一樣

圖片.png