跟着存档教程动手学RNAseq分析（三）：使用DESeq2进行计数标准化
跟着存档教程动手学RNAseq分析（三）：使用DESeq2进行计数标准化
跟着存档教程动手学RNAseq分析（一）
标准化
差异表达分析工作流程的第一步是计数标准化，这是对样本间基因表达进行准确比较所必需的。
img
除了许多不关心因素之外，每个基因的比对reads计数与RNA的表达成正比。标准化是对原始计数值进行缩放以解释无关因素的过程。通过这种方式，表达水平在样本之间和/或样本内部更具有可比性。
标准化过程中经常考虑的主要因素有：
- 测序深度：样品间基因表达的比较需要考虑测序深度。在下面的例子中，每个基因在样本A中的表达似乎比样本B的表达增加了一倍，然而这是样本A测序深度增加一倍的结果。
img
注意:在上图中，每个粉色和绿色的矩形代表一个与基因对齐的read。由虚线连接连接跨越内含子的read。
- 基因长度：为了比较同一样本内不同基因的表达，需要考虑基因长度。在这个例子中，基因X和基因Y有相似的表达水平，但是映射到基因X的reads数会比映射到基因Y的reads数多得多，因为基因X更长。
img
- RNA组成：样品之间有一些高度差异的基因表达，样品之间表达的基因数量的差异，或者存在污染，这些都可能导致某些归一化方法的偏差。为了准确比较样品之间的表达，建议考虑RNA组成，这在进行差异表达分析时尤为重要。在本例中，假设样本A和样本B的测序深度相似，除DE基因外，各基因在样本间的表达水平相似。由于DE基因占了样本B的绝大部分，样本B的计数会有很大的偏倚。因此，B样本中的其他基因似乎比A样本中的相同基因表达更少。
img
虽然标准化对于差异表达分析是必要的，但对于探索性数据分析、数据可视化以及在样本之间或样本内部研究或比较计数时也是必要的。
常见标准化方法
存在几种常见的标准化方法来解释这些差异：
标准化方法 | 描述 | 考虑因素 | 使用推荐 |
|---|---|---|---|
CPM (counts per million) | counts scaled by total number of reads | sequencing depth | gene count comparisons between replicates of the same sample group; NOT for within sample comparisons or DE analysis |
TPM (transcripts per kilobase million) | counts per length of transcript (kb) per million reads mapped | sequencing depth and gene length | gene count comparisons within a sample or between samples of the same sample group; NOT for DE analysis |
RPKM/FPKM (reads/fragments per kilobase of exon per million reads/fragments mapped) | similar to TPM | sequencing depth and gene length | gene count comparisons between genes within a sample; NOT for between sample comparisons or DE analysis |
DESeq2’s median of ratios | counts divided by sample-specific size factors determined by median ratio of gene counts relative to geometric mean per gene | sequencing depth and RNA composition | gene count comparisons between samples and for DE analysis; NOT for within sample comparisons |
EdgeR’s trimmed mean of M values (TMM) | uses a weighted trimmed mean of the log expression ratios between samples | sequencing depth, RNA composition | gene count comparisons between samples and for DE analysis; NOT for within sample comparisons |
RPKM/FPKM (不推荐用于样本间比较）
虽然TPM和RPKM/FPKM归一化方法都考虑了测序深度和基因长度，但RPKM/FPKM不被推荐。原因是RPKM/FPKM方法输出的归一化计数值在样本之间没有可比性。
使用RPKM/FPKM归一化，每个样本的RPKM/FPKM归一化计数总数将不同。因此，你不能在样本之间对每个基因的标准化计数进行相等的比较。
RPKM-normalized计数表
gene | sampleA | sampleB |
|---|---|---|
XCR1 | 5.5 | 5.5 |
WASHC1 | 73.4 | 21.8 |
… | … | … |
Total RPKM-normalized counts | 1,000,000 | 1,500,000 |
例如，在上表中，即使RPKM计数值相同，SampleA与XCR1(5.5/1,000,000)相关的计数比例也比与sampleB(5.5/1,500,000)相关的计数比例更高。因此，我们不能直接比较sampleA和sampleB中XCR1(或其他基因)的计数，因为样本间归一化计数的总数是不同的。
注:StatQuest的这个视频[1]更详细地展示了为什么应该使用TPM来代替RPKM/FPKM，如果需要对测序深度和基因长度进行标准化。
DESeq2-归一化计数：比值中位数法
由于差异表达分析的工具是比较同一基因的样本组之间的计数，分析工具不需要考虑基因长度。然而，分析确实需要考虑测序深度和RNA组成。
为了对测序深度和RNA组成进行归一化，DESeq2使用了比率中位数法。在用户端只有一个步骤，但在后端有多个步骤，如下所述。
注意：下面的步骤详细描述了当你运行单个函数来获取DE基因时DESeq2执行的一些步骤。基本上，对于典型的RNA-seq分析，你不会单独运行这些步骤。
步骤1：创建伪引用样本（行几何平均值）
对于每个基因，创建一个伪参考样本，它等于所有样本的几何平均值。
https://zhuanlan.zhihu.com/p/23809612
gene | sampleA | sampleB | pseudo-reference sample |
|---|---|---|---|
EF2A | 1489 | 906 | sqrt(1489 * 906) = 1161.5 |
ABCD1 | 22 | 13 | sqrt(22 * 13) = 17.7 |
… | … | … | … |
步骤2：计算每个样本与参考样本的比值
对样本中的每个基因，计算其比值(sample/ref)（如下图所示）。这将对数据集中的每个示例执行。由于大多数基因没有差异表达，所以每个样本中的大多数基因在样本内的比例应该是相似的。
gene | sampleA | sampleB | pseudo-reference sample | ratio of sampleA/ref | ratio of sampleB/ref |
|---|---|---|---|---|---|
EF2A | 1489 | 906 | 1161.5 | 1489/1161.5 = 1.28 | 906/1161.5 = 0.78 |
ABCD1 | 22 | 13 | 16.9 | 22/16.9 = 1.30 | 13/16.9 = 0.77 |
MEFV | 793 | 410 | 570.2 | 793/570.2 = 1.39 | 410/570.2 = 0.72 |
BAG1 | 76 | 42 | 56.5 | 76/56.5 = 1.35 | 42/56.5 = 0.74 |
MOV10 | 521 | 1196 | 883.7 | 521/883.7 = 0.590 | 1196/883.7 = 1.35 |
… | … | … | … |
步骤3：计算每个样本的归一化因子（尺度因子，size factor）
将给定样本中所有比率的中值(上表按列计算)作为该样本的标准化因子（尺度因子），如下所计算。注意差异表达的基因不应影响中位数：
normalization_factor_sampleA <- median(c(1.28, 1.3, 1.39, 1.35, 0.59))
normalization_factor_sampleB <- median(c(0.78, 0.77, 0.72, 0.74, 1.35))
下图说明了单个样本中所有基因比率分布的中值（频率在y轴上）。
img
比率中位数法假设并非所有基因都有差异表达；因此，归一化因子应考虑到样本的测序深度和RNA组成（大的异常值基因不会代表中值比率值）。该方法对上调/下调失衡和大量差异表达基因具有较强的对抗作用。
通常这些大小因子在1左右，如果你看到样本之间有很大的差异，注意这一点很重要，因为这可能表明极端离群值的存在。
步骤4：使用归一化因子计算归一化计数值
这是通过将给定样本中的每个原始计数值除以该样本的标准化因子来生成标准化计数值来实现的。这是对所有计数值(每个样本中的每个基因)执行的。例如，如果SampleA的中值比为1.3,SampleB的中值比为0.77，则可以按如下方式计算归一化计数：
SampleA median ratio = 1.3
SampleB median ratio = 0.77
Raw Counts
gene | sampleA | sampleB |
|---|---|---|
EF2A | 1489 | 906 |
ABCD1 | 22 | 13 |
… | … | … |
Normalized Counts
gene | sampleA | sampleB |
|---|---|---|
EF2A | 1489 / 1.3 = 1145.39 | 906 / 0.77 = 1176.62 |
ABCD1 | 22 / 1.3 = 16.92 | 13 / 0.77 = 16.88 |
… | … | … |
请注意，标准化的计数值不是整数。
使用DESeq2对Mov10数据集进行计数标准化
现在我们已经了解了计数归一化理论，接下来我们将使用DESeq2归一化Mov10数据集的计数。这需要几个步骤：
- 确保出现元数据数据框有行名，并且与计数数据框的列名顺序相同。
- 创建一个
DESeqDataSet对象。 - 生成标准化计数
1. 匹配元数据和计数数据
我们应该始终确保示例名称在两个文件之间匹配，并且示例的顺序正确。如果不是这样，DESeq2将输出一个错误。
### Check that sample names match in both files
all(colnames(txi$counts) %in% rownames(meta))
all(colnames(txi$counts) == rownames(meta))
如果数据不匹配，可以使用match()函数重新排列它们以匹配。
2. 创建DESeq2对象
Bioconductor软件包通常在R中定义和使用一个自定义类来存储数据(输入数据、中间数据和结果)。这些自定义数据结构与列表相似，因为它们可以包含多种不同的数据类型/结构。但是，与列表不同的是，它们有预先指定的数据槽，用于存放特定类型/类的数据。存储在这些预先指定槽位中的数据可以通过使用特定的包定义函数来访问。
让我们从创建DESeqDataSet对象开始，然后我们可以更多地讨论它里面存储的内容。为了创建对象，我们将需要计数矩阵和元数据表作为输入。我们还需要指定一个设计公式。设计公式指定元数据表中的列，以及在分析中应该如何使用这些列。对于我们的数据集，我们只对一个列感兴趣，即~sampletype。这一栏有三个因子水平，它告诉DESeq2，对于每个基因，我们想要评估这些不同水平的基因表达变化。
我们的计数矩阵输入存储在txi列表对象中，因此我们使用DESeqDataSetFromTximport()函数传递它，该函数将提取计数部分并将值四舍五入到最接近的整数中。
## Create DESeq2Dataset object
dds <- DESeqDataSetFromTximport(
txi,
colData = meta,
design = ~ sampletype)
注意：由于我们在上一章节中创建了一个包含计数的数据变量，我们也可以使用它作为输入。但是，在这种情况下，我们希望使用DESeqDataSetFromMatrix()函数。
img
如果愿意，可以使用特定于DESeq的函数来访问不同的数据槽和检索信息。例如，假设我们想要原始的计数矩阵，我们将使用counts()(注意:我们将其嵌套在View()函数中，这样我们就可以在脚本编辑器中看到它，而不是在控制台中打印)：
View(counts(dds))
在我们完成工作流的过程中，我们将使用相关函数来检查对象中存储了哪些信息。
3. 生成Mov10标准化计数
下一步是标准化计数数据，以便能够在样本之间进行公平的基因比较。
img
为了执行归一化的比率中值方法，DESeq2有一个estimateSizeFactors()函数，它将为我们生成大小因子。我们将在下面的例子中使用这个函数，但是在一个典型的RNA-seq分析中，这个步骤是由DESeq()函数自动执行的，我们将在后面看到。
dds <- estimateSizeFactors(dds)
通过将结果赋值给dds对象，我们将用适当的信息填充DESeqDataSet对象的槽。我们可以使用以下方法查看每个样本的归一化因子：
sizeFactors(dds)
#Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 Mov10_oe_2 Mov10_oe_3
# 1.1150371 0.9606366 0.7493552 1.5634128 0.9359082 1.2257749 1.1406863 0.6541689
现在，要从dds检索标准化计数矩阵，我们使用counts()函数并添加参数normalized=TRUE。
normalized_counts <- counts(dds, normalized=TRUE)
我们可以将这个标准化的数据矩阵保存到文件中以备以后使用：
write.table(normalized_counts, file="data/normalized_counts.txt",
sep="\t", quote=F, col.names=NA)
注意：DESeq2实际上并不使用标准化计数，而是使用原始计数，并在广义线性模型(GLM)中对标准化进行建模。这些标准化的计数对于结果的下游可视化是有用的，但是不能作为DESeq2或任何其他使用负二项模型进行差异表达分析的工具的输入。
参考资料
[1]
StatQuest的这个视频: http://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/