引言
本系列)将开展全新的转录组分析专栏,主要针对使用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 的四种方法:
- 从转录丰度文件和 tximport 生成
- 从计数矩阵生成
- 从 htseq-count 文件生成
- 从 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。