重生之我在剑桥大学学习单细胞RNA-seq分析——3. 单细胞分析中的R/Bioconductor简介（4）

3.13 ggplot2介绍
3.13.1 什么是ggplot2
ggplot2是一个由Hadley Wickham设计的R包，用于辅助数据绘图。我们将简单介绍该包的一些功能。如果您想了解有关如何使用ggplot2的更多信息，我们建议您阅读Hadley Wickham撰写的“ggplot2 用于数据分析的优雅图形”。
3.13.2 ggplot2的原理

• 如果您想使用ggplot2绘制数据，则数据必须是数据框。
• 使用aes函数指定数据框中的变量如何映射到绘图上的特征。
• 使用geoms指定如何在图表上表示数据，例如散点图、条形图、箱线图等。

3.13.3 使用aes映射函数
aes函数指定数据框中的变量如何映射到绘图上的特征。为了理解它的工作原理，让我们看一个例子：

> library(ggplot2)
> library(tidyverse)
> set.seed(1)
> counts <- as.data.frame(matrix(rpois(100, lambda = 10), ncol=10, nrow=10))
> Gene_ids <- paste("gene", 1:10, sep = "")
> colnames(counts) <- paste("cell", 1:10, sep = "")
> counts<-data.frame(Gene_ids, counts)
> counts
   Gene_ids cell1 cell2 cell3 cell4 cell5 cell6 cell7 cell8 cell9 cell10
1     gene1     8     8     3     5     5     9    11     9    13      6
2     gene2    10     2    11    13    12    12     7    13    12     15
3     gene3     7     8    13     8     9     9     9     5    15     12
4     gene4    11    10     7    13    12    12    12     8    11     12
5     gene5    14     7     8     9    11    10    13    13     5     11
6     gene6    12    12    11    15     8     7    10     9    10     15
7     gene7    11    11    14    11    11     5     9    13    13      7
8     gene8     9    12     9     8     6    14     7    12    12     10
9     gene9    14    12    11     7    10    10     8    14     7     10
10   gene10    11    10     9     7    11    16     8     7     7      4
> ggplot(data = counts, mapping = aes(x = cell1, y = cell2))

让我们仔细看看最后一个命令，ggplot(data = counts, mining = aes(x = cell1, y = cell2))。ggplot()初始化一个ggplot对象并接受参数data和mining。我们将counts数据框传递给data，并使用aes()函数指定我们希望使用变量cell1作为x变量，使用变量cell2作为y变量。
任务1：修改上述命令以初始化ggplot对象，其中cell10是x变量，cell8是y变量。
显然，我们刚刚创建的图表信息量不大，因为其中没有显示任何数据。要显示数据，我们需要使用geoms。
3.13.4 Geoms
我们可以使用geoms来指定如何在图表上显示数据。例如，我们对geom的选择可以指定我们希望将数据显示为散点图、条形图或箱线图。
让我们看看我们的图表作为散点图看起来是什么样子。

> ggplot(data = counts, mapping = aes(x = cell1, y = cell2)) + geom_point()

现在我们可以看到，cell1和cell2中的基因表达似乎没有任何相关性。考虑到我们是随机生成计数的，这并不太令人惊讶。
任务2：修改上述命令以创建线图。提示：执行?ggplot并向下滚动帮助页面。底部是指向ggplot包索引的链接。滚动索引直到找到geom选项。
3.13.5 绘制来自2个以上单元格的数据
到目前为止，我们一直在考虑数据框中2个细胞的基因计数。但实际上我们的数据框中有10个细胞，最好将它们全部进行比较。如果我们想同时绘制所有10个细胞的数据怎么办？
目前我们还不能这样做，因为我们将每个细胞视为一个变量，并将该变量分配给x或y轴。我们可以创建一个10维图来绘制所有10个细胞的数据，但这是a)无法用ggplot来实现的，并且b)不太容易解释。我们可以做的是整理数据，以便我们有一个代表细胞ID的变量和另一个代表基因计数的变量，并将它们相互绘制在一起。在代码中，这看起来像：

> counts<-gather(counts, colnames(counts)[2:11], key = 'Cell_ID', value='Counts')
> head(counts)
  Gene_ids Cell_ID Counts
1    gene1   cell1      8
2    gene2   cell1     10
3    gene3   cell1      7
4    gene4   cell1     11
5    gene5   cell1     14
6    gene6   cell1     12

本质上，之前的问题是由于一个变量（Cell_ID）分散在多个列中，导致数据不整齐。现在我们已经解决了这个问题，我们可以更轻松地在一张图上绘制所有10个细胞的数据。

> ggplot(counts,aes(x=Cell_ID, y=Counts)) + geom_boxplot()

任务3：使用更新后的counts数据框绘制条形图，其中Cell_ID为x变量，Counts为y变量。提示：您可能会发现阅读?geom_bar很有帮助。
任务4：使用更新后的counts数据框绘制散点图，其中Gene_ids作为x变量，Counts作为y变量。
3.13.6 绘制热图
可视化基因表达数据的常用方法是使用热图。在这里，我们将使用R包pheatmap对一些我们称之为test的基因表达数据进行分析。

> library(pheatmap)
> set.seed(2)
> test = matrix(rnorm(200), 20, 10)
> test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
> test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
> test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
> colnames(test) = paste("Cell", 1:10, sep = "")
> rownames(test) = paste("Gene", 1:20, sep = "")
> pheatmap(test)

让我们花点时间弄清楚这张图向我们展示了什么。每一行代表一个基因，每一列代表一个细胞。每个基因在每个细胞中的表达程度由相应框的颜色表示。例如，我们可以从该图中看出Gene18在Cell10中表达较高，但在Cell1中表达较低。
该图还为我们提供了有关聚类算法结果的信息。一般来说，聚类算法旨在将数据点（例如细胞）分成几组，这些组的成员彼此之间的相似性要高于它们与其他数据点的相似性。图表顶部和左侧绘制的树是聚类算法的结果，使我们能够看到，例如，Cell4、8、2、6和10彼此之间的相似性，比Cell7、3、5、1和9之间的相似性要高。图表左侧的树表示对我们数据集中的基因应用聚类算法的结果。
如果我们仔细观察这些树，我们就会发现它们最终具有与细胞和基因相同数量的分支。也就是说，细胞簇的总数与细胞的总数相同，基因簇的总数与基因的总数相同。显然，这并没有提供太多信息，而且当我们查看超过10个细胞和20个基因时，这种方法就不切实际了。幸运的是，我们可以设置我们在图上看到的簇数。让我们尝试将基因簇数设置为2：

> pheatmap(test, kmeans_k = 2)

现在我们可以看到基因分为两类——相对于其他细胞，在Cell2、10、6、4 和8中上调的8个基因簇以及相对于其他细胞，在Cell2、10、6、4和8中下调的12个基因簇。
任务5：尝试将聚类数设置为3。您认为哪个聚类数更具信息量？
3.13.7 主成分分析
主成分分析（PCA）是一种统计过程，它使用变换将一组观测值转换为一组称为主成分的线性不相关变量的值。进行变换的目的是使第一个主成分尽可能多地解释数据中的变异性，并且每个后续主成分在必须与前一个成分正交的约束下，解释最大的可能方差。
PCA图是概览数据的好方法，有时可以帮助识别解释数据中大量变异的混杂因素。我们将在未来更深入地研究如何在单细胞RNA-seq分析中使用PCA图，这里的目的是让您了解PCA图是什么以及它们是如何生成的。
让我们为测试数据制作一个PCA图。我们可以使用ggfortify包让ggplot知道如何解释主成分。

> library(ggfortify)
> Principal_Components<-prcomp(test)
> autoplot(Principal_Components, label=TRUE)

任务6：将您的聚类与pheatmap聚类进行比较。它们有关联吗？（提示：查看我们绘制的第一个 pheatmap 的基因树）
任务7：制作counts的热图和PCA图：

> set.seed(1)
> counts <- as.data.frame(matrix(rpois(100, lambda = 10), ncol=10, nrow=10))
> rownames(counts) <- paste("gene", 1:10, sep = "")
> colnames(counts) <- paste("cell", 1:10, sep = "")

往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——3. 单细胞分析中的R/Bioconductor简介（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——3. 单细胞分析中的R/Bioconductor简介（2）
重生之我在剑桥大学学习单细胞RNA-seq分析——3. 单细胞分析中的R/Bioconductor简介（3）