单细胞转录组高级分析四：scRNA数据推断CNV

上期专题我们介绍了单细胞转录组数据的基础分析，然而那些分析只是揭开了组织异质性的面纱，还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容：多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨!

inferCNV简介

inferCNV是大名鼎鼎的broad研究所开发的，可以使用单细胞转录组数据分析肿瘤细胞CNV。相关文章2014年就发表在了Science上，之后算法不断优化，分析结果也多次刊登在CNS文章上。大家不必担心分析可靠性的问题，反正我是相信品牌的力量?

分析原理

We reasoned that averaging the relative expression of a large number of genomically-adjacent genes would average out gene-specific expression patterns and yield profiles that primarily reflect chromosomal copy number variations (CNVs). We sorted all analyzed genes by their genomic locations, first by chromosome (from chromosome 1 to chromosome 22 and ending with chromosome X), and then by the gene start position. We used a moving average of 100 analyzed genes and the following equation to estimate of chromosomal CNVs in each cell and at each analyzed gene. 大意就是用染色体上相邻的101个基因的平均表达值表示中间位置基因的CNV值，计算公式如下：

这是最核心的计算公式，具体分析比这个公式复杂很多，对原理感兴趣的朋友可以参考文末references列出的文献。

分析流程

流程图所示主要分析过程：

原始数据进行标准化处理，去除测序深度不同造成的差异，结果是流程图中的第1个热图；
依据上文“分析原理”提到的公式把基因表达值转换为CNV值，结果是流程图中的第2个热图；
对每个细胞的CNV值进行中心化处理，使细胞之间相同染色体区域的CNV值具有可比性，结果是流程图中的第3个热图；
用肿瘤细胞的CNV值减去对应区域正常细胞的CNV值，使热图展现的CNV结果更直观，结果是流程图中的第4个热图；
如果infercnv::run函数中的参数denoise=TRUE，则使用算法进一步去除背景噪音凸显CNV区域，结果是流程图中的蓝框左图；
如果infercnv::run函数中的参数HMM=TRUE，则使用隐马尔可夫模型（Hidden Markov Model, HMM）预测CNV区域，并用贝叶斯潜在混合模型（Bayesian Network Latent Mixture Model）对结果进行校正，结果是流程图中的蓝框下图。

安装inferCNV

安装inferCNV之前需要安装JAGS程序，下载地址：

https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/

此程序安装之后，inferCNV依赖的rjags包才能正常安装，否则报错：configuration failed for package ‘rjags’

#安装发行版，作者推荐
if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")
BiocManager::install("infercnv")
#安装github上的最新版
library("devtools")
devtools::install_github("broadinstitute/infercnv")

示例数据演示

inferCNV自带了测试数据，可以先运行看看

library(infercnv)
library(tidyverse)
dir.create("inferCNV")
dir.create("inferCNV/test")
##读取示例数据目录
exprMatrix = system.file("extdata", "oligodendroglioma_expression_downsampled.counts.matrix.gz", package="infercnv")
cellAnnota = system.file("extdata", "oligodendroglioma_annotations_downsampled.txt", package="infercnv")
geneLocate = system.file("extdata", "gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt", package="infercnv")
#创建inferCNV对象，直接给相应的文件路径即可
infercnv_obj = CreateInfercnvObject(delim = 't',
                  raw_counts_matrix = exprMatrix,
                  annotations_file = cellAnnota,
                  gene_order_file = geneLocate,
                  ref_group_names = c("Microglia/Macrophage","Oligodendrocytes (non-malignant)"))

##分析细胞CNV
#cutoff阈值，Smart-seq2数据选“1”, 10x Genomics数据选“0.1”
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1, 
                             out_dir='inferCNV/test', 
                             cluster_by_groups=TRUE, 
                             denoise=TRUE,
                             HMM=TRUE)

10x数据实操

输入文件要求

inferCNV运行需要三个输入文件，原始表达矩阵、细胞注释文件和基因定位文件，具体格式如下：

inferCNV说明要raw_counts_matirx，但它的示例数据是处理后的数据

注意细胞注释文件和基因定位文件都不要表头

GSE149180数据集

inferCNV是针对scRNA数据分析肿瘤细胞CNV开发的，因此我们需要一个包含肿瘤细胞的10x scRNA数据集。这次分析我们使用GSE149180的数据，原文研究对象是小细胞肺癌，只有1个单细胞样本。数据来源文章如下：

根据gse编号可以在GEO数据库下载到cellranger输出文件，原始数据有2000多个细胞，但是数据质量好像不太好，原始数据的相关指标如下：

我使用常规参数过滤之后还有差不多900个细胞，之后降维聚类并用SingleR识别细胞类型，整理好的seurat对象保存为scRNAsclc.rds。

准备表达矩阵和注释文件

library(Seurat)
library(infercnv)
library(tidyverse)
scRNAsclc <- readRDS("inferCNV/scRNAsclc.rds")
cellAnnota <- subset(scRNAsclc@meta.data, select='celltype')
exprMatrix <- as.matrix(GetAssayData(scRNAsclc, slot='counts'))
write.table(exprMatrix, 'inferCNV/exprMatrix.txt', col.names=NA, sep='t')
write.table(cellAnnota, 'inferCNV/cellAnnota.txt', col.names=F, sep='t')

基因定位文件

基因定位文件相对不好准备，有三个途径获得：

如果用自己测序的数据分析，可以找测序公司要，这个最准确;
使用broad准备好的文件，问题是版本太老了，会有一些数据丢失；
下载基因组对应的GTF文件，自己提取基因坐标信息，这个对初学者可能有点困难。

broad准备好的文件：https://data.broadinstitute.org/Trinity/CTAT/cnv/

10x scRNA数据一般比对GRCh38版基因组，下载文件请选择V21版

分析代码

#创建inferCNV对象
#all_celltype: B_cell Monocyte  Neurons  NK_cell  T_cells
infercnv_obj = CreateInfercnvObject(delim = 't',
                  raw_counts_matrix = 'inferCNV/exprMatrix.txt',
                  annotations_file = 'inferCNV/cellAnnota.txt',
                  gene_order_file = 'inferCNV/geneLocate.txt',
                  ref_group_names = c("B_cell","Monocyte","NK_cell","T_cells"))
dir.create("inferCNV/gse149180")
#10x数据cutoff推荐使用0.1
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=0.1, 
                             out_dir='inferCNV/gse149180/', 
                             cluster_by_groups=TRUE, 
                             denoise=TRUE,
                             HMM=TRUE)

分析结果概览

最终会输出很多文件在out_dir的目录下，它们可以分为中间过程数据、初步分析结果、去噪分析结果、HMM预测后结果、最终分析结果5部分。注意文件名称，infercnv开头的是分析结果，其他都是中间数据。剩余的文件注意关键字preliminary、denosied、HMM，分别是初步结果、去噪结果和HMM预测结果，没有这些关键字的是最终结果。我们以最终结果的系列文件来说明一下：

infercnv.png : 去噪之后的最终热图

infercnv.references.txt : 正常细胞的CNV分值矩阵

infercnv.observations.txt : 肿瘤细胞的CNV分值矩阵

infercnv.observation_groupings.txt : 肿瘤细胞的聚类分组

infercnv.observations_dendrogram.txt : 热图的newick格式文件

HMM模型预测的结果不是连续的分值，而是用6种数值代表不同的状态：

0：完全缺失数的变异
0.5：缺失一个拷贝数的变异
1：正常
1.5：增加一个拷贝数的变异
2：增加两个拷贝数的变异
3：所有大于两个拷贝数的变异

References

https://github.com/broadinstitute/inferCNV

1. Anoop P. Patel, Itay Tirosh, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014 Jun 20: 1396-1401

2. Tirosh I et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016 Apr 8;352(6282):189-96

3. Tirosh I et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature. 2016 Nov 10;539(7628):309-313. PubMed PMID: 27806376; PubMed Central PMCID: PMC5465819.

4. Venteicher AS, Tirosh I, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. 2017 Mar 31;355(6332).PubMed PMID: 28360267; PubMed Central PMCID: PMC5519096.

5. Puram SV, Tirosh I, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017 Dec 14;171(7):1611-1624.e24. PubMed PMID: 29198524; PubMed Central PMCID: PMC5878932.

获取帮助

本教程的目的在于把常用的单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以点击 “阅读原文” 找到作者联系方式，获取帮助。

往期回顾