但最常规做法是使用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数量啦。

往期回顾

细胞身份何以在分裂中得以保持？

CNS图表复现12—检查原文的细胞亚群的标记基因

狼来了！聊个天就能做生信分析的人工智能是否要替代一大波生信人员？

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