單細胞分析實錄(12): 如何推斷腫瘤細胞

公眾號后台回復20210418，獲取本文的測試數據

1. 回顧

前面我已經講解了inferCNV的基本用法，這一節繼續學習：如何推斷腫瘤細胞。本文應該是到目前為止，全網最詳細的帖子（之一）。

關於腫瘤細胞的推斷，比如研究上皮細胞起源的腫瘤，我見過的一些做法如下

1. 簡單粗暴型：將來源腫瘤病灶的所有上皮細胞都看成是腫瘤細胞
2. 看CNV：將來源於腫瘤病灶的上皮細胞和參考細胞比較，將那些有明顯CNV的細胞看作腫瘤細胞
3. 看表達特征：考慮到有的腫瘤細胞可能沒有CNV層面的變化，有的文章會根據腫瘤與正常的表達譜來區分是否為腫瘤細胞

如果你假設即使來源於病灶，也可能有正常（惡行程度很低）的上皮細胞，那么第2,3種做法就顯得合理一些。實際數據處理中，我也確實看到過這種情況。這一節內容，我主要介紹第2種思路，第3種思路簡單說一下。

2. 看CNV

既然是看，那就要選擇合適的圖，來輔助你看。inferCNV跑完的結果會顯示一張熱圖

就比如這張圖，相比於上方的參考panel，下方的觀測panel可見一些明顯的CNV，當然也可見一些細胞並沒有明顯的CNV，這些就是惡性程度很低的細胞。現在需要做的就是將這兩種類型區分出來。

2.1 方法一

將參考和觀測panel對應的矩陣提取出來，混在一起，聚類，之后每一個細胞會有一個類別label，挑出是腫瘤細胞的類別即可。為什么要混在一起呢？其實也是提供一個參考。聚成多少類，不需要特別精確，多於你肉眼能看出的類別數就行。不宜過多，操作起來麻煩；不宜過少，不然沒有區分度。

示例數據之前已經跑好了，代碼在try1.R，inferCNV的結果保存在try1文件夾中，利用的是run.final.infercnv_obj文件，里面保存了最終畫熱圖用到的表達矩陣，以及細胞注釋。下面是代碼，有點長，先上圖

library(infercnv)
library(tidyverse)
library(ComplexHeatmap)
library(circlize)
library("RColorBrewer")

infercnv_obj = readRDS("./try1/run.final.infercnv_obj")
expr <- [e-mail]
normal_loc <- infercnv_obj@reference_grouped_cell_indices
normal_loc <- normal_loc$normal
test_loc <- infercnv_obj@observation_grouped_cell_indices
test_loc <- test_loc$test

anno.df=data.frame(
  CB=c(colnames(expr)[normal_loc],colnames(expr)[test_loc]),
  class=c(rep("normal",length(normal_loc)),rep("test",length(test_loc)))
)
head(anno.df)

gn <- rownames(expr)
geneFile <- read.table("EXAMPLE_ONLY_DONT_REUSE.txt",header = F,sep = "\t",stringsAsFactors = F)
rownames(geneFile)=geneFile$V1
sub_geneFile <-  geneFile[intersect(gn,geneFile$V1),]
expr=expr[intersect(gn,geneFile$V1),]
head(sub_geneFile,4)
expr[1:4,1:4]

#聚類，7類，提取結果
set.seed(20210418)
kmeans.result <- kmeans(t(expr), 7)
kmeans_df <- data.frame(kmeans_class=kmeans.result$cluster)
kmeans_df$CB=rownames(kmeans_df)
kmeans_df=kmeans_df%>%inner_join(anno.df,by="CB") #合並
kmeans_df_s=arrange(kmeans_df,kmeans_class) #排序
rownames(kmeans_df_s)=kmeans_df_s$CB
kmeans_df_s$CB=NULL
kmeans_df_s$kmeans_class=as.factor(kmeans_df_s$kmeans_class) #將kmeans_class轉換為因子，作為熱圖的一個注釋，最終在熱圖上就會按照1:7排序
head(kmeans_df_s)

#定義熱圖的注釋，及配色
top_anno <- HeatmapAnnotation(foo = anno_block(gp = gpar(fill = "NA",col="NA"), labels = 1:22,labels_gp = gpar(cex = 1.5)))
color_v=RColorBrewer::brewer.pal(8, "Dark2")[1:7] #類別數
names(color_v)=as.character(1:7)
left_anno <- rowAnnotation(df = kmeans_df_s,col=list(class=c("test"="red","normal" = "blue"),kmeans_class=color_v))

#下面是繪圖
pdf("try1.pdf",width = 15,height = 10)
ht = Heatmap(t(expr)[rownames(kmeans_df_s),], #繪圖數據的CB順序和注釋CB順序保持一致
                col = colorRamp2(c(0.4,1,1.6), c("#377EB8","#F0F0F0","#E41A1C")), #如果是10x的數據，這里的刻度會有所變化
                cluster_rows = F,cluster_columns = F,show_column_names = F,show_row_names = F,
                column_split = factor(sub_geneFile$V2, paste("chr",1:22,sep = "")), #這一步可以控制染色體順序，即使你的基因排序文件順序是錯的
                column_gap = unit(2, "mm"),
                
                heatmap_legend_param = list(title = "Modified expression",direction = "vertical",title_position = "leftcenter-rot",at=c(0.4,1,1.6),legend_height = unit(3, "cm")),
                
                top_annotation = top_anno,left_annotation = left_anno, #添加注釋
                row_title = NULL,column_title = NULL)
draw(ht, heatmap_legend_side = "right")
dev.off()

#每一類對應的CB保存在kmeans_df_s數據框中
write.table(kmeans_df_s, file = "kmeans_df_s.txt", quote = FALSE, sep = '\t', row.names = T, col.names = T)

這張圖中，第1, 5類的細胞包含了normal以及少量的test，從圖中看基本也是沒啥CNV，所以將1，5剔除，剩下的5類定義為腫瘤細胞，對應的CB可以在kmeans_df_s.txt中獲取。

2.2 方法二

除了上面這種純看圖的方法，有的文獻還會對每個細胞的CNV水平做一個定量，我這里參考：Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in pancreatic ductal adenocarcinoma. 原文的CNV score是這么計算的"gene expression of cells was re-standardized and values were limited as −1 to 1. The CNV score of each cell was calculated as quadratic sum of CNV region"。

原文獻的圖如下：

我是這么理解的，將inferCNV的返回值做一些調整，使其背景值為0，然后再求每個細胞的所有基因CNV值的平方和。這里我覺得對inferCNV返回值矯正是不必要的（這些值已經是矯正過的），每個值減1之后求平方和就可以了。另外求和還不夠，得除以基因數，因為每個矩陣inferCNV選用的基因數是不一樣的，所以其實就是看方差，即每個細胞所有基因偏離背景的程度。

我仍然以方法一中的7個類為例，計算它們的CNV level。

expr2=expr-1
expr2=expr2 ^ 2
CNV_score=as.data.frame(colMeans(expr2))
colnames(CNV_score)="CNV_score"
CNV_score$CB=rownames(CNV_score)
kmeans_df_s$CB=rownames(kmeans_df_s)
CNV_score=CNV_score%>%inner_join(kmeans_df_s,by="CB")

CNV_score%>%ggplot(aes(kmeans_class,CNV_score))+geom_violin(aes(fill=kmeans_class),color="NA")+
  scale_fill_manual(values = color_v)+
  theme_bw()
ggsave("CNV level.pdf",width = 10,height = 6,units = "cm")

可以看到，仍然是1，5兩類值最低，且在0附近，加上這兩類大部分是normal cell, 所以將他們划分為非惡性應該是沒有問題的。

網上有幾篇關於計算CNV score的博客，如下，是將inferCNV的矩陣轉換為0 1 2三種CNV程度，再在每個細胞中求和，這種思路和上面是一樣的，但是我覺得實際操作起來不好做，主要是什么范圍被划分為1還是2，不好說，不同的平台不同的數據集，數據范圍有差異，比如我的10X數據最大1.3最小0.7，圖片中的定義明顯就不對。另外，這種方法在文獻中沒怎么看到（知道的讀者可以回復一下）。

3. 看表達特征

這種方法，我看到的不多，主要參考2020年的一篇文章：Dissecting transcriptional heterogeneity in primary gastric adenocarcinoma by single cell RNA sequencing。

在這篇胃癌的單細胞文獻中，作者首先根據TCGA中的癌症和正常配對樣本做差異分析，得到了兩個初始signature，一個癌症樣本高表達，一個正常樣本高表達，各50個基因。

• 用兩個signature對每一個單細胞進行評分，可以得到一個cell_number x 2的矩陣
• 接下來用k-means對所有的單細胞做聚類，只聚成兩類，高表達癌症signature的定義為候選惡性細胞，高表達正常signature的定義為候選非惡性細胞
• 繼續對候選惡性、非惡性細胞做差異分析，得到兩個新的signature，一類表達惡性，一類表示非惡性

重復以上步驟，直到最終的兩個signature的基因穩定。此處穩定的理解，我之前一直覺得是這一次的signature和上一次的signature相同，跑過幾次代碼之后，覺得這樣理解太嚴格，絕大部分相同應該就可以。下面是原文獻的圖：

這種方法之前用過，但對於我的數據集效果不是很好，反映在圖上就是最終得不到文獻圖那么清晰的分界(見上圖)，一堆是惡性，一堆是非惡性，感興趣的讀者可以用這種方法試試自己的數據。下面的演示，我還是使用同上文一樣的數據集，只不過從TCGA找初始signature這一步我省了，使用的是兩個已經可以較好區分惡性、非惡性細胞的signature（比TCGA signature更接近最終signature）。
主體代碼如下，我只展示對於初始signature的操作，后續的迭代過程，代碼類似，我就不展示了，有需要的朋友后台私信小編。

library(Seurat)
library(tidyverse)

testdf=read.table("count_matrix.txt",header = T,row.names = 1,sep = "\t")
test.seu=CreateSeuratObject(counts = testdf)
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
[e-mail]$CB=rownames([e-mail])
[e-mail]=[e-mail] %>% inner_join(anno.df,by="CB") #anno.df有兩列，一列CB，一列class表示細胞來源
rownames([e-mail])=[e-mail]$CB

TCGA=read.table("signature0.txt",header = T,sep = "\t",stringsAsFactors = F) #通過TCGA得到的兩個signature
colnames(TCGA)=c("pathways","gene")
for (i in unique(TCGA$pathways)) {
  TCGA_small=TCGA%>%filter(pathways==i)
  genes.for.scoring <- list(TCGA_small$gene)
  test.seu <- AddModuleScore(object = test.seu, features = genes.for.scoring, name = i) #AddModuleScore計算得分
}

metadf=[e-mail][,c("normal1","tumor1","class","CB")] #初始得分
metadf%>%ggplot(aes(x=normal1,y=tumor1,color=class))+geom_point()
ggsave("try_0.a.png")
#在示例數據class列中，test表示腫瘤病灶取樣，normal表示正常對照；在signature中，normal表示正常樣本，tumor表示癌症樣本

#kmeans聚類
kmeans.result <- kmeans(metadf[,1:2],2)
kmeans_df <- data.frame(
  kmeans_class=kmeans.result$cluster
)
kmeans_df$CB=rownames(kmeans_df)
metadf=merge(metadf,kmeans_df,by="CB")

#kmeans為1 2分別代表什么細胞
small1_metadf=metadf%>%filter(kmeans_class==1)
small2_metadf=metadf%>%filter(kmeans_class==2)
if( (mean(small1_metadf$tumor1) < mean(small2_metadf$tumor1)) & (mean(small1_metadf$normal1)>mean(small2_metadf$normal1)) ) {
  print("1:normal; 2:tumor")
  metadf$classification=ifelse(metadf$kmeans_class==1,"normal","tumor")
}else if (mean(small1_metadf$tumor1) > mean(small2_metadf$tumor1) & mean(small1_metadf$normal1) < mean(small2_metadf$normal1)){
  print("1:tumor; 2:normal")
  metadf$classification=ifelse(metadf$kmeans_class==1,"tumor","normal")
}else{
  print("error")
}
metadf%>%ggplot(aes(x=normal1,y=tumor1,color=classification))+geom_point()
ggsave("try_0.b.png")

metadf=metadf[,c("CB","classification")]
colnames(metadf)[2]="classification_0"
[e-mail]=[e-mail]%>%inner_join(metadf,by = "CB")
rownames([e-mail])=[e-mail]$CB
###以上是針對初始signature的操作

這里展示的圖片說明這兩次的結果很接近，且圖形上看也比較好，所以沒必要再繼續迭代；反之，如果繼續，可能會出現兩堆細胞越來越混的情況（一步錯，步步錯）。就像下面這樣，迭代20次之后，基本分不開了：

---

這一篇的內容就到這里，下一期會寫inferCNV在腫瘤細胞克隆進化方面的應用。

因水平有限，有錯誤的地方，歡迎批評指正！