查看原文
其他

CNS图表复现17—inferCNV结果解读及利用之进阶

生信技能树 单细胞天地 2022-06-06

分享是一种态度

前言

我们的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。



如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程



看完记得顺手点个“在看”哦!


生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存