phyloseq: Explore microbiome pro

最近看了关于介绍phyloseq包的这篇文章，在帮助文档界面上有从数据导入到各种分析的代码，对其代码进行重复、学习，从而更加熟练其代码，以期将来实战分析更得心应手。

1.数据的导入(Data import)

在加载library("phyloseq")后利用 ?import来获取导入数据的帮助文档

1.1 Loading included data

利用R语言中的data命令加载包中包含的预导入数据集data(GloalPatterns),为ggplot图形定义一个默认主题them_set(theme_bw())

1.2 phyloseq-ize Data already in R

在R会话中已经存在的任何数据都可以通过phyloseq的函数和方法进行调涌/识别载入,If you can get the data into R, then you can get it “into” phyloseq.phyloseq提供了构建phyloseq组成数据的工具，以及不同实验水平多组成的数据对象，即phyloseq-class，有以下函数组成,主要就是包含otu的表格导入，样品组成的导入，以及otu分类组成的导入.

?phyloseq
?otu_table
?sample_data
?tax_table

• otu_table-适用于任何数字 matrix,还必须指定物种是行还是列
• sample_data-适用于任何data.frame,如果您计划将它们组合为一个phyloseq-objec，则行名称必须与otu_table中的样品名称匹配
• tax_table--适用于任何字符 matrix,如果您计划将其与一个phyloseq对象组合，则行名称必须与otu_table的OTU名称(taxa_names)匹配
• phyloseq-
• merge-phyloseq-可以获取任意数量的phyloseq对象和/或phyloseq组分，并尝试将它们组合成一个更大的phyloseq对象。这对于将独立导入的component添加到已经创建的phyloseq对象非常有用

Example-在下面的示例中，作者在R中定义随机的示例数据表，然后将它们组合成一个phyloseq对象，使用基本的R代码创建典型的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)
#sample为抽样函数
sample(x, size, replace = FALSE, prob = NULL)

加入样品名和OTU表名称

#利用了paste0函数，连接向量后转换为字符，paste0(..., collapse = NULL)
rownames(otumat) <- paste0("OTU", 1:nrow(otumat))
colnames(otumat) <- paste0("Sample", 1:ncol(otumat))

新建一个taxonomy分类表

#letters为26个小写英文字母，LETTERS为26个大写的英文字母
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如何将上述它们组合成一个phyloseq-object

library("phyloseq")
OTU = otu_table(otumat, taxa_are_rows = TRUE)#注意otu表的行名
TAX = tax_table(taxmat)
#将otu表和tax进行组合
physeq = phyloseq(OTU, TAX)
#对不同样品的**Family**水平进行柱状图展示,plot_bar(physeq, x="Sample", y="Abundance", fill=NULL,title=NULL, facet_grid=NULL)
plot_bar(physeq, fill = "Family")

创建随机样本数据，并将其添加到组合数据集，确保样本名称sample_names和otu_table的行名匹配

#nsamples()得到样本的个数,sample_names()为获得样品名称
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))

现在使用ape包创建一个随机的系统发育树，并将其添加到数据集中,确保它的标签与OTU_table匹配

#ntaxa()获取分类单元/物种的数量，rtree()生成随机的树，taxa_names()获取分类单元/物种的名称
library("ape")
random_tree = rtree(ntaxa(physeq), rooted=TRUE, tip.label=taxa_names(physeq))
plot(random_tree)

使用merge_phyloseq将新的数据component添加到我们已经拥有的phyloseq对象中来.

#合并新数据与当前的phyloseq对象:
physeq1 = merge_phyloseq(physeq, sampledata, random_tree)
#使用我们刚刚生成的所有模拟数据组件从头重建phyloseq数据,与merge_phyloseq()函数功能一样
physeq2 = phyloseq(OTU, TAX, sampledata, random_tree)
#TRUE
identical(physeq1, physeq2)
#用新的组合数据构建几个树图plot_tree()函数,ladderize参数代表是否对树的结构进行梯形化
plot_tree(physeq1, color="Location", label.tips="taxa_names", ladderize="left", plot.margin=0.3)
#进行热图的绘制
plot_heatmap(physeq1)
#添加门类标签
plot_heatmap(physeq1, taxa.label="Phylum")

1.3 The import family of functions import_biom

QIIME的新版本生成了更全面、更正式定义的JSON或HDF5文件格式，被称为biom file format,biom文件格式(通常读作“biome”)是一种通用格式，用于表示一个或多个生物样本中的观察计数,BIOM是地球微生物组项目公认的标准，也是基因组学标准联盟候选项目
下面展示了如何导入四种主要类型的biom文件(实际上，您不需要知道您的文件是哪种类型，只需要知道它是biom文件)，import_biom函数允许您同时导入相关的系统发育树文件和参考序列文件(例如fasta文件)

#system.file()查找包中文件的完整文件名，定义文件路径
rich_dense_biom  = system.file("extdata", "rich_dense_otu_table.biom",  package="phyloseq")
rich_sparse_biom = system.file("extdata", "rich_sparse_otu_table.biom", package="phyloseq")
min_dense_biom   = system.file("extdata", "min_dense_otu_table.biom",   package="phyloseq")
min_sparse_biom  = system.file("extdata", "min_sparse_otu_table.biom",  package="phyloseq")
treefilename = system.file("extdata", "biom-tree.phy",  package="phyloseq")
refseqfilename = system.file("extdata", "biom-refseq.fasta",  package="phyloseq")
#树文件和引用序列文件都适用于任何示例biom文件，利用import_biom导入.parseFunction=parse_taxonomy_greengenes导入期间解决分类法相关错误时建议采取的第一步.
import_biom(rich_dense_biom, treefilename, refseqfilename, parseFunction=parse_taxonomy_greengenes)
#在实际操作中，您将把导入的结果存储为一些变量名，如myData，然后在下游数据操作和分析中使用这个数据对象
myData = import_biom(rich_dense_biom, treefilename, refseqfilename, parseFunction=parse_taxonomy_greengenes)
#查看myData中的对象，因为它是S4 object
mgData@otu_table
#利用myData中的数据绘制进化树
plot_tree(myData, color="Genus", shape="BODY_SITE", size="abundance")
#绘制各个身体位点的微生物丰度
plot_richness(myData, x="BODY_SITE", color="Description")
#包括绘制条形图
plot_bar(myData, fill="Genus")
#refseq()函数Retrieve reference sequences 的缩写，检索参考序列
refseq(myData)

1.4 import_qiime

可以使用phyloseq函数导入这些qiime软件的输入结果格式，特别是包含OTU-abundance和taxonomic identity信息的OTU文件

#其中gz文件也可以直接用read.table()函数读取
otufile = system.file("extdata", "GP_otu_table_rand_short.txt.gz", package="phyloseq")
mapfile = system.file("extdata", "master_map.txt", package="phyloseq")
#其中tree文件可以用ape包中的read.tree()函数读取
trefile = system.file("extdata", "GP_tree_rand_short.newick.gz", package="phyloseq")
rs_file = system.file("extdata", "qiime500-refseq.fasta", package="phyloseq")
#将这些数据整合读取利用import_qiime()函数
qiimedata = import_qiime(otufile, mapfile, trefile, rs_file)
#利用导入的文件进行分析，进行丰度的柱状图分析
plot_bar(qiimedata, x="SampleType", fill="Phylum")
#也包括热图的分析
plot_bar(qiimedata, x="SampleType", fill="Phylum")

其中也包括读取mothur软件的产出格式的文件，与import_qiime函数大同小异，利用函数import_mothur