R语言层次聚类、多维缩放MDS分类RNA测序（RNA-seq）乳腺发育基因数据可视化|附数据代码

全文链接：https://tecdat.cn/?p=35691

在乳腺发育过程中，促生存基因Mcl1被认为是关键的调控因子之一。为了更全面地揭示Mcl1在乳腺发育中的调控作用，本研究帮助客户采用层次聚类和多维缩放（MDS）等方法对RNA-seq数据进行深入分析。层次聚类可以帮助我们识别不同样本之间的相似性和差异性，从而揭示乳腺发育过程中不同阶段的基因表达模式。而MDS则可以将高维的基因表达数据转化为低维空间中的点，便于我们进行可视化和模式识别。

本文数据来自一项关于乳腺发育的研究，表明促生存基因Mcl1是乳腺发育的主要调控因子。这些数据代表了来自6个组（处女、怀孕和哺乳小鼠的基底细胞和腔细胞）的重复RNA-seq测量。

概述

在这里，我们要达到的最终产物是一个层次聚类。

试剂是一个包含每个基因每个重复样本的原始RNAseq计数的文件。我们还有一个包含有关每个重复样本的额外信息的文件……即元数据。我们需要将第二个文件的元数据合并到第一个文件中，即基因表达数据。至于计数，我们需要去除缺失或低值的基因。最终，计数将进行对数转换，从而变成高斯分布的数据。这对于聚类（或其他统计模型）是必要的，因为这些模型假设输入数据是高斯分布的。

导入原始计数数据

到这个阶段，原始的序列数据已经经过了大量的生物信息学处理。序列已经被匹配到它们各自的基因和彼此之间。这些方法超出了我们在这里想要涵盖的范围。

RNAseq核心设施很可能已经为客户完成了生物信息学分析。通常，他们会为研究人员提供一个计数文件，如这个，作为最低限度的交付成果。

该表格中的值是每个基因在每个12个样本中的转录读取数。它们采用定界文件格式，因此我们使用read.delim函数。

将这些数据读入一个对象中。

首先，让我们花些时间检查数据。

查看数据集的维度。

有14列和超过27000行。

head(seqdata)

我们可以看到其中有两列与其他列不同。我们将处理这些列，以创建一个具有行名和列名，并且所有列中的数值都是序列计数的数据对象。

另一个分隔的文本文件包含有关每个样本重复的一些额外的因子信息，包括一些重要的独立变量。这些信息将在后面的聚类分析解释中作为一个重要的方面。但现在我们先将其读入环境并查看一下。

接下来的步骤将处理这些数据特征。

过滤

在继续分析之前，我们通常需要过滤掉一些基因。这可能是出于多种原因，例如：

低表达基因：这些基因在所有样本中的表达量都很低，可能是噪声或无关紧要的。

没有变异的基因：如果某个基因在所有样本中的表达量都是相同的，那么它对于区分样本类型或条件可能没有帮助。

过滤步骤可以帮助我们集中精力在更有意义的基因上，从而提高后续分析的准确性和效率。通常，我们会根据基因的表达水平、变异程度或其他相关指标来设定过滤条件。

接下来，我们需要过滤掉那些没有读数、重复样本间读数不一致或读数较低的基因。

这是一个多步骤的过程。

第一步是选择归一化技术。RPKM（每千碱基每百万读数）和CPM（每百万计数）是常见的选项。我们将使用后者。

我们将使用edgeR来完成这个步骤，这是我们第一个使用的Bioconductor函数。

我们的过滤规则是保留至少在两个样本中CPM > 0.5的转录本。在这个大小的文库中，CPM为0.5大约对应于每个基因10-15个计数。这个阈值的决定是一个判断问题。低计数往往不可靠。将阈值设定为1或2 CPM也是常见的。研究者在选择CPM与RPKM以及阈值水平选项时有一定的自由度。这两方面都没有普遍的共识。

Qing Li

拓端分析师

首先，将原始计数转换为CPM并查看。

myCPM <- edgeR::cpm(countdata)

cpm函数不仅进行了归一化，还将countdata对象从数据框（data frame）转换为了矩阵（matrix）。这样做是因为在后续的基因表达分析中，特别是使用edgeR这样的包时，通常需要矩阵格式的数据。通过转换，我们可以更高效地处理和分析大量的基因表达数据。

接下来，我们要应用阈值。回想一下我们的过滤规则：保留至少在两个样本中CPM > 0.5的转录本。

这个规则包含两个部分。

下面的脚本是一个简单的逻辑判断，用于识别满足过滤规则第一部分的基因和分组。

结果描述

table(rowSums(thrh))

这段输出可以解释如下：

在所有12个样本中，有10857个基因的CPM值≤0.5。有518个基因仅在1个样本中的CPM值大于0.5。544个基因仅在2个样本中的CPM值大于0.5，307个基因在3个样本中的CPM值大于0.5，以此类推。有11433个基因在所有12个样本中的CPM值都大于0.5。

接下来，我们要识别那些在至少两个12个重复样本中满足阈值的基因。这是另一个逻辑判断。它会产生一个长的逻辑向量，其中每个行名对应一个True（真）或False（假）。

以下是仅包含这些经过过滤的基因的更新后的计数数据集。这是将用于统计分析的最终过滤数据集。

请注意，我们是如何使用向量keep作为行索引的。用通俗的话说，countdata[keep,]给我们提供了keep中基因ID值为TRUE的每一行。

同时请注意，counts.keep是一个数据框。

在继续之前，请确保你已经保存了过滤后的数据集，因为这将是你接下来进行差异表达分析或其他统计测试的基础。你可以使用write.csv或类似的函数将数据框保存为CSV文件，以便后续使用。

此外，如果后续步骤需要使用归一化后的数据（如CPM值），请确保也保存了经过过滤的CPM数据集。这将有助于你在分析过程中保持数据的一致性，并避免在多个步骤之间重复进行相同的过滤和归一化操作。

counts.keep数据框下面被转换为名为y的对象，使用了DGEList函数。

DGEList列表对象是我们在课程中之前未见过的一个R类对象。这些对象在一个列表项中携带计数数据，同时在其他列表项中携带其他“元数据”信息。

例如，在y中，列表项y$counts是一个包含计数数据的矩阵。

使用DGEList函数创建对象y是为了适应edgeR包的分析流程。edgeR是一个专门用于差异表达分析的R包，它要求数据以DGEList对象的形式输入。这个对象结构不仅包含了原始的计数数据，还可以包含实验设计信息、样本分组等元数据，这些对于后续的统计分析和模型拟合至关重要。

在创建y对象后，你通常会继续进行诸如TMM（Trimmed Mean of M-values）归一化、差异表达分析、模型拟合等步骤。这些步骤都是edgeR包的核心功能，用于识别不同条件下基因表达水平的变化。

简而言之，将counts.keep转换为DGEList对象是为了准备数据以便进行后续的差异表达分析，这是RNA-seq数据分析中的一个关键步骤。

y$sample这个列表项是一个包含样本元数据的数据框。

在edgeR的DGEList对象中，y$sample用于存储样本相关的信息，如样本名称、分组、条件等。这些信息对于后续的统计分析至关重要，因为它们被用来构建模型，确定哪些样本属于哪个比较组，以及哪些样本应该用于计算差异表达。

DGEList对象被传递给各种Bioconductor函数以完成各种任务。

在某种程度上，就像数据框对象是tidyverse的核心一样，DGEList对象是用于RNA-seq分析函数的核心。

可视化

许多组学可视化使用R基础绘图函数。我们大部分时间都避免使用基础绘图。但是，利用手头的数据，使用这些绘图函数得到一些结果是相当简单的。

现在，让我们制作视图以快速检查RNA-seq数据。

比较文库大小可以检查是否存在异常。例如，我们不希望某个库的读数显著少于其他库。这里没有问题。每个样本的库大小都差不多。

barplot(y$saolnames(y),las=2)

下面是每个基因的原始计数图，以说明它们不是正态分布，这是过度分散的离散计数数据的典型特征。

RNA-seq数据通常具有偏态分布，即计数数据在许多基因中可能很低（接近于零），而在少数基因中可能很高。这种分布模式不适合正态分布假设的许多统计方法。因此，在分析RNA-seq数据时，我们通常会使用专门为计数数据设计的统计模型，如负二项分布模型，这些模型能够处理这种过度分散的特性。

Qing Li

拓端分析师

以下是数据作为CPM（每百万计数）的即时归一化图，它也不是正态分布。这说明了CPM只是计数数据的一个简单线性变换。