phyloseq: Explore microbiome pro

本节主要是在Phyloseq[Tutorials(http://joey711.github.io/phyloseq/gap-statistic.html)学习Ordination Plots的功能

---

plot_ordination函数在很大程度上依赖于distance和ordiante函数。

---

此外，phyloseq包还包含一个convenience function用于从一个排序中的大量点集合中取子集，称为subset_ord_plot.其中也有一个单独的subset_ord_plot_tutorial教程，了解更详细的细节和例子。

Load packages,Prepare Data

library("phyloseq")
data(Globalpatterns)
library("ggplot2)
library("pyr")
theme_set(theme_bw())

我们希望从这些数据中过滤低发生率、低代表性的otu，因为对于本教程的目的来说，它们本质上是噪音变量。在实践中，您可能应该执行并清楚地记录合理的预处理步骤，这些步骤在phyloseq包中提供了示例和关于dedicated preprocessing tutorial。预处理对于以图形方式显示数据中的高级模式特别有用，以及创建在短时间内计算的示例。你的推理和决定在预处理是非常重要的，并取决于你。我在这里使用了几种不同的预处理方法来进行说明，因为数据简化的程度对我们目的很有用。我不能断言这些是您的数据和研究目标的“最佳”方法,但是，我强烈建议你们认真思考你们所做的任何预处理,将其完整地记录下来，只有在您能够决定这些选择并检查它们是否可靠的情况下，才决定将其包含在您的最终分析管道中。为了快速演示和比较不同排序方法的结果，我将首先进一步过滤/预处理otu在GP1.这种预处理的重点将是限制otu的数量，而不是保护数据中的固有模式.

#删除在超过一半的样本中没有显示超过5次的otu,genefilter_sample()用任意函数过滤OTUs
GP = GlobalPatterns
#genefilter_sample(),genefilter_sample(X, flist, A=1)
wh0=genefilter_sample(GP,filterfun_sample(function(x) x > 5), A=0.5*nsamples(GP))
GP1 = prune_taxa(wh0, GP)
#转换为均匀采样深度,函transform_sample_counts(),逐个样本地转换otu_table中的丰富数据,transform_sample_counts(physeq, fun, ...)
#只保留丰度最高的五个门类，tapply()函数，对一个不规则数组的每个单元格应用一个函数tapply(X, INDEX, FUN = NULL, ..., default = NA, simplify = TRUE)
phylum.sum=tapply(taxa_sum(GP1),taxa_table(GP1)[,"Phylum"],sum,na.rm=TRUE)
top5phyla=names(sort(phylum.sum,TRUE))[1:5]
GP1=prune_taxa((taxa_table(GP1)[,"Phylum"]%in%top5phylum),GP1)
#定义一个人相关的和非人相关的分类变量
human = get_variable(GP1, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
sample_data(GP1)$human <- factor(human)

####Four main ordination plots
`plot_ordination`函数支持一个排序的四种基本表示.
#####(1)Just OTUs
让我们从绘制otu开始，并通过门水平来为这些点着色.
注意，即使在我们的“修剪”数据集中，也有ntaxa(GP1)= 204 OTUs
GP.ord <- ordinate(GP1, "NMDS", "bray")
p1 = plot_ordination(GP1, GP.ord, type="taxa", color="Phylum", title="taxa")
print(p1)
#facet_wrap(),将一维的带状面板包装成二维
p1+facet_wrap(~Phylum,3)

(2)Just samples

接下来，我们只绘制样本，并使用“SampleType”对这些点进行着色，同时根据它们是否与人类相关来修改形状,还添加了一些额外的ggplot2层，以使图片更加美观

p2 = plot_ordination(GP1, GP.ord, type="samples", color="SampleType", shape="human")
#ggplot2的美化，添加点间的连线geom_polygon
p2 + geom_polygon(aes(fill=SampleType)) + geom_point(size=5) + ggtitle("samples")

(3)biplot graphic

plot_ordination函数还可以自动创建两个不同的图形布局，其中样品和OTUs一起被绘制在"bioplot"中。请注意，这需要本质上不是仅限采样的排序的方法，例如，这不适用于UniFrac/PCoA

#不同的颜色代表不同的样品，不同的形状代表不同的门类
p3 = plot_ordination(GP1, GP.ord, type="biplot", color="SampleType", shape="Phylum", title="biplot")
#修改了一些自动形状缩放的东西
p3 = plot_ordination(GP1, GP.ord, type="biplot", color="SampleType", shape="Phylum", title="biplot")
# Some stuff to modify the automatic shape scale
GP1.shape.names = get_taxa_unique(GP1, "Phylum")
GP1.shape <- 15:(15 + length(GP1.shape.names) - 1)
names(GP1.shape) <- GP1.shape.names
GP1.shape["samples"] <- 16
p3 + scale_shape_manual(values=GP1.shape)

(4)split graphic

在前面的图形中，遮挡问题相当严重,在这种情况下type="split"选项可能会有帮助,其中样品/ otu被分开在两个并排的面板上…

p4=plot_ordination(GP1, GP.ord, type="split", color="Phylum", shape="human",label="SampleType", title="split")
p4
#下面的函数重现ggplot2的默认颜色比例
gg_color_hue<-function(n){
    hues = seq(15, 375, length=n+1)
    hcl(h=hues, l=65, c=100)[1:n]
}
color.names <- levels(p4$data$Phylum)
p4cols <- gg_color_hue(length(color.names))
names(p4cols) <- color.names
p4cols["samples"] <- "black"
#scale_colour_manual访问指定的颜色,[R语言中颜色使用方法参考汇总](https://blog.csdn.net/adan_journal_of_data/article/details/79242721)
p4 + scale_color_manual(values=p4cols)

Supported Ordination Methods

在本节中，我将通过不同的方法参数选项循环到plot_ordination函数，将绘图结果存储在一个列表中，然后使用ggplot2将这些结果绘制到一个组合图形中.

dist = "bray"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
#llply拆分列表、应用函数并在列表中返回结果
plist = llply(as.list(ord_meths), function(i, physeq, dist){
        ordi = ordinate(physeq, method=i, distance=dist)
        plot_ordination(physeq, ordi, "samples", color="SampleType")
}, GP1, dist)
#重命名plist对象
names(plist) <- ord_meths

前面的代码块执行每个排序方法，根据每个排序结果的前两个轴创建相应的图形，然后将每个ggplot2绘图对象存储在名为plist的列表的不同命名元素中。下面的数据块将从每个单独的图中提取数据，并将它们重新放在一个大数据中。

#ldply拆分列表，应用函数，并在数据框中返回结果,其中llply中的第二字母代表返回list
pdataframe = ldply(plist, function(x){
    df = x$data[, 1:2]
    colnames(df) = c("Axis_1", "Axis_2")
    return(cbind(df, x$data))
})
names(pdataframe)[1] = "method"

现在所有的排序结果都合并在data.frame中了，称为pdataframe,我们可以使用它来制作一个标准的分面ggplot散点图

p = ggplot(pdataframe, aes(Axis_1, Axis_2, color=SampleType, shape=human, fill=SampleType))
p = p + geom_point(size=4) + geom_polygon()
p = p + facet_wrap(~method, scales="free")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + scale_colour_brewer(type="qual", palette="Set1")
p

如果您想要重绘单个plot的更大版本，可以通过打印原始plist (pdataframe就是从这个plist生成的)来实现。plist的每个元素已经是一个ggplot2图形。例如，我们可以通过打印列表的第二个元素来重新绘制(DCA)plist[[2]]
现在添加一些额外的层，使它看起来更好

p=plist[[2]] + scale_colour_brewer(type="qual", palette="Set1")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + geom_point(size=5) + geom_polygon(aes(fill=SampleType))
p

MDS (“PCoA”) on Unifrac Distances

使用ordinate函数，同时执行weightd UniFrac，然后对距离矩阵(第一行)进行主坐标分析.接下来将数据和排序结果传递给plot_ordination使用缺省ggplot2设置创建ggplot2输出图形

ordu=ordinate(GP1,"PCoA", "unifrac", weighted=TRUE)
plot_ordination(GP1, ordu, color="SampleType", shape="human")

现在使用一些额外的ggplot2层使图形看起来更好

p=plot_ordination(GP1,ordu,color="SampleType",shape="human")
p=p+geom_point(size=7,alpha=0.75)
p = p + scale_colour_brewer(type="qual", palette="Set1")
p + ggtitle("MDS/PCoA on weighted-UniFrac distance, GlobalPatterns")