单细胞分析实录(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 <- infercnv_obj@expr.data
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 = "	",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 = '	', 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 = "	")
test.seu=CreateSeuratObject(counts = testdf)
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
test.seu@meta.data$CB=rownames(test.seu@meta.data)
test.seu@meta.data=test.seu@meta.data %>% inner_join(anno.df,by="CB") #anno.df有两列，一列CB，一列class表示细胞来源
rownames(test.seu@meta.data)=test.seu@meta.data$CB

TCGA=read.table("signature0.txt",header = T,sep = "	",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=test.seu@meta.data[,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"
test.seu@meta.data=test.seu@meta.data%>%inner_join(metadf,by = "CB")
rownames(test.seu@meta.data)=test.seu@meta.data$CB
###以上是针对初始signature的操作

这里展示的图片说明这两次的结果很接近，且图形上看也比较好，所以没必要再继续迭代；反之，如果继续，可能会出现两堆细胞越来越混的情况（一步错，步步错）。就像下面这样，迭代20次之后，基本分不开了：

---

这一篇的内容就到这里，下一期会写inferCNV在肿瘤细胞克隆进化方面的应用。

因水平有限，有错误的地方，欢迎批评指正！