phyloseq 包，集OTU 数据导入，存储，分析和图形可视化于一体。不但利用了 R 中许多经典的工具进行生态学和系统发育分析（例如：vegan，ade4，ape， picante），同时还结合 ggplot2 以轻松生成发表级别的可视化结果。phyloseq 使用的S4类将一个研究所有相关的测序数据及元数据存储为单个对象，从而更容易共享数据并重复结果。
- GitHub 地址：https://github.com/joey711/phyloseq
安装
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("phyloseq", version = "3.8")
数据导入
一个 phyloseq 类，通常由以下几个部分组成：
-
otu_table：一个数字矩阵matrix，包含了 OTU 在每个样本中的丰度信息； -
sample_data：一个data.frame，包含了所有样本的表型信息，行名必须匹配otu_table中的样本名； -
tax_table：一个字符矩阵matrix，包含了 OTU 的物种信息，行名必须匹配otu_table中的 OTU 名。
1.jpg
在下面的例子中，我用 R 随机创造了一些模拟数据，要是有真实数据的话，也可以直接导入。
# Create a pretend OTU table that you read from a file, called otumat
otumat = matrix(sample(1:100, 100, replace = TRUE), nrow = 10, ncol = 10)
rownames(otumat) <- paste0("OTU", 1:nrow(otumat))
colnames(otumat) <- paste0("Sample", 1:ncol(otumat))
# Create a pretend taxonomy table
taxmat = matrix(sample(letters, 70, replace = TRUE), nrow = nrow(otumat), ncol = 7)
rownames(taxmat) <- rownames(otumat)
colnames(taxmat) <- c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species")
接下来，我们需要将他们组合成一个 phyloseq 对象：
library("phyloseq")
OTU = otu_table(otumat, taxa_are_rows = TRUE)
TAX = tax_table(taxmat)
physeq = phyloseq(OTU, TAX)
plot_bar(physeq, fill = "Family")
image.png
下面导入样本的表型信息：
> # Create random sample data
> sampledata = sample_data(data.frame(
+ Location = sample(LETTERS[1:4], size=nsamples(physeq), replace=TRUE),
+ Depth = sample(50:1000, size=nsamples(physeq), replace=TRUE),
+ row.names=sample_names(physeq),
+ stringsAsFactors=FALSE
+ ))
> sampledata
Sample Data: [10 samples by 2 sample variables]:
Location Depth
Sample1 C 722
Sample2 C 409
Sample3 A 154
Sample4 C 705
Sample5 B 464
Sample6 C 815
Sample7 D 841
Sample8 C 681
Sample9 A 773
Sample10 D 205
使用ape包建立 OTU 系统发育树并导入 phyloseq 对象：
library("ape")
random_tree = rtree(ntaxa(physeq), rooted=TRUE, tip.label=taxa_names(physeq))
plot(random_tree)
3.jpg
现在，我们有了otu_table，sample_data，tax_table，phy_tree这四类数据，以下两种方法都可以将他们合并为一个 phyloseq 对象：1. 使用merge_phyloseq函数在之前创建的physeq对象中加入sample_data和phy_tree数据；2. 使用physeq函数重新创建一个physeq对象。
physeq1 = merge_phyloseq(physeq, sampledata, random_tree)
使用physeq函数重新创建一个physeq对象：
physeq2 = phyloseq(OTU, TAX, sampledata, random_tree)
看看用这两种方法创建的physeq对象是否一样：
identical(physeq1, physeq2)
## [1] TRUE
尝试对这些数据进行一些可视化：
plot_tree(physeq1, color="Location", label.tips="taxa_names", ladderize="left", plot.margin=0.3)
image.png
plot_tree(physeq1, color="Depth", shape="Location", label.tips="taxa_names", ladderize="right", plot.margin=0.3)
image.png
再做一些热图：
plot_heatmap(physeq1)
image.png
plot_heatmap(physeq1, taxa.label="Phylum")
7.jpg
Reference
phyloseq | 用 ******R 分析微*********生物组数据及可视化(二)***
Phyloseq（一）简单介绍了phyloseq包的功能与数据导入，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))
image.png
还可以往上面加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")
image.png
条形图
导入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)
image.png
上面的例子，我们用到了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" "Project"
[6] "Nationality" "Gender" "Age" "ClinicalStatus"
在下面的例子中，我们根据每个排名前十的菌属来分面（facet），在每个子图中以三种不同的测序方法分成三组，并以颜色表示每个菌属来源的样品的肠型标签。可以发现，不同的测序技术对检测到哪个属和菌属的丰度水平具有很大的影响。
plot_bar(ent10, "SeqTech", fill="Enterotype", facet_grid=~Genus)
image.png
从上图还可以看出，肠型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"))
image.png
如果需要自定义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")
image.png
Reference
- http://www.bioconductor.org/packages/release/bioc/vignettes/phyloseq/inst/doc/phyloseq-analysis.html
phyloseq | 用 R 分析微生物组数据及可视化(三）
Phyloseq（三）主要介绍如何使用phyloseq包进行排序分析（Ordination analysis）。
排序分析
排序分析（Ordination analysis）是探索复杂的系统发育测序数据的得力工具。简而言之，排序（ordination）的过程就是在一个可视化的低维空间重新排列这些样方，使得样方之间的距离最大程度地反映出平面散点图内样方之间的关系信息。phyloseq分别通过ordinate()和plot_ordination()函数，来进行排序分析及其可视化。虽然排序分析有很多方法和参数，但最基本的代码就是这样：
my.physeq <- import("Biom", BIOMfilename="myBiomFile.biom")my.ord <- ordinate(my.physeq)plot_ordination(my.physeq, my.ord, color="myFavoriteVarible")
PCoA 主坐标分析
PCoA（principal co-ordinates analysis）主坐标分析法是一种与PCA 类似的降维排序方法。PCoA与PCA的区别在于PCA是基于原始的物种组成矩阵所做的分析，仅仅比较的是物种丰度的不同，而PCoA首先根据不同的距离算法计算样品之间的距离，然后对距离矩阵进行处理，使图中点间的距离正好等于原来的差异数据，实现定性数据的定量转换。
下面以 Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample 这篇文章的Figure 5为例，作者展示了 unweighted-UniFrac 距离矩阵 PCoA 的前三个轴。
15.jpg
下面将在完整的数据集上再现unweighted-UniFrac距离矩阵的计算。由于OTU数量众多，所以可能需要很长的时间，建议对大型数据集进行并行化运算，有关并行化的详细信息，请参阅UniFrac()的文档及示例。
data(GlobalPatterns)# GPUF <- UniFrac(GlobalPatterns)# Load the pre-computed distance matrix, GPUFload(system.file("doc", "Unweighted_UniFrac.RData", package="phyloseq"))# Calculate the PCoA on this distance matrix, GPUFGloPa.pcoa = ordinate(GlobalPatterns, method="PCoA", distance=GPUF)
在查看结果之前，我们可以用碎石图来查看排名前几的轴的贡献率：
plot_scree(GloPa.pcoa, "Scree plot for Global Patterns, UniFrac/PCoA")
16.jpg
通过碎石图，可以看出前三个轴代表距离总变化的43％，第四轴代表另外9％，因此也可能需要进行探索。碎石图是任何排序方法的重要工具，因为轴的相对贡献率在不同的数据集之间可能变化很大。
下面开始重复Figure 5，但我们将原来三维的图，转化为两张图：
(p12 <- plot_ordination(GlobalPatterns, GloPa.pcoa, "samples", color="SampleType") + geom_point(size=5) + geom_path() + scale_colour_hue(guide = FALSE) )
17.jpg
(p13 <- plot_ordination(GlobalPatterns, GloPa.pcoa, "samples", axes=c(1, 3), color="SampleType") + geom_line() + geom_point(size=5) )
18.jpg
NMDS 非度量多维尺度分析
非度量多维尺度分析是一种将多维空间的研究对象（样本或变量）简化到低维空间进行定位、分析和归类，同时又保留对象间原始关系的数据分析方法。适用于无法获得研究对象间精确的相似性或相异性数据，仅能得到他们之间等级关系数据的情形。其基本特征是将对象间的相似性或相异性数据看成点间距离的单调函数，在保持原始数据次序关系的基础上，用新的相同次序的数据列替换原始数据进行度量型多维尺度分析。
其特点是根据样品中包含的物种信息，以点的形式反映在多维空间上，而对不同样品间的差异程度，则是通过点与点间的距离体现的，最终获得样品的空间定位点图。
# (Re)load UniFrac distance matrix and GlobalPatterns datadata(GlobalPatterns)load(system.file("doc", "Unweighted_UniFrac.RData", package="phyloseq"))# perform NMDS, set to 2 axesGP.NMDS <- ordinate(GlobalPatterns, "NMDS", GPUF)
## Run 0 stress 0.1432774 ## Run 1 stress 0.1432625 ## ... New best solution## ... Procrustes: rmse 0.003499558 max resid 0.01265589 ## Run 2 stress 0.2467057 ## Run 3 stress 0.2234687 ## Run 4 stress 0.1432774 ## ... Procrustes: rmse 0.003483686 max resid 0.01260292 ## Run 5 stress 0.1432774 ## ... Procrustes: rmse 0.003487121 max resid 0.01261713 ## Run 6 stress 0.1856306 ## Run 7 stress 0.1432625 ## ... New best solution## ... Procrustes: rmse 1.405265e-05 max resid 3.151096e-05 ## ... Similar to previous best## Run 8 stress 0.2574296 ## Run 9 stress 0.1809282 ## Run 10 stress 0.1809285 ## Run 11 stress 0.1432625 ## ... New best solution## ... Procrustes: rmse 4.788746e-06 max resid 1.24343e-05 ## ... Similar to previous best## Run 12 stress 0.1432774 ## ... Procrustes: rmse 0.003470474 max resid 0.01255052 ## Run 13 stress 0.1432625 ## ... Procrustes: rmse 1.152937e-05 max resid 3.31644e-05 ## ... Similar to previous best## Run 14 stress 0.1432625 ## ... New best solution## ... Procrustes: rmse 2.462092e-06 max resid 6.050872e-06 ## ... Similar to previous best## Run 15 stress 0.1856315 ## Run 16 stress 0.1432625 ## ... Procrustes: rmse 9.826556e-06 max resid 3.260556e-05 ## ... Similar to previous best## Run 17 stress 0.1809284 ## Run 18 stress 0.1856306 ## Run 19 stress 0.1432774 ## ... Procrustes: rmse 0.003490316 max resid 0.01261738 ## Run 20 stress 0.1432774 ## ... Procrustes: rmse 0.003471665 max resid 0.01255283 ## *** Solution reached
(p <- plot_ordination(GlobalPatterns, GP.NMDS, "samples", color="SampleType") + geom_line() + geom_point(size=5) )
19.jpg
使用非度量多维尺度分析对GlobalPatterns数据集的样本之间的 unweighted-UniFrac 距离矩阵进行排序。每个样本点根据环境类型着色，如果它们属于同一类型，则由一条线连接。与GlobalPatterns文章中的图5进行比较。该图很好地显示了来自不同栖息地的微生物群落之间的相对差异。但是，它未能表明群落之间的差异。对于提供解释样本（或样本组）之间差异的群落信息的排序方法，我们使用对应分析（CA）。
CA 对应分析
在这一节中，我们使用 Correspondence Analysis 的排序方法的各种功能继续探索GlobalPatterns数据集。
让我们首先进行对应分析并查看碎石图。碎石图表明GlobalPatterns数据集具有相当高的维度，前两个 CA 轴占总变异的不到17％，但随着轴数的增加，贡献度并没有急剧下降，每个轴仅比前一个下降一点点。
首先，为了节约运行时间，我们只取GlobalPatterns数据集的子集：
data(GlobalPatterns)# Take a subset of the GP dataset, top 200 speciestopsp <- names(sort(taxa_sums(GlobalPatterns), TRUE)[1:200])GP <- prune_taxa(topsp, GlobalPatterns)# Subset further to top 5 phyla, among the top 200 OTUs.top5ph <- sort(tapply(taxa_sums(GP), tax_table(GP)[, "Phylum"], sum), decreasing=TRUE)[1:5]GP <- subset_taxa(GP, Phylum %in% names(top5ph))# 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)
进行 CA 对应分析：
# Now perform a unconstrained correspondence analysisgpca <- ordinate(GP, "CCA")# Scree plotplot_scree(gpca, "Scree Plot for Global Patterns Correspondence Analysis")
20.jpg
(p12 <- plot_ordination(GP, gpca, "samples", color="SampleType") + geom_line() + geom_point(size=5) )
21.jpg
(p34 <- plot_ordination(GP, gpca, "samples", axes=c(3, 4), color="SampleType") + geom_line() + geom_point(size=5) )
22.jpg
第一张图说明粪便和模拟数据紧密地聚集在一起，远离CA1上的其他样本。在CA2上，皮肤和舌头样本也类似地分开。总而言之，前两个轴已经解释了人类相关“环境”与数据集中其他非人类环境的分离，以及舌头和皮肤样本与粪便的分离。
我们现在将进一步研究这种数据结构，使我们能够比较单个分类群在相同图形空间上的贡献度：
p1 <- plot_ordination(GP, gpca, "species", color="Phylum")(p1 <- ggplot(p1$data, p1$mapping) + geom_point(size=5, alpha=0.5) + facet_wrap(~Phylum) + scale_colour_hue(guide = FALSE) )
23.jpg
我们还可以将物种点总结为 2D 密度图：
(p3 <- ggplot(p1$data, p1$mapping) + geom_density2d() + facet_wrap(~Phylum) + scale_colour_hue(guide = FALSE) )
24.jpg
上面的图表展示了一些有用的模式和有趣的异常值，但如果我们想要具体了解这些门在多大程度上影响每个轴呢？下面的代码使用了箱线图来展示这个结果。
library("reshape2")# Melt the species-data.frame, DF, to facet each CA axis separatelymdf <- melt(p1$data[, c("CA1", "CA2", "Phylum", "Family", "Genus")], id=c("Phylum", "Family", "Genus") )# Select some special outliers for labellingLF <- subset(mdf, variable=="CA2" & value < -1.0)# build plot: boxplot summaries of each CA-axis, with labelsp <- ggplot(mdf, aes(Phylum, value, color=Phylum)) + geom_boxplot() + facet_wrap(~variable, 2) + scale_colour_hue(guide = FALSE) + theme_bw() + theme( axis.text.x = element_text(angle = -90, vjust = 0.5) )# Add the text label layer, and render ggplot graphic(p <- p + geom_text(data=subset(LF, !is.na(Family)), mapping = aes(Phylum, value+0.1, color=Phylum, label=Family), vjust=0, size=2))
25.jpg
通过这种方法，更容易看到相对于其他门的异常聚集的特定物种，例如上图中的拟杆菌（Bacteroidetes）属的Prevotellaceae，在舌/皮肤样品中负CA2方向贡献最多。
另一种方法是直接展示物种丰度来说明菌群的差异变化，使用前面提及的plot_bar函数：
plot_bar(GP, x="human", fill="SampleType", facet_grid= ~ Phylum)
26.jpg
从上面这个图中也可以看出一些模式：
-
Cyanobacteria 和 Actinobacteria 在人类样品较少；
-
Firmicutes 在人类样本中较多；
-
Acidobacteria 和 Verrucomicrobia 在粪便样本中较多。
根据之前的研究，这些结果并不是非常令人惊讶，但希望它可以应用于你所拥有的其他数据集中。
DPCoA 双主坐标分析
下面的示例展示了使用了phyloseq中的plot_ordination()函数来进行双主坐标分析 ：
# Perform ordinationGP.dpcoa <- ordinate(GP, "DPCoA") # Generate default ordination bi-plotpdpcoa <- plot_ordination( physeq = GP, ordination = GP.dpcoa, type="biplot", color="SampleType", shape="Phylum")# Adjust the shape scale manually # to make taxa hollow and samples filled (advanced)shape.fac <- pdpcoa$data$Phylumman.shapes <- c(19, 21:25)names(man.shapes) <- c("Samples", levels(shape.fac)[levels(shape.fac)!="Samples"])p2dpcoa <- pdpcoa + scale_shape_manual(values=man.shapes)p2dpcoa
27.jpg
# Show just Samples or just Taxaplot_ordination(GP, GP.dpcoa, type="taxa", shape="Phylum")
28.jpg
plot_ordination(GP, GP.dpcoa, type="samples", color="SampleType")
29.jpg
# Splitplot_ordination(GP, GP.dpcoa, type="split", color="SampleType", shape="Phylum") + ggplot2::scale_colour_discrete()
30.jpg
聚类分析
层级聚类（例如，hclust）是可视化样本距离矩阵的另一种流行的方式。在下面的示例中，我们使用了 unweighted UniFrac 距离矩阵和 UPGMA 方法（hclust参数： method="average"）：
# (Re)load UniFrac distance matrix and GlobalPatterns datadata(GlobalPatterns)load(system.file("doc", "Unweighted_UniFrac.RData", package="phyloseq"))# Manually define color-shading vector based on sample type.colorScale <- rainbow(length(levels(get_variable(GlobalPatterns, "SampleType"))))cols <- colorScale[get_variable(GlobalPatterns, "SampleType")] GP.tip.labels <- as(get_variable(GlobalPatterns, "SampleType"), "character")# This is the actual hierarchical clustering call, specifying average-link clusteringGP.hclust <- hclust(GPUF, method="average")plot(GP.hclust, col=cols)
31.jpg
Reference
网友评论