本文是參考學(xué)習(xí)CNS圖表復(fù)現(xiàn)16—inferCNV結(jié)果解讀及利用的學(xué)習(xí)筆記?？赡芨鶕?jù)學(xué)習(xí)情況有所改動(dòng)。
前面我們提到了因?yàn)榧?xì)胞數(shù)量比較多，運(yùn)行infercnv::run的時(shí)候，下面兩個(gè)參數(shù)，都是默認(rèn)值即可：

HMM參數(shù) when set to True, runs HMM to predict CNV level (default: FALSE)
denoise If True, turns on denoising according to options below (default: FALSE)

這樣速度就超級(jí)快，可以得到如下圖表：

圖片

可以看到， 部分Fibroblasts和Endothelial_cells細(xì)胞我拿它們作為ref，理論上它們是不可能有CNV事件的，所以上面的熱圖的上半部分可以看到部分Fibroblasts和Endothelial_cells細(xì)胞都是比較整齊。而且，我在需要被檢驗(yàn)CNV的細(xì)胞里面也摻入了部分Fibroblasts和Endothelial_cells細(xì)胞，它們也是沒有CNV事件的，僅僅是上皮細(xì)胞可以看到?jīng)芪挤置鞯腃NV有無(wú)的差異。

現(xiàn)在最重要的目標(biāo)就是根據(jù)這個(gè)圖表或者說(shuō)inferCNV的結(jié)果文件，把上皮細(xì)胞里面那些惡性的部分挑選出來(lái)，可以跟文章挑選好的進(jìn)行對(duì)比?；蛘吒覀兦懊娴纳掀ぜ?xì)胞聚類的結(jié)果進(jìn)行對(duì)比！

首先查看inferCNV結(jié)果文件夾，可以看到每個(gè)步驟的中間文件，都是保存下來(lái)了的：

01_incoming_data.infercnv_obj
02_reduced_by_cutoff.infercnv_obj
03_normalized_by_depth.infercnv_obj
04_logtransformed.infercnv_obj
08_remove_ref_avg_from_obs_logFC.infercnv_obj
09_apply_max_centered_expr_threshold.infercnv_obj
10_smoothed_by_chr.infercnv_obj
11_recentered_cells_by_chr.infercnv_obj
12_remove_ref_avg_from_obs_adjust.infercnv_obj
14_invert_log_transform.infercnv_obj
15_no_subclustering.infercnv_obj

但是最重要的文件是：

infercnv.observations_dendrogram.txt
infercnv.observations.txt

解析熱圖附帶的層次聚類結(jié)果
主要是該文章的GitHub代碼，讀取自己走完inferCNV流程后的結(jié)果文件夾里面的 infercnv.observations_dendrogram.txt 文件，里面存儲(chǔ)著inferCNV的CNV熱圖的細(xì)胞的層次聚類情況：

rm(list=ls())
options(stringsAsFactors = F)
library(phylogram)
library(gridExtra)
library(grid)
require(dendextend)
require(ggthemes)
library(tidyverse)
library(Seurat)
library(infercnv)
library(miscTools)

#  Import inferCNV dendrogram
infercnv.dend <- read.dendrogram(file = "plot_out/inferCNV_output2/infercnv.observations_dendrogram.txt")
# Cut tree 
infercnv.labels <- cutree(infercnv.dend, k = 6, order_clusters_as_data = FALSE)
table(infercnv.labels)
# Color labels
the_bars <- as.data.frame(tableau_color_pal("Tableau 20")(20)[infercnv.labels])
colnames(the_bars) <- "inferCNV_tree"
the_bars$inferCNV_tree <- as.character(the_bars$inferCNV_tree)
 
infercnv.dend %>% set("labels",rep("", nobs(infercnv.dend)) )  %>% plot(main="inferCNV dendrogram") %>%
  colored_bars(colors = as.data.frame(the_bars), dend = infercnv.dend, sort_by_labels_order = FALSE, add = T, y_scale=100 , y_shift )

如下所示：

圖片

可以很清楚的看到，最大的那一類 就是第二類，有3035個(gè)細(xì)胞的群體是非惡性細(xì)胞，不僅僅是從前面的熱圖可以看到它們這些細(xì)胞基本上沒有CNV事件，而且這一群細(xì)胞里面有我們自己摻入的已知的非惡性細(xì)胞，就是各300個(gè)的Fibroblasts和Endothelial_cells細(xì)胞。

> table(infercnv.labels)
infercnv.labels
   1    2    3    4    5    6 
 954 3035  285  644  686  440 

而且可以仍然是跟之前的spike-in細(xì)胞對(duì)比：

infercnv.labels=as.data.frame(infercnv.labels)
groupFiles='groupFiles.txt'   
meta=read.table(groupFiles)
infercnv.labels$V1=rownames(infercnv.labels)
meta=merge(meta,infercnv.labels,by='V1')
table(meta[,2:3])

得到的結(jié)果如下：

> table(meta[,2:3])
            infercnv.labels
V2              1    2    3    4    5    6
  epi         954 2436  285  644  685  440
  spike-endo    0  300    0    0    0    0
  spike-fib     0  299    0    0    1    0

自己人為摻入的Fibroblasts和Endothelial_cells細(xì)胞都是在第2群，也就是非惡性細(xì)胞。而剩余的其它細(xì)胞亞群，都是惡性細(xì)胞！

根據(jù)inferCNV結(jié)果判定的細(xì)胞惡性與否的結(jié)果和普通聚類的差異

在CNS圖表復(fù)現(xiàn)09—上皮細(xì)胞可以區(qū)分為惡性與否，我分享過(guò)第1,2,7,14,21,23,25 是跨越病人的聚類情況，所以先暫時(shí)認(rèn)為他們是非惡性細(xì)胞?，F(xiàn)在我們有了inferCNV結(jié)果，就可以看看兩個(gè)策略判斷細(xì)胞惡性與否的差異：

load(file = 'first_sce.Rdata')  
load(file = 'phe-of-first-anno.Rdata')
sce=sce.first
table(phe$immune_annotation)
sce@meta.data=phe 
cells.use <- row.names(sce@meta.data)[which(phe$immune_annotation=='epi')]
length(cells.use)
sce <-subset(sce, cells=cells.use)  
sce
load(file = 'phe-of-cancer-or-not.Rdata')
sce@meta.data=phe 
colnames(sce@meta.data)
table(sce@meta.data$cancer)
sce@meta.data$inferCNV=meta[match(rownames(sce@meta.data),meta$V1),'infercnv.labels']
table(sce@meta.data$cancer,sce@meta.data$inferCNV) 
 

結(jié)果如下：

> table(sce@meta.data$cancer,sce@meta.data$inferCNV)
            # 其中inferCNV得到的亞群里面，第2群細(xì)胞是 非惡性 
                1    2    3    4    5    6
  cancer      948  956  284  640  466  438
  non-cancer    6 1480    1    4  219    2

可以看到，通過(guò)聚類是否跨越病人來(lái)區(qū)分惡性與否，跟inferCNV的一致性尚可。

因?yàn)橹澳莻€(gè)非常重要的文件是：cnv_scores.csv ，在新版的inferCNV里面是沒有的，所以我們僅僅是能根據(jù)層次聚類的劃分情況，來(lái)粗略的把第2群的3000多個(gè)細(xì)胞統(tǒng)一劃分成為非惡性。但實(shí)際上，很容易看出來(lái)：

圖片

這個(gè)第2群的3000多個(gè)細(xì)胞是可以繼續(xù)細(xì)分后，重新判定細(xì)胞的惡性與否，這樣就可以提高兩個(gè)技術(shù)的吻合性。

但是這樣仍然是“治標(biāo)不治本”，無(wú)論我們的肉眼多么厲害，僅僅是從這個(gè)CNV熱圖去判斷具體的某個(gè)層次聚類的亞群細(xì)胞是惡性與否，實(shí)在是太主觀了。

雖然文章作者是靠這樣的策略來(lái)判斷細(xì)胞惡性與否

雖然我和文章都是取全部的上皮細(xì)胞，以及部分Fibroblasts和Endothelial_cells細(xì)胞來(lái)一起運(yùn)行inferCNV流程。上皮細(xì)胞肯定都是一樣的，但是因?yàn)槭请S機(jī)函數(shù)取部分Fibroblasts和Endothelial_cells細(xì)胞，所以沒辦法保證我跟文章后面的inferCNV結(jié)果一模一樣。文章的確可以直接層次聚類后的6類群，就區(qū)分出來(lái)了惡性細(xì)胞，全部代碼如下：

#replace inferCNV.class 1 as "tumor"
inferCNV.annotation.malignant$inferCNV.class <- gsub(pattern = 1, replacement = "tumor", x = inferCNV.annotation.malignant$inferCNV.class)
#replace inferCNV.class 2 as "nontumor"
inferCNV.annotation.malignant$inferCNV.class <- gsub(pattern = 2, replacement = "nontumor", x = inferCNV.annotation.malignant$inferCNV.class)
#replace inferCNV.class 3 as "tumor"
inferCNV.annotation.malignant$inferCNV.class <- gsub(pattern = 3, replacement = "tumor", x = inferCNV.annotation.malignant$inferCNV.class)
#replace inferCNV.class 4 as "tumor"
inferCNV.annotation.malignant$inferCNV.class <- gsub(pattern = 4, replacement = "tumor", x = inferCNV.annotation.malignant$inferCNV.class)
#replace inferCNV.class 5 as "tumor"
inferCNV.annotation.malignant$inferCNV.class <- gsub(pattern = 5, replacement = "tumor", x = inferCNV.annotation.malignant$inferCNV.class)
#replace inferCNV.class 6 as "tumor"
inferCNV.annotation.malignant$inferCNV.class <- gsub(pattern = 6, replacement = "tumor", x = inferCNV.annotation.malignant$inferCNV.class)
#update colnames of inferCNV.annotation.malignant 
colnames(inferCNV.annotation.malignant) <- c("cell_id", "Epithelial_cluster", "inferCNV_annotation")

table(inferCNV.annotation.malignant$inferCNV_annotation)

但是很明顯，這個(gè)代碼并不適合我自己走一波流程的數(shù)據(jù)結(jié)果解讀。我們需要一些量化指標(biāo)來(lái)具體判定每個(gè)細(xì)胞亞群是否惡性，比如計(jì)算具體每個(gè)細(xì)胞的CNV score。

后面我們來(lái)繼續(xù)分享進(jìn)階方案！