其他
CNS图表复现17—inferCNV结果解读及利用之进阶
分享是一种态度
前言
我们的CNS图表复现之旅已经开始,你可以点击图表复现话题回顾。如果你感兴趣也想加入交流群,可以去:你要的rmarkdown文献图表复现全套代码来了(单细胞)找到我们的拉群小助手哈。
前面我们提到的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。
如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦!
长按扫码可关注