phyloseq | 用 R 分析微生物组数据及可视化（二）

上期回顾

上一篇中我简单介绍了phyloseq包的功能与数据导入，在这一期中将详细介绍phyloseq包中的各种绘图函数。

首先载入phyloseq和ggplot：

library("phyloseq")
library("ggplot2")
# 设置 ggplot2 主题
theme_set(theme_bw())

在这次的教程中以phyloseq自带的GlobalPatterns数据集和enterotype数据集作为例子。GlobalPatterns数据来自一篇2011年的 PNAS 文章 （Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample），比较了 25 个环境样本和 3 个已知的“模拟微生物群落”，包含了 9 种样本类型。enterotype数据集即人类肠型数据集来自一篇2011 年的 Nature 文章（Enterotypes of the human gut microbiome），该文章使用宏基因组鸟枪法测序比较了22名受试者的粪便微生物群落，作者进一步将这些微生物群落与来自其他研究的受试者的粪便群落进行比较。

物种多样性分析

导入GlobalPatterns数据集：

data(GlobalPatterns)

首先对这些数据做一些预处理：

# prune OTUs that are not present in at least one sample
GP <- prune_taxa(taxa_sums(GlobalPatterns) > 0, GlobalPatterns)
# Define a human-associated versus non-human categorical variable:
human <- get_variable(GP, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
# Add new human variable to sample data:
sample_data(GP)$human <- factor(human)

我们可以使用plot_richness函数轻松创建一个复杂的图形，用于比较GlobalPatterns数据集中不同环境类型的样本的物种多样性。

alpha_meas = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson")
(p <- plot_richness(GP, "human", "SampleType", measures=alpha_meas))

还可以往上面加ggplot图层：

p + geom_boxplot(data=p$data, aes(x=human, y=value, color=NULL), alpha=0.1)

以上就是GlobalPatterns数据集中样本的 Alpha 多样性分析。每幅图都展示了不同的物种多样性指数，不同的颜色表示不同的样本类型。在每个小组中，将样本进一步分为与人类相关的（TRUE）或不相关的（FALSE），并且在这两个组的顶部叠加箱图，说明这些与人类相关样本不如环境样本物种丰富。

绘制进化树

phyloseq还可以很容易地绘制带注释的系统发育树，并加入特定分类的信息，以及观察到的个体数量。

这一节中，我们先取只包含Chlamydiae门的数据：

GP.chl <- subset_taxa(GP, Phylum=="Chlamydiae")

现在我们将创建这个GlobalPatterns数据集的子集的进化树并用SampleType变量进行着色，该变量表示微生物组样本的来源。以下命令还可以标记每个分类中每个样本（如果有的话）中观察到的菌群数量，菌群丰度越高，形状就越大。

plot_tree(GP.chl, color="SampleType", shape="Family", label.tips="Genus", size="Abundance")

条形图

导入enterotype数据集：

data(enterotype)

在enterotype数据集中，可用的数据是菌群的相对丰度，而不是原始的数量，所以就不能用来做物种多样性分析。而且 OTU 仅有属级别的注释信息。

我们先从一个简单的秩丰度条形图开始，计算数据集中每个OTU的累积分数丰度。在这个条形图中，我们通过样本总数（280）进一步标准化。

par(mar = c(10, 4, 4, 2) + 0.1) # make more room on bottom margin
N <- 30
barplot(sort(taxa_sums(enterotype), TRUE)[1:N]/nsamples(enterotype), las=2)

上面的例子，我们用到了R最基础的barplot函数和phyloseq提供的taxa_sums和nsamples函数。你可以看到大约在排名第十的菌属后菌群的相对丰度出现了剧烈下滑。因此我们可以尝试在最多的十个菌属中来分析肠型数据。

TopNOTUs <- names(sort(taxa_sums(enterotype), TRUE)[1:10]) 
ent10   <- prune_taxa(TopNOTUs, enterotype)
print(ent10)

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 280 samples ]
## sample_data() Sample Data:       [ 280 samples by 9 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 1 taxonomic ranks ]

这个数据集中共有 280 个样本。接着，看看原始的metadata中包含了哪些表型信息：

sample_variables(ent10)

## [1] "Enterotype"     "Sample_ID"      "SeqTech"        "SampleID"      
## [5] "Project"        "Nationality"    "Gender"         "Age"           
## [9] "ClinicalStatus"

在下面的例子中，我们根据每个排名前十的菌属来分面（facet），在每个子图中以三种不同的测序方法分成三组，并以颜色表示每个菌属来源的样品的肠型标签。可以发现，不同的测序技术对检测到哪个属和菌属的丰度水平具有很大的影响。

plot_bar(ent10, "SeqTech", fill="Enterotype", facet_grid=~Genus)

从上图还可以看出，肠型1 中有较高丰度的 Bacteroides， 肠型2 中有较高丰度的 Prevotella；而对于肠型3 可以观察到高丰度的 Blautia 仅在454-焦磷酸测序的数据中出现，在 Illumina 或 Sanger 的数据中丰度都很低，说明高丰度的 Blautia 可能是引入的误差。

热图

下面以GlobalPatterns数据集的子集（仅包含Crenarchaeota门的数据）来展示热图的功能：

data("GlobalPatterns")
gpac <- subset_taxa(GlobalPatterns, Phylum=="Crenarchaeota")
(p <- plot_heatmap(gpac, "NMDS", "bray", "SampleType", "Family"))

如果需要自定义X轴和Y轴的标签，只要这样就好啦：

p$scales$scales[[1]]$name <- "My X-Axis"
p$scales$scales[[2]]$name <- "My Y-Axis"
print(p)

除此之外，你还可以很方便地选择不同的距离矩阵（例子中使用了Bray Curtis距离矩阵）来进行排序，是不是很简单~

微生物网络图

这一节中，我们使用enterotype数据集，来创建网络图：

data(enterotype)
plot_net(enterotype, maxdist=0.4, color="SeqTech", shape="Enterotype")

Reference

• http://www.bioconductor.org/packages/release/bioc/vignettes/phyloseq/inst/doc/phyloseq-analysis.html

---

猜你喜欢

还有更多文章，请移步公众号阅读

▼ 如果你生信基本技能已经入门，需要提高自己，请关注下面的生信技能树，看我们是如何完善生信技能，成为一个生信全栈工程师。

▼ 如果你是初学者，请关注下面的生信菜鸟团，了解生信基础名词，概念，扎实的打好基础，争取早日入门。