跟着存档教程动手学RNAseq分析（四）：使用DESeq2进行DE分析的QC方法

跟着存档教程动手学RNAseq分析（一）

跟着存档教程动手学RNAseq分析（二）

跟着存档教程动手学RNAseq分析（三）：使用DESeq2进行计数标准化

质量控制QC

DESeq2工作流程中的下一个步骤是QC，它包括对计数数据执行样本级和基因级QC检查的步骤，以帮助我们确保样本/重复看起来良好。

img

样本水平QC

RNA-seq分析的一个有用的初始步骤通常是评估样本之间的整体相似性：

• 哪些样本相似，哪些不同？
• 这符合实验设计的期望吗？
• 数据集中变异的主要来源是什么？

为了探索我们的样本的相似性，我们将使用主成分分析（PCA）和层次聚类方法来执行样本级QC。我们的样本水平QC让我们可以看到我们的重复聚在一起的情况，以及观察我们的实验条件是否代表数据中变化的主要来源。执行样本级QC还可以识别任何样本离群值，这可能需要进一步研究，以确定它们是否需要在DE分析之前删除。

img

当使用这些非监督聚类方法时，标准化计数的log2转换可以提高可视化的距离/聚类。DESeq2对样本级QC使用标准化计数的正则化对数变换(rlog)，因为它调节了均值间的方差，从而改进了聚类。

img

注意:DESeq2文档建议大数据集（100个样本）使用方差稳定转换(vst)而不是rlog来进行计数转换，因为rlog函数可能运行太长时间，而vst()函数具有与rlog相似的属性，速度更快。

主成分分析PCA[1]

主成分分析（PCA）是一种技术，用于强调变化，并提出数据集中强大的模式（降维）。关于PCA的细节如下所示（基于来自StatQuest的材料[2]，如果你想要更全面的描述，我们鼓励你探索StatQuest的视频[3]和我们更长的课程[4]）。

假设我们有一个包含两个样本和四个基因的数据集。基于这个表达式数据，我们想要评估这些样本之间的关系。我们可以绘制一个样本与另一个样本的计数关系，样本1在x轴上，样本2在y轴上，如下所示：

img

对于PCA分析，第一步是绘制这个图，并在代表变化最多的方向上通过数据画一条线。在本例中，沿对角线变化最多。也就是说，数据中最大的分布在这条线的两个端点之间。这被称为第一个主成分，或PC1。这条线两端的基因(基因B和基因C)对这条线的方向影响最大。

img

在绘制这条线并确定每个基因的影响量之后，PCA将计算每个样本的得分。每个样本的PC1评分是通过将影响和标准化计数的乘积以及所有基因的总和来计算的。我们可以通过表示数据(PC2)中第二大变化量的数据绘制另一条线，然后计算分数，然后是第三条线，以此类推，直到数据集中的样本总数。

Sample1 PC1 score = (read count Gene A * influence Gene A) + 
                    (read count Gene B * influence Gene B) + 
                    .. for all genes

实际上，你的实际数据集维度更大（更多的样本和更多的基因）。初始的样本-样本图，将在n维空间中n个轴代表样本的总数。最终结果是一个二维矩阵，其中行表示样本，列反映每个主成分的分数。为了评估主成分分析的结果，我们通常将主成分相互比拼，从解释数据中最大量变化的pc开始。

img

如果两个样本中对PC1所代表的变异有显著贡献的基因表达水平相似，那么它们将在PC1轴上紧密地绘制在一起。因此，我们预计生物学重复具有相似的得分(因为相同的基因发生改变)，并聚集在PC1和/或PC2上，来自不同处理组的样本具有不同的得分。这是最容易理解的可视化示例PCA图。

解释PCA图

下面我们有一个示例数据集和一些相关的PCA图，以了解如何解释它们。实验的元数据显示在下面。主要感兴趣的条件是treatment。

img

当在PC1和PC2上可视化时，我们没有看到通过处理分离的样本，所以我们决定探索数据中存在的其他变异来源。我们希望我们已经在元数据表中包含了所有可能的已知变异源，并且我们可以使用这些因素来为PCA图着色。

img

我们从因子cage开始，但cage因子似乎不能解释PC1或PC2上的变化。

img

然后，我们根据性别因素着色，这似乎是在PC2上分离样品。这是值得注意的良好信息，因为我们可以使用它来解释模型中由于性别而产生的变化，并将其回归建模出来。

img

接下来我们探讨strain因子，发现它可以解释PC1上的变化。

img

我们已经能够确定PC1和PC2的变异来源，这是非常棒的。通过在我们的模型中考虑到这一点，我们应该能够检测到更多因治疗而有差异表达的基因。

令人担忧的是，我们看到两个样本没有与正确的菌株聚集在一起。这将表明可能进行样品交换，并应进行调查，以确定这些样品是否确实是标记的菌株。如果我们发现存在（错误的）交换，我们可以交换元数据中的样本。然而，如果我们认为它们被正确标记或不确定，我们可以从数据集中删除样本。

但我们仍然没有发现，治疗是否是strain和性别后变异的主要来源。所以，我们探索PC3和PC4来看看治疗是否驱动了这些PC所代表的变化。

img

我们发现在PC3上处理分离的样品，我们对我们的DE分析很乐观，因为我们感兴趣的条件，处理，是在PC3上分离的，我们可以回归出驱动PC1和PC2的变异。

根据前几个主成分解释了多少变化，你可能想要探索更多(即考虑更多成分并绘制成对组合)。即使你的样本不能被实验变量清楚地分开，你仍然可以从DE分析中得到生物学上相关的结果。如果你期望的效应量非常小，那么信号可能会被外来的变化源所淹没。在你可以确定这些来源的情况下，在你的模型中考虑这些来源是很重要的，因为它为检测DE基因的工具提供了更强大的功能。

层次聚类的热图

与主成分分析相似，层次聚类是另一种用于识别数据集中的强模式和潜在异常值的补充方法。热图显示了数据集中所有成对组合的样本的基因表达的相关性。由于大多数基因没有差异表达，所以样本之间的相关性一般较高(值大于0.80)。低于0.80的样品可能表明你们的数据中存在异常值和/或样品污染。

层次树可以根据归一化的基因表达值指出哪些样本彼此更相似。颜色块表示数据中的子结构，您可能会看到每个示例组的复制聚在一起作为一个块。此外，我们希望看到聚集的样本类似于在PCA图中观察到的分组。在下面的图中，我们将非常关注‘Wt_3’和‘KO_3’的样本与其他重复的样本没有聚类。我们想要探索主成分分析，看看我们是否看到了相同的样本聚类。

img

基因水平QC

除了检查样品/重复聚类是否良好，还有更多的QC步骤。在进行差异表达分析之前，省略很少或没有机会被检测出差异表达的基因是有益的。这将提高检测差异表达基因的能力。被遗漏的基因分为三类：

• 在所有样本中计数为0的基因
• 有极端异常值的基因
• 低均值标准化计数的基因

img

默认情况下，DESeq2将执行此过滤：而其他DE工具，如edgeR则不会。即使使用limma-voom和/或edgeR的准似然方法，过滤也是必要的一步。在使用其他工具时，请确保遵循预过滤步骤，正如Bioconductor上的用户指南中概述的那样，因为它们通常性能更好。

使用DESeq2进行Mov10质量评估和探索性分析

现在我们已经很好地理解了通常用于RNA-seq的QC步骤，让我们为将要使用的Mov10数据集实现它们。

使用rlog转换标准化计数

为了改进PCA和分层聚类可视化方法的距离/聚类，我们需要通过对标准化计数应用rlog变换来调节均值方差。

在质量评估期间，标准化计数的rlog转换仅对这些可视化方法是必要的。我们不会在下游使用这些标准计数。

### Transform counts for data visualization
rld <- rlog(dds, blind=TRUE)

blind=TRUE参数将导致对样本条件信息无偏倚的转换。在执行质量评估时，包含此选项是很重要的。DESeq2[5]文档有更多的细节。

rlog函数返回一个DESeqTransform对象，另一种特定的DESeq的对象类型。你不只是得到一个转换后的值的矩阵的原因是，计算rlog转换的所有参数(即大小因子)都存储在该对象中。我们使用这个对象来绘制质量评估的主成分分析和层次聚类图。

注意:当你有例如> 20个样本时，rlog()函数可能会有点慢。在这些情况下，vst()函数要快得多，并执行类似的转换，适合与plotPCA()一起使用。由于优化和转换的性质，使用vst()通常只需要几秒钟。

主成分分析(PCA)

DESeq2有一个用于绘制PCA图的内置函数，它在底层使用ggplot2。这是非常棒的，因为它节省了我们输入代码行和摆弄不同ggplot2层的时间。此外，它直接将rlog对象作为输入，从而省去了从其中提取相关信息的麻烦。

### Plot PCA 
plotPCA(rld, intgroup="sampletype")

plotPCA()函数需要两个参数作为输入rlog对象和intgroup(我们感兴趣的元数据中的列)。

img

关于样本的相似性，这张图告诉了你什么?它符合实验设计的期望吗?默认情况下，该函数使用前500个最可变的基因。您可以通过添加ntop参数并指定要使用多少个基因来绘制图表来改变这一点。

注意：plotPCA()函数将只返回PC1和PC2的值。如果你想在数据中探索其他的pc，或者如果你想确定对这些pc起主要作用的基因，你可以使用prcomp()函数。例如，要绘制任意一个pc，我们可以运行以下代码：

# Input is a matrix of log transformed values
rld <- rlog(dds, blind=T)
rld_mat <- assay(rld)
pca <- prcomp(t(rld_mat))

# Create data frame with metadata and PC3 and PC4 values for input to ggplot
df <- cbind(meta, pca$x)
ggplot(df) + geom_point(aes(x=PC3, y=PC4, color = sampletype))

可以使用资源[6]学习如何使用PC进行更复杂的查询。

分层聚类

由于在DESeq2中没有针对热图的内置函数，我们将使用pheatmap包中的pheatmap()函数。这个函数需要一个由数值组成的矩阵/数据帧作为输入，因此我们需要做的第一件事就是从rld对象中获取信息：

### Extract the rlog matrix from the object
rld_mat <- assay(rld)    
## assay() is function from the "SummarizedExperiment" package that was loaded when you loaded DESeq2

然后我们需要计算样本的两两相关值。我们可以使用cor()函数来做到这一点：

### Compute pairwise correlation values
rld_cor <- cor(rld_mat)    ## cor() is a base R function

head(rld_cor)   ## check the output of cor(), make note of the rownames and colnames

现在将相关值绘制成热图：

### Load pheatmap package
library(pheatmap)

### Plot heatmap
pheatmap(rld_cor, annotation = meta)

img

总体而言，我们观察到整体相关性相当高(> 0.999)，表明没有离群样本。此外，与PCA图类似，你可以看到样本按样本组聚类在一起。总之，这些图向我们表明数据质量良好，我们可以进行差异表达分析。

注意：pheatmap函数有许多不同的参数，我们可以通过改变默认值来增强图形的美观性。如果你感到好奇并想了解更多，请尝试运行下面的代码。你的图形是如何变化的？查看帮助页面(?pheatmap)并确定每个添加的参数对图形的贡献。

heat.colors <- brewer.pal(6, "Blues")
pheatmap(rld_cor, annotation = meta, color = heat.colors, border_color=NA, fontsize = 10, 
   fontsize_row = 10, height=20)

想知道RColorBrewer包提供的所有可用的调色板吗？试着在控制台中输入display.brewer.all()，看看会发生什么！

参考资料

[1]

主成分分析PCA: https://hbctraining.github.io/DGE_workshop/lessons/principal_component_analysis.html

[2]

StatQuest的材料: https://www.youtube.com/watch?v=_UVHneBUBW0

[3]

StatQuest的视频: https://www.youtube.com/watch?v=_UVHneBUBW0

[4]

更长的课程: https://hbctraining.github.io/scRNA-seq/

[5]

DESeq2: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#blind-dispersion-estimation

[6]

资源: http://www.sthda.com/english/wiki/principal-component-analysis-in-r-prcomp-vs-princomp-r-software-and-data-mining