上期专题我们介绍了单细胞转录组数据的基础分析,然而那些分析只是揭开了组织异质性的面纱,还有更多的生命奥秘隐藏在数据中等待我们发掘。本专题将介绍一些单细胞转录组的高级分析内容:多样本批次校正、转录因子分析、细胞通讯分析、基因集变异分析和更全面的基因集富集分析。不足之处请大家批评指正,欢迎添加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列出的文献。
分析流程
流程图所示主要分析过程:
安装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准备好的文件: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种数值代表不同的状态:
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语言基础的朋友参考。分析原理和代码我没有详细解释,官网有详细的教程和权威的解释,翻译好的中文教程也有多个版本,有兴趣的可以搜索一下。如果您学习本教程有一定困难,可以点击 “阅读原文” 找到作者联系方式,获取帮助。
往期回顾