RNA-seq 差异分析的点点滴滴（1）

引言

本系列^[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。