基于phyloseq的排序分析

wentao

2019年5月16日

载入本次分析需要的包

phyloseq的排序函数为我们的排序分析进行了极大的简化，保证计算的同时兼顾可视化工作。

library("phyloseq"); packageVersion("phyloseq")## [1] '1.26.1'data(GlobalPatterns)
library("ggplot2"); packageVersion("ggplot2")## [1] '3.1.0'library("cluster"); packageVersion("cluster")## [1] '2.0.7.1'library("plyr"); packageVersion("plyr")## [1] '1.8.4'

设置主题 准备数据

theme_set(theme_bw())

数据得到 并进行过滤

count数量在一半样品中出现抵低于5的OTU

GP = GlobalPatterns
wh0 = genefilter_sample(GP, filterfun_sample(function(x) x > 5), A=0.5*nsamples(GP))
GP1 = prune_taxa(wh0, GP)

相对丰度标准化，这里使用10的六次方

GP1 = transform_sample_counts(GP1, function(x) 1E6 * x/sum(x))

仅仅保留主要的五个门类的OTU

phylum.sum = tapply(taxa_sums(GP1), tax_table(GP1)[, "Phylum"], sum, na.rm=TRUE)
top5phyla = names(sort(phylum.sum, TRUE))[1:5]
GP1 = prune_taxa((tax_table(GP1)[, "Phylum"] %in% top5phyla), GP1)

get_variable 函数返回mapping文件中的目标列

get_variable函数结合需要判断的变量，将分组分为两类，并添加到mapping文件中

human = get_variable(GP1, "SampleType") %in% c("Feces", "Mock", "Skin", "Tongue")
sample_data(GP1)$human <- factor(human)

phyloaeq封装的排序脚本使用

ordinate函数，我们只需要定义使用的排序方法：method，和使用的距离参数 distance,下面是基于brayZ距离的NMDS排序

下面展示loading矩阵信息

GP.ord <- ordinate(GP1, "NMDS", "bray")## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1333468
## Run 1 stress 0.1333468
## ... Procrustes: rmse 4.32893e-06 max resid 7.783055e-06
## ... Similar to previous best
## Run 2 stress 0.1637672
## Run 3 stress 0.1844742
## Run 4 stress 0.1529686
## Run 5 stress 0.1498725
## Run 6 stress 0.1469145
## Run 7 stress 0.1505794
## Run 8 stress 0.1891198
## Run 9 stress 0.1452642
## Run 10 stress 0.1488978
## Run 11 stress 0.1469145
## Run 12 stress 0.1685942
## Run 13 stress 0.1333468
## ... Procrustes: rmse 5.818375e-06 max resid 1.112148e-05
## ... Similar to previous best
## Run 14 stress 0.1675647
## Run 15 stress 0.171298
## Run 16 stress 0.1385323
## Run 17 stress 0.1710784
## Run 18 stress 0.1385323
## Run 19 stress 0.1498724
## Run 20 stress 0.1385323
## *** Solution reachedp1 = plot_ordination(GP1, GP.ord, type="taxa", color="Phylum", title="taxa")
print(p1)

### 可选分面loading展示，更加清晰直观p1 + facet_wrap(~Phylum, 3)

### 对样品进行排序

type 选择samples 对样品进行排序，这也是我们使用更多的一种方法。

p2 = plot_ordination(GP1, GP.ord, type="samples", color="SampleType", shape="human")
p2 + geom_polygon(aes(fill=SampleType)) + geom_point(size=5) + ggtitle("samples")

### biplot 选项将样品和loading一同展示

但是重叠很严重

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)

### split 通过将两者分开分，面展示解决问题p4 = plot_ordination(GP1, GP.ord, type="split", color="Phylum", shape="human", label="SampleType", title="split")
p4

### 自定义颜色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"
p4 + scale_color_manual(values=p4cols)

### 下面使用多种排序方法进行排序，并比对

使用llply函数

dist = "bray"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
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)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1333468
## Run 1 stress 0.1528906
## Run 2 stress 0.1827625
## Run 3 stress 0.1333468
## ... New best solution
## ... Procrustes: rmse 3.155156e-06 max resid 7.816854e-06
## ... Similar to previous best
## Run 4 stress 0.1468477
## Run 5 stress 0.1333468
## ... New best solution
## ... Procrustes: rmse 2.975489e-06 max resid 7.840301e-06
## ... Similar to previous best
## Run 6 stress 0.1681228
## Run 7 stress 0.16608
## Run 8 stress 0.1521722
## Run 9 stress 0.1385323
## Run 10 stress 0.1570722
## Run 11 stress 0.1455574
## Run 12 stress 0.1510167
## Run 13 stress 0.1653386
## Run 14 stress 0.1658681
## Run 15 stress 0.1385327
## Run 16 stress 0.1507161
## Run 17 stress 0.169588
## Run 18 stress 0.1469145
## Run 19 stress 0.180481
## Run 20 stress 0.1518735
## *** Solution reached

提取作图所需数据

names(plist) <- ord_meths
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"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

### 选择其中一个展示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

### 使用unifrac距离进行PCoa分析ordu = ordinate(GP1, "PCoA", "unifrac", weighted=TRUE)
plot_ordination(GP1, ordu, color="SampleType", shape="human")

展示排序结果

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")

当我们将排序做完之后

这里仅仅做了排序的分析，其实我们的排序需要进行检验。