前面我们提到的inferCNV结果里面的CNV热图，虽然说可以肉眼简单看看不同的细胞亚群是否具有大面积的CNV事件来判定它是否是恶性细胞，但是这样的判定具有很大程度的主观性，所以理论上需要有更好的方法，也就是计算具体每个细胞的CNV score 。

读取infercnv.observations.txt文件来计算CNV score

代码如下：

cnv_table <- read.table("plot_out/inferCNV_output2/infercnv.observations.txt", header=T)
# Score cells based on their CNV scores 
# Replicate the table 
cnv_score_table <- as.matrix(cnv_table)
cnv_score_mat <- as.matrix(cnv_table)
# Scoring
# CNV的5分类系统
cnv_score_table[cnv_score_mat > 0 & cnv_score_mat < 0.3] <- "A" #complete loss. 2pts
cnv_score_table[cnv_score_mat >= 0.3 & cnv_score_mat < 0.7] <- "B" #loss of one copy. 1pts
cnv_score_table[cnv_score_mat >= 0.7 & cnv_score_mat < 1.3] <- "C" #Neutral. 0pts
cnv_score_table[cnv_score_mat >= 1.3 & cnv_score_mat <= 1.5] <- "D" #addition of one copy. 1pts
cnv_score_table[cnv_score_mat > 1.5 & cnv_score_mat <= 2] <- "E" #addition of two copies. 2pts
cnv_score_table[cnv_score_mat > 2] <- "F" #addition of more than two copies. 2pts

# Check
table(cnv_score_table[,1])
# Replace with score 
cnv_score_table_pts <- cnv_table
rm(cnv_score_mat)
# 把5分类调整为 3分类系统
cnv_score_table_pts[cnv_score_table == "A"] <- 2
cnv_score_table_pts[cnv_score_table == "B"] <- 1
cnv_score_table_pts[cnv_score_table == "C"] <- 0
cnv_score_table_pts[cnv_score_table == "D"] <- 1
cnv_score_table_pts[cnv_score_table == "E"] <- 2
cnv_score_table_pts[cnv_score_table == "F"] <- 2

# Scores are stored in “cnv_score_table_pts”. Use colSums to add up scores for each cell and store as vector 
cell_scores_CNV <- as.data.frame(colSums(cnv_score_table_pts))
colnames(cell_scores_CNV) <- "cnv_score"
head(cell_scores_CNV)
write.csv(x = cell_scores_CNV, file = "cnv_scores.csv")

因为infercnv.observations.txt文件超级大，而且上面的代码很多地方并不高效，所以运行耗时很长哦！

不同的细胞恶性与否判定方法的量化

载入前面靠细胞分群是否跨越病人，以及全局CNV层次聚类，这两个策略判定的细胞恶性与否，目前就可以看CNV score啦

load(file = 'phe-of-cancer-or-not.Rdata')
phe=phe[rownames(phe) %in% rownames(cell_scores_CNV),]
infercnv.labels <- cutree(infercnv.dend, k = 6, order_clusters_as_data = FALSE)
phe$inferCNV= infercnv.labels[match(rownames(phe), names(infercnv.labels) )]
phe$cnv_scores = cell_scores_CNV[rownames(phe),]
table(phe$cancer,phe$inferCNV) 
library(ggpubr)
p1=ggboxplot(phe,'cancer','cnv_scores', fill = "cancer")
p2=ggboxplot(phe,'inferCNV','cnv_scores', fill = "inferCNV")
library(patchwork)
p1+p2

出图如下：

这个时候就很明显了，这个CNV score可以代替我们的肉眼来对细胞亚群的恶性与否进行判定。

 1 2 3 4 5 6
 cancer 948 956 284 640 466 438
 non-cancer 6 1480 1 4 219 2

也就是说，我能判定为恶性细胞的，应该是即属于我们的inferCNV热图层次聚类的第二细胞亚群，也要属于我们的跨越病人聚类的细胞亚群。

最后仅仅是那 1480 个细胞被我判定为非恶性细胞。文章也是说epithelial cells (n = 5,581), 里面有3,754 cancer cells。