参考资料：
https://www.yunbios.net/phyloseq.html
https://joey711.github.io/phyloseq/import-data.html
phyloseq 包，集OTU 数据导入，存储，分析和图形可视化于一体。不但利用了 R 中许多经典的工具进行生态学和系统发育分析（例如：vegan，ade4，ape， picante），同时还结合 ggplot2 以轻松生成发表级别的可视化结果。phyloseq 使用的S4类将一个研究所有相关的测序数据及元数据存储为单个对象，从而更容易共享数据并重复结果。
# 安装，已经安装的忽略。
> if (!requireNamespace("BiocManager", quietly = TRUE))
+ install.packages("BiocManager")
> BiocManager::install("phyloseq", version = "3.8")
# 在包的安装过程中出现了一个错误信息Error in install.packages : cannot remove prior installation of package ‘backports’
# 解决方案：关闭Rstudio，在R包安装路径中找到backports文件夹，全部删除，重启Rstudio重新安装就可以了
```
一个 phyloseq 类，通常由以下几个部分组成：
- `otu_table` ：一个数字矩阵 `matrix`，包含了 OTU 在每个样本中的丰度信息；
- `sample_data` ：一个`data.frame`，包含了所有样本的表型信息，行名必须匹配`otu_table` 中的样本名；
- `tax_table` ：一个字符矩阵 `matrix`，包含了 OTU 的物种信息，行名必须匹配`otu_table` 中的 OTU 名。
```R
# 导入out_table和tax_table，数据一定要转化成矩阵形式，不然后面使用phyloseq包的函数的时候会报错。
> rm(list = ls())
> otu_p <- as.matrix(read.table("otu.csv", header = TRUE, sep = ",", row.names = 1))
> tax_p <- as.matrix(read.table("taxon.csv", header = TRUE, sep = ",", row.names = 1))
# 将数据组合成一个 phyloseq 对象
> library("phyloseq")
> OTU = otu_table(otu_p, taxa_are_rows = TRUE)
> TAX = tax_table(tax_p)
> physeq_p = phyloseq(OTU, TAX)
> physeq_p
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 9702 taxa and 18 samples ]
## tax_table() Taxonomy Table: [ 9702 taxa by 7 taxonomic ranks ]
> plot_bar(physeq_p, fill = "Phylum")
# 导入sample_data
> sampledata = sample_data(data.frame(
+ Times = c(rep("BL", times=6), rep("Y10", times=6), rep("Y4", times=6)),
+ Depth = c(rep(c("top", "top", "top", "sub", "sub", "sub"), times = 3)),
+ row.names = sample_names(physeq_p),
+ stringsAsFactors = FALSE
+))
# 使用ape包建立 OTU 系统发育树并导入 phyloseq 对象
> library("ape")
> random_tree = rtree(ntaxa(physeq_p), rooted=TRUE, tip.label=taxa_names(physeq_p))
> plot(random_tree)
# 太密集了，数据太多，后期可能需要筛选一些OTU
# 现在，我们有了`otu_table`，`sample_data`，`tax_table`，`phy_tree`这四类数据，以下两种方法都可以将他们合并为一个 phyloseq 对象：1. 使用`merge_phyloseq`函数在之前创建的`physeq`对象中加入`sample_data`和`phy_tree`数据；2. 使用`physeq`函数重新创建一个`physeq`对象。（两种方式没有区别，我选择用第二种）
> physeq = phyloseq(OTU, TAX, sampledata, random_tree)
> physeq
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 9702 taxa and 18 samples ]
## sample_data() Sample Data: [ 18 samples by 2 sample variables ]
## tax_table() Taxonomy Table: [ 9702 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 9702 tips and 9701 internal nodes ]
# 数据的简单可视化，数据量太大没有跑
> plot_tree(physeq, color="Times", label.tips="taxa_names", ladderize="left", plot.margin=0.3)
> plot_tree(physeq, color="Depth", shape="Times", label.tips="taxa_names", ladderize="right", plot.margin=0.3)
> plot_heatmap(physeq)
> plot_heatmap(physeq, taxa.label="Phylum")
```
