RNA-seq 差异分析的点点滴滴(1)

简介: RNA-seq 差异分析的点点滴滴(1)

引言

本系列)将开展全新的转录组分析专栏,主要针对使用DESeq2时可能出现的问题和方法进行展开。

为何使用未经标准化的计数数据?

DESeq2 工具包在接收输入时,期望得到的是未经处理的原始计数数据,比如从 RNA-seq 或其他高通量测序实验中获得的,这些数据以整数值矩阵的形式呈现。在这个矩阵中,第 i 行第 j 列的数值表示在样本 j 中可以归属于基因 i 的读段数。同样地,对于其他类型的实验,矩阵的行可能代表结合区域(例如 ChIP-Seq 实验)或肽序列(例如定量质谱实验)。

矩阵中的数值应当是未经标准化的读段计数(对于单端 RNA-seq)或片段计数(对于双端 RNA-seq)。RNA-seq 的工作流程中描述了多种制备此类计数矩阵的技术。为 DESeq2 的统计模型提供计数矩阵作为输入非常关键,因为只有原始的计数数据才能准确评估测量的精确度。DESeq2 模型在内部会校正文库大小的影响,因此不应该使用经过转换或标准化的数值,比如按文库大小调整后的计数,作为输入数据。

DESeqDataSet 对象

DESeq2 工具包中,用于存储读取计数和统计分析过程中的中间估计量的类对象是 DESeqDataSet,通常在代码中以 dds 表示。

技术细节上,DESeqDataSet 类扩展了 SummarizedExperiment 包中的 RangedSummarizedExperiment 类。“Ranged” 指的是测定数据的行(即计数)可以与基因组的特定区域(如基因的外显子)相关联。

DESeqDataSet 对象必须关联一个设计公式。这个公式描述了将在模型中使用的变量,通常以波浪号 (~) 开始,后跟用加号 (+) 分隔的变量(如果不是公式形式,系统会自动转换)。设计公式可以在后续更改,但需要重新执行所有差异分析步骤,因为设计公式用于估计离散度和模型的 log2 倍数变化。

注意:为了利用包的默认设置,应将感兴趣的变量放在公式的末尾,并确保对照组水平是第一水平。

接下来,将展示根据在 DESeq2 之前使用的管道不同,构建 DESeqDataSet 的四种方法:

  1. 从转录丰度文件和 tximport 生成
  2. 从计数矩阵生成
  3. 从 htseq-count 文件生成
  4. 从 SummarizedExperiment 对象生成

转录本丰度数据

建议在使用 DESeq2 之前,先采用快速的转录本丰度定量工具,然后通过 tximport导入这些定量数据来创建 DESeq2 所需的基因水平计数矩阵。这种方法允许用户从多种外部软件中导入转录本丰度估计值,包括以下方法:Salmon; Sailfish; kallisto ;RSEM

采用上述方法进行转录本丰度估计的好处包括:(i)这种方法能够校正样本间可能的基因长度变化(例如,由于异构体的不同使用),(ii)其中一些方法(Salmon, Sailfish, kallisto)相比需要创建和存储 BAM 文件的基于比对的方法,速度显著更快,且对内存和磁盘空间的需求更少,以及(iii)可以避免丢弃那些能够与多个具有同源序列的基因对齐的片段,从而提高检测的灵敏度。

请注意,tximport-to-DESeq2 方法使用的是转录本丰度定量器估计的基因计数,而不是标准化计数。

在这里,将展示如何从存储在 tximportData 包中的 Salmon quant.sf 文件导入转录本丰度,并构建一个基因水平的 DESeqDataSet 对象。

library("tximport")
library("readr")
library("tximportData")
dir <- system.file("extdata", package="tximportData")
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
samples$condition <- factor(rep(c("A","B"),each=3))
rownames(samples) <- samples$run
samples[,c("pop","center","run","condition")]

##           pop center       run condition
## ERR188297 TSI  UNIGE ERR188297         A
## ERR188088 TSI  UNIGE ERR188088         A
## ERR188329 TSI  UNIGE ERR188329         A
## ERR188288 TSI  UNIGE ERR188288         B
## ERR188021 TSI  UNIGE ERR188021         B
## ERR188356 TSI  UNIGE ERR188356         B

接下来,使用适当的样本列指定文件的路径,并读取一个将转录本与该数据集的基因链接起来的表。

files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))

使用 tximport 函数导入 DESeq2 所需的量化数据。

txi <- tximport(files, type="salmon", tx2gene=tx2gene)

最后,可以根据样本中的 txi 对象和样本信息构造一个 DESeqDataSet。

library("DESeq2")
ddsTxi <- DESeqDataSetFromTximport(txi,
                                   colData = samples,
                                   design = ~ condition)

这里的ddsTxi对象就可以在下面的分析步骤中用作dds。

相关文章
|
7天前
|
Python
RNA-seq 差异分析的点点滴滴(2)
RNA-seq 差异分析的点点滴滴(2)
26 10
RNA-seq 差异分析的点点滴滴(2)
|
6天前
|
SQL 数据挖掘 Python
R中单细胞RNA-seq数据分析教程 (1)
R中单细胞RNA-seq数据分析教程 (1)
23 5
R中单细胞RNA-seq数据分析教程 (1)
|
5月前
|
数据可视化 Java 数据处理
单细胞|RNA-seq & ATAC-seq 联合分析
单细胞|RNA-seq & ATAC-seq 联合分析
62 3
|
6月前
|
机器学习/深度学习 SQL 数据可视化
单细胞分析(Signac): PBMC scATAC-seq 整合
单细胞分析(Signac): PBMC scATAC-seq 整合
73 0
|
6月前
|
数据可视化 数据挖掘 Serverless
单细胞分析(Signac): PBMC scATAC-seq 聚类
单细胞分析(Signac): PBMC scATAC-seq 聚类
52 0
|
6月前
|
存储 移动开发 Shell
单细胞分析(Signac): PBMC scATAC-seq 预处理
单细胞分析(Signac): PBMC scATAC-seq 预处理
81 2
|
6月前
|
数据可视化 数据挖掘 索引
R语言层次聚类、多维缩放MDS分类RNA测序(RNA-seq)乳腺发育基因数据可视化|附数据代码2
R语言层次聚类、多维缩放MDS分类RNA测序(RNA-seq)乳腺发育基因数据可视化|附数据代码
|
6月前
|
存储 数据可视化 数据挖掘
R语言层次聚类、多维缩放MDS分类RNA测序(RNA-seq)乳腺发育基因数据可视化|附数据代码1
R语言层次聚类、多维缩放MDS分类RNA测序(RNA-seq)乳腺发育基因数据可视化|附数据代码
|
6月前
|
Java 数据挖掘 Go
JCR一区7.7分|单细胞联合bulk-seq的线粒体自噬,分析方法都挺好
这篇文章介绍了研究者通过分析单细胞和Bulk RNA测序数据,鉴定出18个与胃癌(GC)进展相关的线粒体自噬相关基因(MRG),并建立了基于这些基因的预后模型。研究发现GABARAPL2和CDC37可能是GC的预后标志物和潜在治疗靶点。此外,分析揭示了细胞间通讯模式和免疫浸润状态,暗示MRG可能影响GC的免疫治疗响应。整体而言,这项工作为GC的诊断和治疗提供了新见解。
120 0
|
数据可视化 关系型数据库 数据挖掘
scRNA分析|一(尽)文(力)解决你的单细胞火山图问题
scRNA分析|一(尽)文(力)解决你的单细胞火山图问题
671 0