R语言层次聚类、多维缩放MDS分类RNA测序（RNA-seq）乳腺发育基因数据可视化|附数据代码1:https://developer.aliyun.com/article/1501241
可视化
许多组学可视化使用R基础绘图函数。我们大部分时间都避免使用基础绘图。但是，利用手头的数据，使用这些绘图函数得到一些结果是相当简单的。
现在，让我们制作视图以快速检查RNA-seq数据。
比较文库大小可以检查是否存在异常。例如，我们不希望某个库的读数显著少于其他库。这里没有问题。每个样本的库大小都差不多。
barplot(y$saolnames(y),las=2)
下面是每个基因的原始计数图，以说明它们不是正态分布，这是过度分散的离散计数数据的典型特征。
RNA-seq数据通常具有偏态分布，即计数数据在许多基因中可能很低（接近于零），而在少数基因中可能很高。这种分布模式不适合正态分布假设的许多统计方法。因此，在分析RNA-seq数据时，我们通常会使用专门为计数数据设计的统计模型，如负二项分布模型，这些模型能够处理这种过度分散的特性。
Qing Li
拓端分析师
以下是数据作为CPM（每百万计数）的即时归一化图，它也不是正态分布。这说明了CPM只是计数数据的一个简单线性变换。
在RNA-seq分析中，对原始计数数据进行归一化是一个重要的步骤，因为它可以帮助消除由于测序深度、文库大小或批次效应等因素导致的差异。CPM（每百万计数）是一种简单的归一化方法，它将每个样本的原始计数除以该样本中所有基因计数的总和，并乘以一百万，以得到每个基因在每个样本中的相对表达量。
现在让我们绘制CPM数据的对数变换图。这种变换使得每个样本内的计数值大致呈现高斯分布，尽管显然存在离群值。
由于这种高斯特性，对数变换后的CPM数据将在下面的许多统计建模中使用。因此，我们将继续并为其创建一个对象。
首先，我们将创建一个对数变换后的对象。
现在，让我们查看logcount数据，通过绘制箱线图来进行。
# 使用箱线图检查样本的分布 boxplot(logcounts, xlab="",
多维缩放分类
现在，我们可以进行一些统计建模了。
多维缩放是主成分分析的一个类似物。它提供了一种简单的方法，可以查看样本组是否沿着其第一和第二维度分开。这是基于主要差异倍数度量的。
这检查了在样本之间表现出最大倍数差异的基因子集。
下面的图显示了12个样本的分离和聚类，但除非你还记得每个样本代表什么，否则很难看出这向我们展示了什么。
下面我们将使用颜色来制作一些更好的图，重点关注预测变量的分类。
回想一下，有12个样本，每个样本都是基细胞和腔细胞中的基因表达重复，分别对应于以下三种条件：细胞来源于处女、怀孕或哺乳的小鼠。
head(sampleinfo, n=8)
首先，我们将根据涉及基细胞和腔细胞的特征进行着色。
复制代码 # 为 CellType 设置颜色方案 # 有多少种细胞类型，以及它们是以什么顺序存储的？ levels(as.fafo$CellType))
下面的脚本中的第二行只是检查第一行是否正确创建了col.cell向量。
现在我们将为细胞类型在MDS图上添加颜色方案。
多维缩放检测最大的变异来源。在这种情况下，它检测与在细胞类型之间表现出最大倍数差异的基因相关的变异。
因此，下面的图清晰地表明，负责基细胞/腔细胞特征的基因代表了样本之间分离的第一维度。
换句话说，我们可以说基细胞和腔细胞基因表达之间的差异占据了数据集中大部分变异。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%的高度可变的基因决定了细胞类型的分化，而剩下的基因则区分了内分泌状态（处女、怀孕、哺乳）。此外，细胞类型和状态之间也存在明显的相互作用。
