专栏首页单细胞天地单细胞转录组高级分析四:scRNA数据推断CNV

单细胞转录组高级分析四: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. 原始数据进行标准化处理,去除测序深度不同造成的差异,结果是流程图中的第1个热图;
  2. 依据上文“分析原理”提到的公式把基因表达值转换为CNV值,结果是流程图中的第2个热图;
  3. 对每个细胞的CNV值进行中心化处理,使细胞之间相同染色体区域的CNV值具有可比性,结果是流程图中的第3个热图;
  4. 用肿瘤细胞的CNV值减去对应区域正常细胞的CNV值,使热图展现的CNV结果更直观,结果是流程图中的第4个热图;
  5. 如果infercnv::run函数中的参数denoise=TRUE,则使用算法进一步去除背景噪音凸显CNV区域,结果是流程图中的蓝框左图;
  6. 如果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')

基因定位文件

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

  1. 如果用自己测序的数据分析,可以找测序公司要,这个最准确;
  2. 使用broad准备好的文件,问题是版本太老了,会有一些数据丢失;
  3. 下载基因组对应的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语言基础的朋友参考。分析原理和代码我没有详细解释,官网有详细的教程和权威的解释,翻译好的中文教程也有多个版本,有兴趣的可以搜索一下。如果您学习本教程有一定困难,可以点击 “阅读原文” 找到作者联系方式,获取帮助。

往期回顾

单细胞转录组基础分析一:分析环境搭建

单细胞转录组基础分析二:数据质控与标准化

单细胞转录组基础分析三:降维与聚类

单细胞转录组基础分析四:细胞类型鉴定

单细胞转录组基础分析五:细胞再聚类

单细胞转录组基础分析六:伪时间分析

单细胞转录组基础分析七:差异基因富集分析

单细胞转录组基础分析八:可视化工具总结

欢迎加入生信技能树小圈子

期待单细胞工具的大浪淘沙,洗尽铅华




本文分享自微信公众号 - 单细胞天地(sc-ngs)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-09-01

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 使用inferCNV分析单细胞转录组中拷贝数变异

    inferCNV用与探索肿瘤单细胞RNA-seq数据,分析其中的体细胞大规模染色体拷贝数变化(copy number alterations, CNA), 例...

    生信技能树jimmy
  • 单细胞全基因组测序—直肠癌的异质性

    当你的才华还撑不起你的野心时,请潜下心来,脚踏实地,跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了,通过文献速递这个栏目很幸运聚集了一些小伙伴...

    生信技能树jimmy
  • 桑基图在单细胞数据探索中的应用

    桑基图(Sankey diagram),即桑基能量分流图,也叫桑基能量平衡图。它是一种特定类型的流程图,图中延伸的分支的宽度对应数据流量的大小,比较适用于用户流...

    生信技能树jimmy
  • 云计算第一股UCloud:生死博弈刚刚开始

    这些年,优刻得营收规模上去了,但成本的增长更甚一步,使其难以实现盈利,甚至出现大幅度的亏损。而在优刻得上市之后最新发布的首份半年报里,亏损乌云同样笼罩。“云计算...

    刘旷
  • UCloud上市后的新苦海:净利一路下坡

    IDC发布的报告显示,受疫情影响,云计算技术将从互联网走向非互联网,深刻影响整个社会的生产生活方式,并指出2018年至2022年我国公有云市场复合增长率达 39...

    刘旷
  • [782]AttributeError: module 'tabula' has no attribute 'read_pdf'

    报错:ImportError: cannot import name ‘read_pdf’

    周小董
  • Android 图片选择到裁剪之步步深坑

    前言 最近在自己的项目里实现了一个头像选择的功能,就是先从相册里选取一张图片再调用系统的裁剪功能来制作头像,效果就像下面这样: ? 本以为很小的一个功能,却远远...

    非著名程序员
  • 不到 3 天截稿!NeurIPS 2020 新要求提交的“影响陈述”还不会写怎么办?

    今年 2 月份,NeurIPS 组委会发布了NeurIPS 2020 在提交和评审机制上做出的一些重大更改,其中一项便是要求作者在投稿论文中单独拟一个“影响陈述...

    AI科技评论
  • 为什么计算机最小的存储单位是字节?而最小到的传输单位是bit?

    数据存储是以“字节”(Byte)为单位,数据传输是以大多是以“位”(bit,又名“比特”)为单位,一个位就代表一个0或1(即二进制),每8个位(bit,简写为b...

    wuweixiang
  • 数据传输 | mysqldiff/mysqldbcompare 实现 DTLE 自动化测试

    本文来源:原创投稿 *爱可生开源社区出品,原创内容未经授权不得随意使用,转载请联系小编并注明来源。

    爱可生开源社区

扫码关注云+社区

领取腾讯云代金券