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

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

原文出处：拓端数据部落公众号

分析师：Qing Li

在生物学和医学研究中，乳腺发育是一个复杂而精细的过程，涉及众多基因的表达调控。近年来，随着高通量测序技术的发展，RNA测序（RNA-seq）技术已经成为研究基因表达模式的有力工具。通过RNA-seq技术，我们可以获得大量关于基因表达的定量信息，从而更深入地理解乳腺发育的分子机制。

在乳腺发育过程中，促生存基因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以及阈值水平选项时有一定的自由度。这两方面都没有普遍的共识。

首先，将原始计数转换为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是一个数据框。

另外，请注意我们又回到了计数数据！我们使用CPM归一化作为过滤低表达基因的手段。接下来，我们将开始分析的下一步骤。

在继续之前，请确保你已经保存了过滤后的数据集，因为这将是你接下来进行差异表达分析或其他统计测试的基础。你可以使用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数据时，我们通常会使用专门为计数数据设计的统计模型，如负二项分布模型，这些模型能够处理这种过度分散的特性。

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

在RNA-seq分析中，对原始计数数据进行归一化是一个重要的步骤，因为它可以帮助消除由于测序深度、文库大小或批次效应等因素导致的差异。CPM（每百万计数）是一种简单的归一化方法，它将每个样本的原始计数除以该样本中所有基因计数的总和，并乘以一百万，以得到每个基因在每个样本中的相对表达量。

现在让我们绘制CPM数据的对数变换图。这种变换使得每个样本内的计数值大致呈现高斯分布，尽管显然存在离群值。

由于这种高斯特性，对数变换后的CPM数据将在下面的许多统计建模中使用。因此，我们将继续并为其创建一个对象。

首先，我们将创建一个对数变换后的对象。

现在，让我们查看logcount数据，通过绘制箱线图来进行。

    # 使用箱线图检查样本的分布  

    boxplot(logcounts, xlab="", 

多维缩放分类

现在，我们可以进行一些统计建模了。

多维缩放是主成分分析的一个类似物。它提供了一种简单的方法，可以查看样本组是否沿着其第一和第二维度分开。这是基于主要差异倍数度量的。

这检查了在样本之间表现出最大倍数差异的基因子集。

下面的图显示了12个样本的分离和聚类，但除非你还记得每个样本代表什么，否则很难看出这向我们展示了什么。

下面我们将使用颜色来制作一些更好的图，重点关注预测变量的分类。

回想一下，有12个样本，每个样本都是基细胞和腔细胞中的基因表达重复，分别对应于以下三种条件：细胞来源于处女、怀孕或哺乳的小鼠。

head(sampleinfo, n=8)

首先，我们将根据涉及基细胞和腔细胞的特征进行着色。

复制代码
    # 为 CellType 设置颜色方案  

    # 有多少种细胞类型，以及它们是以什么顺序存储的？  

    levels(as.fafo$CellType))

下面的脚本中的第二行只是检查第一行是否正确创建了col.cell向量。

多维缩放检测最大的变异来源。在这种情况下，它检测与在细胞类型之间表现出最大倍数差异的基因相关的变异。

因此，下面的图清晰地表明，负责基细胞/腔细胞特征的基因代表了样本之间分离的第一维度。

换句话说，我们可以说基细胞和腔细胞基因表达之间的差异占据了数据集中大部分变异。RNA-seq数据中有信号可以将这两种细胞类型区分开。

MDS图显示了由细胞类型导致的变异解释了数据的第一个维度。

首先，Status是什么？

当然，从这三种内分泌状态中的每一种，研究人员都分别收集了基细胞和腔细胞的两个样本。

下面的图现在根据样本是来自处女、怀孕还是哺乳的小鼠进行着色。

它们在第二维度上有所分离。请注意，同一条件的重复样本非常相似（MCL1.DL和MCL1.DK几乎重叠在一起）。

因此，基因表达的最大变异维度是基于细胞类型进行分离的，而第二大的变异维度则是基于内分泌状态进行分离的。

尽管这些维度是潜在的，但我们可以解释哪些变量最负责观察到的变异，从而解释前两个主成分。这里的变量选择得很好：细胞类型和内分泌状态。这确实是一个很好的结果。

层次聚类

层次聚类是一种在聚类时可视化特定基因与其他聚类之间关系的方法。

在说明这种技术之前，需要进行一些处理。

我们感兴趣的是查看驱动上述维度模式的基因。哪些基因在两种细胞类型之间存在差异？哪些基因在三种状态条件之间存在差异？细胞类型和状态之间是否存在交互作用？

回答这类问题的最有用方法是将焦点放在表达方差最高的基因上。这些基因在整个12个样本中的表达变化最大。

与其对所有15000+个基因进行聚类（这只会给分析增加噪音），我们不如对变化最大的500个基因进行聚类。这是一个任意的截止值，可以根据需要进行调整。

下面的脚本创建了一个与每个GeneID相关联的方差向量。就像ANOVA一样，当行方差较高时，我们有迹象表明分组因子之间的差异将最大。

这个函数只是简单地计算logcounts数据框中每一行的方差。这只是一个数值向量。

var_gnes <- apply(logouts, 1, var)
head(var_gees)

现在，我们创建一个字符向量，其中包含方差从最大到最小的那些GeneID的名称。在这个例子中，我们将选择前50个。

在查看热图后，改变这个数字以选择更高或更低的值可能是有指导意义的。

与其他的相比，处女基细胞中的Wap表达量要低得多。只要记住这些是log counts（对数计数），这个结果就通过了显而易见的检验。

要从logcounts中选择这500个变化最大的基因对应的行，你可以使用前面创建的select_var字符向量来索引数据框。这里我们实际上只需要前50个，但如果你想选择前500个，只需将select_var的长度调整为500即可。

现在，创建层次聚类。

接下来，我们只需将这组500个变化最大的基因选入heatmap.2函数中。

这里表示的值是logCPM值。

 

表达量从低（深色）到高（浅色）不等。

观察水平聚类可以很好地揭示实验设计。存在两个主要组（对应于腔细胞和基底细胞类型）。在每个主要组内，还分别存在3个小组，对应于内分泌状态。

内分泌状态的重复行排列得非常好。这说明数据非常紧密。

垂直聚类很有意思。大约90%的高度可变的基因决定了细胞类型的分化，而剩下的基因则区分了内分泌状态（处女、怀孕、哺乳）。此外，细胞类型和状态之间也存在明显的相互作用。

---

最受欢迎的见解

1.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e7%bb%98%e5%88%b6%e7%94...)R语言分布式滞后非线性模型（DLNM）空气污染研究温度对死亡率影响建模

2.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e7%94%9f%e5%ad%98%e5%88...)R语言分布滞后线性和非线性模型（DLNM）分析空气污染（臭氧）、温度对死亡率时间序列数据

3.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e5%a6%82%e4%bd%95%e5%9c...)R语言群组变量选择、组惩罚GROUP LASSO套索模型预测分析新生儿出生体重风险因素数据和交叉验证、可视化

4.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e4%b8%ad%e4%bd%bf%e7%94...)R语言逻辑回归、随机森林、SVM支持向量机预测FRAMINGHAM心脏病风险和模型诊断可视化

5.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e7%94%9f%e5%ad%98%e5%88...)R语言非线性混合效应 NLME模型(固定效应&随机效应)对抗哮喘药物茶碱动力学研究

6.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80ggplot2%e8%af%af%e5%b7%a...)R语言使用限制平均生存时间RMST比较两条生存曲线分析肝硬化患者

7.[](https://tecdat.cn/r-%e8%af%ad%e8%a8%80%e7%bb%98%e5%88%b6%e5%8...)分类回归决策树交互式修剪和更美观地可视化分析细胞图像分割数据集

8.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e5%a6%82%e4%bd%95%e6%89...)PYTHON深度学习实现自编码器AUTOENCODER神经网络异常检测心电图ECG时间序列

9.[](https://tecdat.cn/r%e8%af%ad%e8%a8%80%e4%b8%ad%e7%9a%84%e7%94...)R语言如何在生存分析与Cox回归中计算IDI，NRI指标