CNS图表复现13—使用inferCNV来区分肿瘤细胞的恶性与否

生信技能树jimmy

发布于 2020-11-02 02:24:45

10.2K10

代码可运行

文章被收录于专栏：单细胞天地单细胞天地

运行总次数：0

代码可运行

分享是一种态度

前言

CNS图表复现之旅前面我们已经进行了12讲，你可以点击图表复现话题回顾。如果你感兴趣也想加入交流群，自己去：你要的rmarkdown文献图表复现全套代码来了（单细胞）找到我们的拉群小助手哈。

我在CNS图表复现09—上皮细胞可以区分为恶性与否，简单演示了，第一次分群后选择上皮细胞后继续分群，第1,2,7,14,21,23,25 是跨越病人的聚类情况，所以初步判定它们这些亚群是正常的上皮细胞。因为目前主流看法是每个病人的恶性肿瘤细胞都是没办法跨越病人进行聚类的。

但最常规做法是使用inferCNV算法可以区分细胞恶性与否。比如online 29 April 2020的文章《Single-Cell Transcriptome Analysis Reveals Intratumoral Heterogeneity in ccRCC, which Results in Different Clinical Outcomes》，就是选取ccRCC的3个病例的21个样本（12个肿瘤，9个对照），质控后总计24550个细胞，使用inferCNV算法可以区分成为7786个非恶性，16764个恶性细胞。

回到我们的这个文章，首先看看文章对的描述：

inferCNV方法描述

文章运行inferCNV算法后的结果图表在附件：

inferCNV算法结果图表

我们在《单细胞天地》多次分享过inferCNV分析的教程：使用inferCNV分析单细胞转录组中拷贝数变异 ，还专门强调了：infercnv输入文件的制作

制作3个文件

下面的代码，基于很多前面步骤产生的文件，所以请务必查看前面的10讲教程哈。

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
load(file = 'first_sce.Rdata')  
load(file = 'phe-of-first-anno.Rdata')
sce=sce.first
table(phe$immune_annotation)
sce@meta.data=phe 
table(phe$immune_annotation,phe$seurat_clusters) 
# BiocManager::install("infercnv")
library(infercnv)

epi.cells  <- row.names(sce@meta.data)[which(phe$immune_annotation=='epi')]
length(epi.cells)
epiMat=as.data.frame(GetAssayData(subset(sce, cells=epi.cells)))

cells.use <- row.names(sce@meta.data)[which(phe$immune_annotation=='stromal')]
length(cells.use)
sce <-subset(sce, cells=cells.use)  
sce 
load(file = 'phe-of-subtypes-stromal.Rdata')
sce@meta.data=phe
DimPlot(sce, reduction = "tsne", group.by = "singleR")
table(phe$singleR)
fib.cells  <-  row.names(sce@meta.data)[phe$singleR=='Fibroblasts']
endo.cells  <-  row.names(sce@meta.data)[phe$singleR=='Endothelial_cells']
fib.cells=sample(fib.cells,800)
endo.cells=sample(endo.cells,800)
fibMat=as.data.frame(GetAssayData(subset(sce, cells=fib.cells)))
endoMat=as.data.frame(GetAssayData(subset(sce, cells=endo.cells)))

dat=cbind(epiMat,fibMat,endoMat)
groupinfo=data.frame(v1=colnames(dat),
                     v2=c(rep('epi',ncol(epiMat)),
                          rep('spike-fib',300),
                          rep('ref-fib',500),
                          rep('spike-endo',300),
                          rep('ref-endo',500)))

library(AnnoProbe)
geneInfor=annoGene(rownames(dat),"SYMBOL",'human')
colnames(geneInfor)
geneInfor=geneInfor[with(geneInfor, order(chr, start)),c(1,4:6)]
geneInfor=geneInfor[!duplicated(geneInfor[,1]),]
length(unique(geneInfor[,1]))
head(geneInfor)
## 这里可以去除性染色体
# 也可以把染色体排序方式改变
dat=dat[rownames(dat) %in% geneInfor[,1],]
dat=dat[match( geneInfor[,1], rownames(dat) ),] 
dim(dat)
expFile='expFile.txt'
write.table(dat,file = expFile,sep = '\t',quote = F)
groupFiles='groupFiles.txt'
head(groupinfo)
write.table(groupinfo,file = groupFiles,sep = '\t',quote = F,col.names = F,row.names = F)
head(geneInfor)
geneFile='geneFile.txt'
write.table(geneInfor,file = geneFile,sep = '\t',quote = F,col.names = F,row.names = F)

table(groupinfo[,2])

 
       epi   ref-endo    ref-fib spike-endo  spike-fib 
      5444        500        500        300        300

代码挺简单的，如果你的R语言还不错，很容易看懂，就是之前我们把单细胞已经分好群了，这个时候取全部的上皮细胞，以及部分Fibroblasts和Endothelial_cells细胞来一起运行inferCNV流程。

这个时候Fibroblasts和Endothelial_cells细胞是正常的二倍体细胞，而全部的上皮细胞里面就会根据CNV情况来区分成为恶性与否的肿瘤细胞。

有一个包（annoprobe）需要安装

上面的代码里面我的方法，采用了annoprobe获取基因坐标，安装annoprobe的代码超级简单：

library(remotes)
# https://git-scm.com/downloads 
# 居然有些人电脑里面没有git软件，可怕!!!
url='https://gitee.com/jmzeng/annoprobe.git'
install_git(url)

当然，你也可以自行前往下载gtf文件或者其它方法获取基因的坐标信息，自己制作 geneInfor 变量后写入基因坐标文件。

安装annoprobe的时候不要更新如何其它R包哈：

基本上都可以安装成功：

但是我看有学生反馈，安装的时候提示他的电脑缺乏git软件，我就很纳闷了，一个搞生物信息学数据分析的电脑里面，怎么可能没有git软件呢？就算是没有，马上去下载git安装也可以啊！！！

运行infercnv

因为前面已经做好了3个文件，所以下面的代码无需改变，直接运行即可：

infercnv_obj = CreateInfercnvObject(raw_counts_matrix=expFile,
                                    annotations_file=groupFiles,
                                    delim="\t",
                                    gene_order_file= geneFile,
                                    ref_group_names=c("WT"))  ## 这个取决于自己的分组信息里面的

infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
                             out_dir=tempfile(), 
                             cluster_by_groups=TRUE, 
                             denoise=TRUE,
                             HMM=TRUE)

运行起来还是蛮快的，20min搞定全部分析：

INFO [2020-10-15 17:31:04] ::process_data:Start
INFO [2020-10-15 17:31:04] Creating output path /var/folders/2x/l4fzfvy17gz94j1sz6klhyjh0000gn/T//Rtmph8chdX/file4bcc7985a2a0
INFO [2020-10-15 17:31:04] Checking for saved results.
INFO [2020-10-15 17:31:04] 

 STEP 1: incoming data
 # 此处省略----
 STEP 21: Denoising
INFO [2020-10-15 17:53:44] Quantiles of plotted data range: 0.598977210060562,1.00469525380699,1.00469525380699,1.00469525380699,1.40102278993944
INFO [2020-10-15 17:53:46] plot_cnv_observations:Writing observation data to /var/folders/2x/l4fzfvy17gz94j1sz6klhyjh0000gn/T//Rtmph8chdX/file4bcc7985a2a0/infercnv.observations.txt
INFO [2020-10-15 17:53:48] plot_cnv_references:Start
INFO [2020-10-15 17:53:48] Reference data size: Cells= 1000 Genes= 366
INFO [2020-10-15 17:53:48] plot_cnv_references:Number reference groups= 2
INFO [2020-10-15 17:53:48] plot_cnv_references:Plotting heatmap.
INFO [2020-10-15 17:53:48] Colors for breaks:  #00008B,#24249B,#4848AB,#6D6DBC,#9191CC,#B6B6DD,#DADAEE,#FFFFFF,#EEDADA,#DDB6B6,#CC9191,#BC6D6D,#AB4848,#9B2424,#8B0000
INFO [2020-10-15 17:53:48] Quantiles of plotted data range: 0.598977210060562,1.00469525380699,1.00469525380699,1.00469525380699,1.40102278993944
INFO [2020-10-15 17:53:48] plot_cnv_references:Writing reference data to /var/folders/2x/l4fzfvy17gz94j1sz6klhyjh0000gn/T//Rtmph8chdX/file4bcc7985a2a0/infercnv.references.txt

因为前面我没有修改输出路径，仍然是在 tempfile() ，所以需要自己进入这个临时文件夹，把全部的分析结果拷贝出来。

有意思的是，我的结果有点诡异，明明是作为二倍体正常细胞参考集的Fibroblasts和Endothelial_cells细胞居然也是在某些染色体上面有明显的CNV情况。

如下所示：

infercnv

不知道为什么我的ref居然不是很干净的样子？后面我们来慢慢探讨第一次运行inferCNV得到这个诡异的结果的原因。

10X单细胞转录组数据是否可以inferCNV

这个文章的单细胞转录组数据来源于smart-seq2技术，但是现在单细胞转录组领域其实是10X数据的天下，那么10X转录组数据是否可以inferCNV呢？早在2018年就有文章，是 Comprehensive analysis of immune evasion in breast cancer by single-cell RNA-seq ，链接是. doi: http://dx.doi.org/10.1101/368605 bioRxiv preprint first posted online Jul. 13, 2018; 就是使用10X转录组数据来推断CNV信息。因为**10X技术出来的单个细胞的reads数量太少，检测到的基因数量太少，**但是inferCNV更新后有一个参数是cutoff ，就很清楚的指出来了：cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics 说明它应用于10X单细胞转录组数据是没有问题的。

10X单细胞转录组数据处理流程在我们单细胞天地有详细介绍：

如果你查看10X单细胞转录组的cellranger报告，会发现显示平均每个细胞的测序数据量是45K条reads。当然，并不是10x一个技术是这样单个细胞的reads数量太少，检测到的基因数量太少。比如文章：Li et al., Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma, Cell (2019), https://doi.org/10.1016/j.cell.2018.11.043 ：同样的，平均每个细胞也就40K左右的reads数量啦。