Package 'phyloseq'

Description

There are already several ecology and phylogenetic packages available in R, including the adephylo, vegan, ade4, picante, ape, phangorn, phylobase, and OTUbase packages. These can already take advantage of many of the powerful statistical and graphics tools available in R. However, prior to phyloseq a user must devise their own methods for parsing the output of their favorite OTU clustering application, and, as a consequence, there is also no standard within Bioconductor (or R generally) for storing or sharing the suite of related data objects that describe a phylogenetic sequencing project. The phyloseq package seeks to address these issues by providing a related set of S4 classes that internally manage the handling tasks associated with organizing, linking, storing, and analyzing phylogenetic sequencing data. phyloseq additionally provides some convenience wrappers for input from common clustering applications, common analysis pipelines, and native implementation of methods that are not available in other R packages.

Author(s)

Paul J. McMurdie II [email protected]

References

www.stanford.edu/~mcmurdie

Description

See the documentation for the Extract generic, defined in the R base-package for the expected behavior.

Usage

## S4 method for signature 'otu_table,ANY,ANY,ANY'
x[i, j, ..., drop = TRUE]

## S4 method for signature 'sample_data,ANY,ANY,ANY'
x[i, j, ..., drop = TRUE]

## S4 method for signature 'taxonomyTable,ANY,ANY,ANY'
x[i, j, ..., drop = TRUE]

## S4 method for signature 'XStringSet,character,ANY,ANY'
x[i]

Arguments

Details

One special exception to standard behavior of these methods in phyloseq is that the drop argument is set internally to FALSE. This helps avoid bugs during complicated subsetting with multiple components, where it is necessary to be able to use a two dimensional indexing even if one of those dimensions has only 1 rank. Put another way, these phyloseq-defined extractions never collapse their result into a vector. See the documentation of Extract for more information about the drop argument.

See Also

Extract

Examples

data(esophagus)
nrow(otu_table(esophagus))
nrow(otu_table(esophagus)[1:5, ])

Description

This function is used internally by many accessors and in many functions/methods that need to access a particular type of component data. If something is wrong, or the slot is missing, the expected behavior is that this function will return NULL. Thus, the output can be tested by is.null as verification of the presence of a particular data component. Unlike the component-specific accessors (e.g. otu_table, or phy_tree), the default behavior is not to stop with an error if the desired slot is empty. In all cases this is controlled by the errorIfNULL argument, which can be set to TRUE if an error is desired.

Usage

access(physeq, slot, errorIfNULL=FALSE)

Arguments

Value

Returns the component object specified by the argument slot. Returns NULL if slot does not exist. Returns physeq as-is if it is a component class that already matches the slot name.

See Also

getslots.phyloseq, merge_phyloseq

Examples

#
## data(GlobalPatterns)
## access(GlobalPatterns, "tax_table")
## access(GlobalPatterns, "phy_tree")
## access(otu_table(GlobalPatterns), "otu_table")
## # Should return NULL:
## access(otu_table(GlobalPatterns), "sample_data")
## access(otuTree(GlobalPatterns), "sample_data")
## access(otuSam(GlobalPatterns), "phy_tree")

Build a tax_table from a named possibly-jagged list

Description

Build a tax_table from a named possibly-jagged list

Usage

build_tax_table(taxlist)

Arguments

Value

A tax_table (taxonomyTable-class) that has been built from taxlist. The OTU names of this output will be the element names of taxlist, and a separate taxonomic rank (column) will be included for each unique rank found among the element names of each vector in the list. NA_character_ is the default value of elements in the tax_table for which there is no corresponding information in taxlist.

See Also

import_biom import_qiime

Examples

taxvec1 = c("Root", "k__Bacteria", "p__Firmicutes", "c__Bacilli", "o__Bacillales", "f__Staphylococcaceae")
 parse_taxonomy_default(taxvec1)
 parse_taxonomy_greengenes(taxvec1)
 taxvec2 = c("Root;k__Bacteria;p__Firmicutes;c__Bacilli;o__Bacillales;f__Staphylococcaceae")
 parse_taxonomy_qiime(taxvec2)
 taxlist1 = list(OTU1=parse_taxonomy_greengenes(taxvec1), OTU2=parse_taxonomy_qiime(taxvec2))
 taxlist2 = list(OTU1=parse_taxonomy_default(taxvec1), OTU2=parse_taxonomy_qiime(taxvec2))
 build_tax_table(taxlist1)
 build_tax_table(taxlist2)

Description

Published in Nature in early 2011, this work compared (among other things), the faecal microbial communities from 22 subjects using complete shotgun DNA sequencing. Authors further compared these microbial communities with the faecal communities of subjects from other studies. A total of 280 faecal samples / subjects are represented in this dataset, and 553 genera. The authors claim that the data naturally clumps into three community-level clusters, or “enterotypes”, that are not immediately explained by sequencing technology or demographic features of the subjects, but with potential relevance to understanding human gut microbiota.

Details

abstract from research article (quoted):

Our knowledge of species and functional composition of the human gut microbiome is rapidly increasing, but it is still based on very few cohorts and little is known about variation across the world. By combining 22 newly sequenced faecal metagenomes of individuals from four countries with previously published data sets, here we identify three robust clusters (referred to as enterotypes hereafter) that are not nation or continent specific. We also confirmed the enterotypes in two published, larger cohorts, indicating that intestinal microbiota variation is generally stratified, not continuous. This indicates further the existence of a limited number of well-balanced host-microbial symbiotic states that might respond differently to diet and drug intake. The enterotypes are mostly driven by species composition, but abundant molecular functions are not necessarily provided by abundant species, highlighting the importance of a functional analysis to understand microbial communities. Although individual host properties such as body mass index, age, or gender cannot explain the observed enterotypes, data-driven marker genes or functional modules can be identified for each of these host properties. For example, twelve genes significantly correlate with age and three functional modules with the body mass index, hinting at a diagnostic potential of microbial markers.

(end quote)

Author(s)

Arumugam, M., Raes, J., et al.

References

Arumugam, M., et al. (2011). Enterotypes of the human gut microbiome.

Nature, 473(7346), 174-180.

http://www.nature.com/doifinder/10.1038/nature09944 See supplemental information for subject data.

OTU-clustered data was downloaded from the publicly-accessible:

http://www.bork.embl.de/Docu/Arumugam_et_al_2011/downloads.html

Examples

data(enterotype)
ig <- make_network(enterotype, "samples", max.dist=0.3)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.3, label=NULL)

Description

Includes just 3 samples, 1 each from 3 subjects. Although the research article mentions 4 subjects, only 3 are included in this dataset.

Details

abstract from research article (quoted):

The esophagus, like other luminal organs of the digestive system, provides a potential environment for bacterial colonization, but little is known about the presence of a bacterial biota or its nature. By using broad-range 16S rDNA PCR, biopsies were examined from the normal esophagus of four human adults. The 900 PCR products cloned represented 833 unique sequences belonging to 41 genera, or 95 species-level operational taxonomic units (SLOTU); 59 SLOTU were homologous with culture-defined bacterial species, 34 with 16S rDNA clones, and two were not homologous with any known bacterial 16S rDNA. Members of six phyla, Firmicutes, Bacteroides, Actinobacteria, Proteobacteria, Fusobacteria, and TM7, were represented. A large majority of clones belong to 13 of the 41 genera (783/900, 87%), or 14 SLOTU (574/900, 64%) that were shared by all four persons. Streptococcus (39%), Prevotella (17%), and Veilonella (14%) were most prevalent. The present study identified 56-79% of SLOTU in this bacterial ecosystem. Most SLOTU of esophageal biota are similar or identical to residents of the upstream oral biota, but the major distinction is that a large majority (82%) of the esophageal bacteria are known and cultivable. These findings provide evidence for a complex but conserved bacterial population in the normal distal esophagus.

(end quote)

A description of the 16S rRNA sequence processing can be found on the mothur-wiki at the link below. A cutoff of 0.10 was used for OTU clustering in that example, and it is taken here as well to create example data, esophagus, which was easily imported with the import_mothur() function.

Author(s)

Pei et al. [email protected]

References

Pei, Z., Bini, E. J., Yang, L., Zhou, M., Francois, F., & Blaser, M. J. (2004). Bacterial biota in the human distal esophagus. Proceedings of the National Academy of Sciences of the United States of America, 101(12), 4250-4255. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC384727

mothur-processed files and the sequence data can be downloaded from a zip-file, along with additional description, from the following URL: http://www.mothur.org/wiki/Esophageal_community_analysis

Examples

data(esophagus)
UniFrac(esophagus, weighted=TRUE)
# How to re-create the esophagus dataset using import_mothur function
mothlist  <- system.file("extdata", "esophagus.fn.list.gz", package="phyloseq")
mothgroup <- system.file("extdata", "esophagus.good.groups.gz", package="phyloseq")
mothtree  <- system.file("extdata", "esophagus.tree.gz", package="phyloseq")
show_mothur_cutoffs(mothlist)
cutoff    <- "0.10"
esophman  <- import_mothur(mothlist, mothgroup, mothtree, cutoff)

Description

Published in PNAS in early 2011. This work compared the microbial communities from 25 environmental samples and three known “mock communities” – a total of 9 sample types – at a depth averaging 3.1 million reads per sample. Authors were able to reproduce diversity patterns seen in many other published studies, while also invesitigating technical issues/bias by applying the same techniques to simulated microbial communities of known composition.

Details

abstract from research article (quoted):

The ongoing revolution in high-throughput sequencing continues to democratize the ability of small groups of investigators to map the microbial component of the biosphere. In particular, the coevolution of new sequencing platforms and new software tools allows data acquisition and analysis on an unprecedented scale. Here we report the next stage in this coevolutionary arms race, using the Illumina GAIIx platform to sequence a diverse array of 25 environmental samples and three known “mock communities” at a depth averaging 3.1 million reads per sample. We demonstrate excellent consistency in taxonomic recovery and recapture diversity patterns that were previously reported on the basis of metaanalysis of many studies from the literature (notably, the saline/nonsaline split in environmental samples and the split between host-associated and free-living communities). We also demonstrate that 2,000 Illumina single-end reads are sufficient to recapture the same relationships among samples that we observe with the full dataset. The results thus open up the possibility of conducting large-scale studies analyzing thousands of samples simultaneously to survey microbial communities at an unprecedented spatial and temporal resolution.

(end quote)

Many thanks to J. Gregory Caporaso for directly providing the OTU-clustered data files for inclusion in this package.

Author(s)

Caporaso, J. G., et al.

References

Caporaso, J. G., et al. (2011). Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. PNAS, 108, 4516-4522. PMCID: PMC3063599

The primary article can be viewed/downloaded at: http://www.pnas.org/content/108/suppl.1/4516.short

See Also

The examples on the phyloseq wiki page for plot_ordination show many more examples:

https://github.com/joey711/phyloseq/wiki/plot_ordination

Examples

data(GlobalPatterns)
plot_richness(GlobalPatterns, x="SampleType", measures=c("Observed", "Chao1", "Shannon"))

Description

Published in early 2011, this work compared 24 separate soil microbial communities under four treatment conditions via multiplexed/barcoded 454-pyrosequencing of PCR-amplified 16S rRNA gene fragments. The authors found differences in the composition and structure of microbial communities between soil treatments. As expected, the soil microbial communities were highly diverse, with a staggering 16,825 different OTUs (species) observed in the included dataset. Interestingly, this study used a larger number of replicates than previous studies of this type, for a total of 56 samples, and the putatively low resampling rate of species between replicated sequencing trials (“OTU overlap”) was a major concern by the authors.

Details

This dataset contains an experiment-level (phyloseq-class) object, which in turn contains the taxa-contingency table and soil-treatment table as otu_table-class and sample_data-class components, respectively.

This data was imported from raw files supplied directly by the authors via personal communication for the purposes of including as an example in the phyloseq-package. As this data is sensitive to choices in OTU-clustering parameters, attempts to recreate the otu_table from the raw sequencing data may give slightly different results than the table provided here.

abstract from research article (quoted):

To determine the reproducibility and quantitation of the amplicon sequencing-based detection approach for analyzing microbial community structure, a total of 24 microbial communities from a long-term global change experimental site were examined. Genomic DNA obtained from each community was used to amplify 16S rRNA genes with two or three barcode tags as technical replicates in the presence of a small quantity (0.1% wt/wt) of genomic DNA from Shewanella oneidensis MR-1 as the control. The technical reproducibility of the amplicon sequencing-based detection approach is quite low, with an average operational taxonomic unit (OTU) overlap of 17.2%+/-2.3% between two technical replicates, and 8.2%+/-2.3% among three technical replicates, which is most likely due to problems associated with random sampling processes. Such variations in technical replicates could have substantial effects on estimating beta-diversity but less on alpha-diversity. A high variation was also observed in the control across different samples (for example, 66.7-fold for the forward primer), suggesting that the amplicon sequencing-based detection approach could not be quantitative. In addition, various strategies were examined to improve the comparability of amplicon sequencing data, such as increasing biological replicates, and removing singleton sequences and less-representative OTUs across biological replicates. Finally, as expected, various statistical analyses with preprocessed experimental data revealed clear differences in the composition and structure of microbial communities between warming and non-warming, or between clipping and non-clipping. Taken together, these results suggest that amplicon sequencing-based detection is useful in analyzing microbial community structure even though it is not reproducible and quantitative. However, great caution should be taken in experimental design and data interpretation when the amplicon sequencing-based detection approach is used for quantitative analysis of the beta-diversity of microbial communities.

(end quote)

Author(s)

Jizhong Zhou, et al.

References

Zhou, J., Wu, L., Deng, Y., Zhi, X., Jiang, Y.-H., Tu, Q., Xie, J., et al. Reproducibility and quantitation of amplicon sequencing-based detection. The ISME Journal. (2011) 5(8):1303-1313. doi:10.1038/ismej.2011.11

The article can be accessed online at http://www.nature.com/ismej/journal/v5/n8/full/ismej201111a.html

Examples

# Load the data
data(soilrep)
################################################################################
# Alpha diversity (richness) example. Accept null hypothesis: 
# No convincing difference in species richness between warmed/unwarmed soils.
################################################################################
# Graphically compare richness between the different treatments.
man.col <- c(WC="red", WU="brown", UC="blue", UU="darkgreen")
plot_richness(soilrep, x="Treatment", color="Treatment", measures=c("Observed", "Chao1", "Shannon"))

An S4 placeholder for the dist class.

Description

See dist for details about this type of a distance matrix object.

See Also

dist, setOldClass

Description

Takes a phyloseq-class object and method option, and returns a distance object suitable for certain ordination methods and other distance-based analyses. Only sample-wise distances are currently supported (the type argument), but eventually species-wise (OTU-wise) distances may be supported as well.

Usage

distance(physeq, method, type = "samples", ...)

## S4 method for signature 'phyloseq,ANY'
distance(physeq, method)

## S4 method for signature 'otu_table,character'
distance(physeq, method,
  type = "samples", ...)

## S4 method for signature 'phyloseq,character'
distance(physeq, method, type = "samples",
  ...)

Arguments

Details

Depending on the method argument, distance() wraps one of UniFrac, DPCoA, JSD, vegdist, betadiver, designdist, or dist.

Value

An object of class “dist” suitable for certain ordination methods and other distance-based analyses.

See Also

plot_ordination, UniFrac, DPCoA, JSD, vegdist, betadiver, designdist, dist.

Examples

data(esophagus)
distance(esophagus, "uunifrac") # Unweighted UniFrac
distance(esophagus, "wunifrac") # weighted UniFrac
distance(esophagus, "jaccard", binary = TRUE) # vegdist jaccard
distance(esophagus, "gower") # vegdist option "gower"
distance(esophagus, "g") # designdist method option "g"
distance(esophagus, "minkowski") # invokes a method from the base dist() function.
distance(esophagus, "(A+B-2*J)/(A+B)") # designdist custom distance
distanceMethodList
help("distance")

List of distance method keys supported in distance

Description

Distance methods should be specified by exact string match. Cannot do partial matching for all options, because too many similar options in downstream method dispatch.

Usage

distanceMethodList

Format

A list of character vectors. Every entry specifies a supported distance method. Names in the list indicate which downstream function is being utilized for further details. Same functions are linked in the itemized list below.

unifrac

UniFrac

wunifrac

UniFrac

dpcoa

DPCoA

jsd

JSD

manhattan

vegdist

euclidean

vegdist

canberra

vegdist

bray

vegdist

kulczynski

vegdist

jaccard

vegdist

gower

vegdist

altGower

vegdist

morisita

vegdist

horn

vegdist

mountford

vegdist

raup

vegdist

binomial

vegdist

chao

vegdist

cao

vegdist

w

betadiver

-

betadiver

c

betadiver

wb

betadiver

r

betadiver

I

betadiver

e

betadiver

t

betadiver

me

betadiver

j

betadiver

sor

betadiver

m

betadiver

-

betadiver

co

betadiver

cc

betadiver

g

betadiver

-

betadiver

l

betadiver

hk

betadiver

rlb

betadiver

sim

betadiver

gl

betadiver

z

betadiver

maximum

dist

binary

dist

minkowski

dist

ANY

designdist

See Also

distance

Examples

distanceMethodList

Description

Function uses abundance (otu_table-class) and phylogenetic (phylo) components of a phyloseq-class experiment-level object to perform a Double Principle Coordinate Analysis (DPCoA), relying heavily on the underlying (and more general) function, dpcoa. The distance object ultimately provided is the square root of the cophenetic/patristic (cophenetic.phylo) distance between the species, which is always Euclidean.

Although this distance is Euclidean, for numerical reasons it will sometimes look non-Euclidean, and a correction will be performed. See correction argument.

Usage

DPCoA(physeq, correction = cailliez, scannf = FALSE, ...)

Arguments

Value

A dpcoa-class object (see dpcoa).

Author(s)

Julia Fukuyama [email protected]. Adapted for phyloseq by Paul J. McMurdie.

References

Pavoine, S., Dufour, A.B. and Chessel, D. (2004) From dissimilarities among species to dissimilarities among communities: a double principal coordinate analysis. Journal of Theoretical Biology, 228, 523-537.

See Also

dpcoa

Examples

# # # # # # Esophagus
data(esophagus)
eso.dpcoa <- DPCoA(esophagus)
eso.dpcoa
plot_ordination(esophagus, eso.dpcoa, "samples")
plot_ordination(esophagus, eso.dpcoa, "species")
plot_ordination(esophagus, eso.dpcoa, "biplot")
#
#
# # # # # # GlobalPatterns
data(GlobalPatterns)
# subset GP to top-150 taxa (to save computation time in example)
keepTaxa <- names(sort(taxa_sums(GlobalPatterns), TRUE)[1:150])
GP       <- prune_taxa(keepTaxa, GlobalPatterns)
# Perform DPCoA
GP.dpcoa <- DPCoA(GP)
plot_ordination(GP, GP.dpcoa, color="SampleType")

Description

Performs a number of standard alpha diversity estimates, and returns the results as a data.frame. Strictly speaking, this function is not only estimating richness, despite its name. It can operate on the cumulative population of all samples in the dataset, or by repeating the richness estimates for each sample individually. NOTE: You must use untrimmed datasets for meaningful results, as these estimates (and even the “observed” richness) are highly dependent on the number of singletons. You can always trim the data later on if needed, just not before using this function.

Usage

estimate_richness(physeq, split = TRUE, measures = NULL)

Arguments

Value

A data.frame of the richness estimates, and their standard error.

See Also

Check out the custom plotting function, plot_richness, for easily showing the results of different estimates, with method-specific error-bars. Also check out the internal functions borrowed from the vegan package:

estimateR

diversity

fisherfit

Examples

## There are many more interesting examples at the phyloseq online tutorials.
## http://joey711.github.com/phyloseq/plot_richness-examples
 data("esophagus")
 # Default is all available measures
 estimate_richness(esophagus)
 # Specify just one:
 estimate_richness(esophagus, measures="Observed")
 # Specify a few:
 estimate_richness(esophagus, measures=c("Observed", "InvSimpson", "Shannon", "Chao1"))

Description

Creates the environment table that is needed for the original UniFrac algorithm. Useful for cross-checking, or if want to use UniFrac server. Optionally the ENV-formatted table can be returned to the R workspace, and the tree component can be exported as Nexus format (Recommended).

Usage

export_env_file(physeq, file = "", writeTree = TRUE, return = FALSE)

Arguments

Examples

# # Load example data
# data(esophagus)
# export_env_file(esophagus, "~/Desktop/esophagus.txt")

Description

The purpose of this function is to allow a user to easily export a distance object as a pair of files that can be immediately imported by mothur for OTU clustering and related analysis. A distance object can be created in R in a number of ways, including via cataloguing the cophentic distances of a tree object.

Usage

export_mothur_dist(x, out=NULL, makeTrivialNamesFile=NULL)

Arguments

Value

A character vector of the different cutoff values contained in the file. For a given set of arguments to the cluster() command from within mothur, a number of OTU-clustering results are returned in the same list file. The exact cutoff values used by mothur can vary depending on the input data. This simple function returns the cutoffs that were actually included in the mothur output. This an important extra step prior to importing the OTUs with the import_mothur_otulist() function.

Examples

#
data(esophagus) 
myDistObject <- as.dist(ape::cophenetic.phylo(phy_tree(esophagus)))
export_mothur_dist(myDistObject)

Description

This function is directly analogous to the genefilter function for microarray filtering, but is used for filtering OTUs from phyloseq objects. It applies an arbitrary set of functions — as a function list, for instance, created by filterfun — as across-sample criteria, one OTU at a time. It takes as input a phyloseq object, and returns a logical vector indicating whether or not each OTU passed the criteria. Alternatively, if the "prune" option is set to FALSE, it returns the already-trimmed version of the phyloseq object.

Usage

filter_taxa(physeq, flist, prune=FALSE)

Arguments

Value

A logical vector equal to the number of taxa in physeq. This can be provided directly to prune_taxa as first argument. Alternatively, if prune==TRUE, the pruned phyloseq-class object is returned instead.

See Also

filterfun, genefilter_sample, filterfun_sample

Examples

data("enterotype")
 require("genefilter")
 flist    <- filterfun(kOverA(5, 2e-05))
 ent.logi <- filter_taxa(enterotype, flist)
 ent.trim <- filter_taxa(enterotype, flist, TRUE)
 identical(ent.trim, prune_taxa(ent.logi, enterotype)) 
 identical(sum(ent.logi), ntaxa(ent.trim))
 filter_taxa(enterotype, flist, TRUE)

A sample-wise filter function builder analogous to filterfun.

Description

See the filterfun, from the Bioconductor repository, for a taxa-/gene-wise filter (and further examples).

Usage

filterfun_sample(...)

Arguments

Value

An enclosure (function) that itself will return a logical vector, according to the functions provided in the argument list, evaluated in order. The output of filterfun_sample is appropriate for the ‘flist’ argument to the genefilter_sample method.

See Also

filterfun, genefilter_sample

Examples

# Use simulated abundance matrix
set.seed(711)
testOTU <- otu_table(matrix(sample(1:50, 25, replace=TRUE), 5, 5), taxa_are_rows=FALSE)
f1  <- filterfun_sample(topk(2))
wh1 <- genefilter_sample(testOTU, f1, A=2)
wh2 <- c(TRUE, TRUE, TRUE, FALSE, FALSE)
prune_taxa(wh1, testOTU)
prune_taxa(wh2, testOTU)

Description

This is a wrapper for the clusGap function, expecting an ordination result as the main data argument.

Usage

gapstat_ord(ord, axes = c(1:2), type = "sites",
  FUNcluster = function(x, k) {     list(cluster = pam(x, k, cluster.only
  = TRUE)) }, K.max = 8, ...)

Arguments

Value

An object of S3 class "clusGap", basically a list with components. See the clusGap documentation for more details.

Examples

data("soilrep")
sord  = ordinate(soilrep, "PCoA", "bray")
# Evaluate axes with scree plot
plot_scree(sord)
# Gap Statistic
gs = gapstat_ord(sord, axes=1:3, verbose=FALSE)
# plot_ordination(soilrep, sord,  color="Treatment")
plot_clusgap(gs)
print(gs, method="Tibs2001SEmax")

Description

A general OTU trimming function for selecting OTUs that satisfy some criteria within the distribution of each sample, and then also an additional criteria for number of samples that must pass. This is a genefilter-like function that only considers sample-wise criteria. The number of acceptable samples is used as the final criteria (set by the argument A) to determine whether or not the taxa should be retained (TRUE) or not (FALSE). Just like with genefilter, a logical having length equal to nrow()/ntaxa is returned, indicating which should be kept. This output can be provided directly to OTU trimming function, prune_taxa. By contrast, genefilter, of the genefilter package in Bioconductor, works only on the rows of a matrix. Note that, because otu_table-class inherits directly from the matrix-class, an unmodified otu_table can be provided to genefilter, but be mindful of the orientation of the otu_table (use taxa_are_rows), and transpose (t) if needed.

Usage

genefilter_sample(X, flist, A=1)

## S4 method for signature 'matrix'
genefilter_sample(X, flist, A = 1)

## S4 method for signature 'otu_table'
genefilter_sample(X, flist, A = 1)

## S4 method for signature 'phyloseq'
genefilter_sample(X, flist, A = 1)

Arguments

Value

A logical vector with names equal to taxa_names (or rownames, if matrix).

See Also

genefilter, filterfun_sample, t, prune_taxa

Examples

#
## testOTU <- otu_table(matrix(sample(1:50, 25, replace=TRUE), 5, 5), taxa_are_rows=FALSE)
## f1  <- filterfun_sample(topk(2))
## wh1 <- genefilter_sample(testOTU, f1, A=2)
## wh2 <- c(TRUE, TRUE, TRUE, FALSE, FALSE)
## prune_taxa(wh1, testOTU)
## prune_taxa(wh2, testOTU)
## 
## tax_table1 <- tax_table(matrix("abc", 5, 5))
## prune_taxa(wh1, tax_table1)
## prune_taxa(wh2, tax_table1)

Description

This is a simple accessor function for investigating a single species-of-interest.

Usage

get_sample(physeq, i)

## S4 method for signature 'otu_table'
get_sample(physeq, i)

## S4 method for signature 'phyloseq'
get_sample(physeq, i)

Arguments

Value

An integer vector of the abundance values for each sample in physeq for species i

See Also

get_taxa taxa_names sample_names

Examples

data(esophagus)
taxa_names(esophagus)
get_sample(esophagus, "59_5_19")

Description

This is a simple accessor function for investigating a single sample-of-interest.

Usage

get_taxa(physeq, i)

## S4 method for signature 'otu_table'
get_taxa(physeq, i)

## S4 method for signature 'phyloseq'
get_taxa(physeq, i)

Arguments

Value

An integer vector of the abundance values for each species in physeq for sample i

See Also

get_sample taxa_names sample_names

Examples

data(esophagus)
sample_names(esophagus)
get_taxa(esophagus, "B")

Description

This is a simple accessor function to make it more convenient to determine the different taxa present for a particular taxonomic rank in a given phyloseq-class object.

Usage

get_taxa_unique(physeq, taxonomic.rank=rank_names(physeq)[1], errorIfNULL=TRUE)

Arguments

Value

Character vector. Unique vector of the observed taxa at a particular taxonomic rank

See Also

get_taxa taxa_names sample_names

Examples

data(enterotype)
get_taxa_unique(enterotype)
data(GlobalPatterns)
get_taxa_unique(GlobalPatterns, "Family")

Description

This is a simple accessor function for streamlining access to values/vectors/factors/etc contained in the sample_data.

Usage

get_variable(physeq, varName)

Arguments

Value

Data. The clas of the data depends on what the contents of sample_data.

See Also

get_taxa taxa_names sample_names

sample_variables

Examples

# Load the GlobalPatterns dataset into the workspace environment
data(GlobalPatterns)
# Look at the different values for SampleType 
get_variable(GlobalPatterns, "SampleType")

Description

Like getSlots, but returns the class name if argument is component data object.

Usage

getslots.phyloseq(physeq)

Arguments

Value

identical to getSlots. A named character vector of the slot classes of a particular S4 class, where each element is named by the slot name it represents. If physeq is a component data object, then a vector of length (1) is returned, named according to its slot name in the phyloseq-class.

See Also

merge_phyloseq

Examples

#
 data(GlobalPatterns)
 getslots.phyloseq(GlobalPatterns)
 data(esophagus)
 getslots.phyloseq(esophagus)

Description

A user must still understand the additional arguments required for each type of import data. Those arguments are described in detail at the tool-specific import_* links below. Each clustering tool / package / pipeline has its own idiosyncratic set of file names / types, and it remains the responsibility of the user to understand which file-path should be provided to each argument for the particular importing submethod. This method merely provides a central documentation and method-name, and the arguments are passed along as-is.

Usage

import(pipelineName, ...)

Arguments

Value

In most cases a phyloseq-class will be returned, though the included component data will vary by pipeline/tool, and also by the types of data files provided. The expected behavior is to return the most-comprehensive object possible, given the provided arguments and pipeline/tool.

References

BIOM: http://www.biom-format.org/

mothur: http://www.mothur.org/wiki/Main_Page

PyroTagger: http://pyrotagger.jgi-psf.org/

QIIME: http://qiime.org/

RDP pipeline: http://pyro.cme.msu.edu/index.jsp

See Also

For BIOM format, see: import_biom

For mothur, see: import_mothur

Separate tools for mothur are also: show_mothur_cutoffs import_mothur_dist export_mothur_dist

For PyroTagger, see: import_pyrotagger_tab

For QIIME legacy format, see: import_qiime

For RDP pipeline, see: import_RDP_cluster

import_RDP_otu

Examples

## See documentation of a specific import function

Description

New versions of QIIME produce a more-comprehensive and formally-defined JSON file format, called biom file format:

Usage

import_biom(BIOMfilename, 
 treefilename=NULL, refseqfilename=NULL, refseqFunction=readDNAStringSet, refseqArgs=NULL,
 parseFunction=parse_taxonomy_default, parallel=FALSE, version=1.0, ...)

Arguments

Details

“The biom file format (canonically pronounced ‘biome’) is designed to be a general-use format for representing counts of observations in one or more biological samples. BIOM is a recognized standard for the Earth Microbiome Project and is a Genomics Standards Consortium candidate project.”

http://biom-format.org/

Value

A phyloseq-class object.

References

biom-format

See Also

import

import_qiime

read_tree

read_tree_greengenes

read_biom

biom_data

sample_metadata

observation_metadata

XStringSet-io

Examples

# An included example of a rich dense biom file
rich_dense_biom  <- system.file("extdata", "rich_dense_otu_table.biom",  package="phyloseq")
import_biom(rich_dense_biom,  parseFunction=parse_taxonomy_greengenes)
# An included example of a sparse dense biom file
rich_sparse_biom <- system.file("extdata", "rich_sparse_otu_table.biom", package="phyloseq")
import_biom(rich_sparse_biom, parseFunction=parse_taxonomy_greengenes)
# # # Example code for importing large file with parallel backend
# library("doParallel")
# registerDoParallel(cores=6)
# import_biom("my/file/path/file.biom", parseFunction=parse_taxonomy_greengenes, parallel=TRUE)

Description

Convenience wrapper function to read the environment-file, as formatted for input to the UniFrac server (http://bmf2.colorado.edu/unifrac/). The official format of these files is that each row specifies (in order) the sequence name, source sample, and (optionally) the number of times the sequence was observed.

Usage

import_env_file(envfilename, tree=NULL, sep="\t", ...)

Arguments

Value

An otu_table-class, or phyloseq-class if a phylo-class argument is provided to tree.

References

http://bmf2.colorado.edu/unifrac/

See Also

import

tip_glom

Examples

# import_env_file(myEnvFile, myTree)

Description

Technically all parameters are optional, but if you don't provide any file connections, then nothing will be returned. While the list and group files are the first two arguments for legacy-compatibility reasons, we don't recommend that you use these file types with modern (large) datasets. They are comically inefficient, as they store the name of every sequencing read in both files. The mothur package provides conversions utilities to create other more-efficient formats, which we recommend, like the shared file for an OTU table. Alternatively, mothur also provides a utility to create a biom-format file that is independent of OTU clustering platform. Biom-format files should be imported not with this function, but with import_biom. The resulting objects after import should be identical in R.

Usage

import_mothur(mothur_list_file = NULL, mothur_group_file = NULL,
  mothur_tree_file = NULL, cutoff = NULL, mothur_shared_file = NULL,
  mothur_constaxonomy_file = NULL,
  parseFunction = parse_taxonomy_default)

Arguments

Value

The object class depends on the provided arguments. A phyloseq object is returned if enough data types are provided. If only one data component can be created from the data, it is returned.

FASTER (recommended for larger data sizes):

If only a mothur_constaxonomy_file is provided, then a taxonomyTable-class object is returned.

If only a mothur_shared_file is provided, then an otu_table object is returned.

SLOWER (but fine for small file sizes):

The list and group file formats are extremely inefficient for large datasets, and they are not recommended. The mothur software provides tools for converting to other file formats, such as a so-called “shared” file. You should provide a shared file, or group/list files, but not both at the same time. If only a list and group file are provided, then an otu_table object is returned. Similarly, if only a list and tree file are provided, then only a tree is returned (phylo-class).

References

http://www.mothur.org/wiki/Main_Page

Schloss, P.D., et al., Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol, 2009. 75(23):7537-41.

Examples

# # The following example assumes you have downloaded the esophagus example
# # dataset from the mothur wiki:
# # "http://www.mothur.org/wiki/Esophageal_community_analysis"
# # "http://www.mothur.org/w/images/5/55/Esophagus.zip"
# # The path on your machine may (probably will) vary
# mothur_list_file  <- "~/Downloads/mothur/Esophagus/esophagus.an.list"
# mothur_group_file <- "~/Downloads/mothur/Esophagus/esophagus.good.groups"
# mothur_tree_file  <- "~/Downloads/mothur/Esophagus/esophagus.tree"
# # # Actual examples follow:
# show_mothur_cutoffs(mothur_list_file)
# test1 <- import_mothur(mothur_list_file, mothur_group_file, mothur_tree_file)
# test2 <- import_mothur(mothur_list_file, mothur_group_file, mothur_tree_file, cutoff="0.02")
# # Returns just a tree
# import_mothur(mothur_list_file, mothur_tree_file=mothur_tree_file)
# # Returns just an otu_table
# import_mothur(mothur_list_file, mothur_group_file=mothur_group_file)
# # Returns an error
# import_mothur(mothur_list_file)
# # Should return an "OMG, you must provide the list file" error
# import_mothur()

Description

The mothur application will produce a file containing the pairwise distances between all sequences in a dataset. This distance matrix can be the basis for OTU cluster designations. R also has many built-in or off-the-shelf tools for dealing with distance matrices.

Usage

import_mothur_dist(mothur_dist_file)

Arguments

Value

A distance matrix object describing all sequences in a dataset.

See Also

import_mothur

Examples

# # Take a look at the dataset shown here as an example:
# # "http://www.mothur.org/wiki/Esophageal_community_analysis"
# # find the file ending with extension ".dist", download to your system
# # The location of your file may vary
# mothur_dist_file <- "~/Downloads/mothur/Esophagus/esophagus.dist"
# myNewDistObject  <- import_mothur_dist(mothur_dist_file)

Description

PyroTagger is a web-server that takes raw, barcoded 16S rRNA amplicon sequences and returns an excel spreadsheet (".xls") with both abundance and taxonomy data. It also includes some confidence information related to the taxonomic assignment.

Usage

import_pyrotagger_tab(pyrotagger_tab_file, 
strict_taxonomy=FALSE, keep_potential_chimeras=FALSE)

Arguments

Details

PyroTagger is created and maintained by the Joint Genome Institute at "http://pyrotagger.jgi-psf.org/"

The typical output form PyroTagger is a spreadsheet format ".xls", which poses additional import challenges. However, virtually all spreadsheet applications support the ".xls" format, and can further export this file in a tab-delimited format. It is recommended that you convert the xls-file without any modification (as tempting as it might be once you have loaded it) into a tab-delimited text file. Deselect any options to encapsulate fields in quotes, as extra quotes around each cell's contents might cause problems during file processing. These quotes will also inflate the file-size, so leave them out as much as possible, while also resisting any temptation to modify the xls-file “by hand”.

A highly-functional and free spreadsheet application can be obtained as part of the cross-platform OpenOffice suite. It works for the above required conversion. Go to "http://www.openoffice.org/".

It is regrettable that this importer does not take the xls-file directly as input. However, because of the moving-target nature of spreadsheet file formats, there is limited support for direct import of these formats into R. Rather than add to the dependency requirements of emphphyloseq and the relative support of these xls-support packages, it seems more efficient to choose an arbitrary delimited text format, and focus on the data structure in the PyroTagger output. This will be easier to support in the long-run.

Value

An otuTax object containing both the otu_table and TaxonomyTable data components, parsed from the pyrotagger output.

References

http://pyrotagger.jgi-psf.org/

Examples

## New_otuTaxObject <- import_pyrotagger_tab(pyrotagger_tab_file)

Description

QIIME produces several files that can be directly imported by the phyloseq-package. Originally, QIIME produced its own custom format table that contained both OTU-abundance and taxonomic identity information. This function is still included in phyloseq mainly to accommodate these now-outdated files. Recent versions of QIIME store output in the biom-format, an emerging file format standard for microbiome data. If your data is in the biom-format, if it ends with a .biom file name extension, then you should use the import_biom function instead.

Usage

import_qiime(otufilename = NULL, mapfilename = NULL,
  treefilename = NULL, refseqfilename = NULL,
  refseqFunction = readDNAStringSet, refseqArgs = NULL,
  parseFunction = parse_taxonomy_qiime, verbose = TRUE, ...)

Arguments

Details

Other related files include the mapping-file that typically stores sample covariates, converted naturally to the sample_data-class component data type in the phyloseq-package. QIIME may also produce a phylogenetic tree with a tip for each OTU, which can also be imported specified here or imported separately using read_tree.

See "http://www.qiime.org/" for details on using QIIME. While there are many complex dependencies, QIIME can be downloaded as a pre-installed linux virtual machine that runs “off the shelf”.

The different files useful for import to phyloseq are not collocated in a typical run of the QIIME pipeline. See the main phyloseq vignette for an example of where ot find the relevant files in the output directory.

Value

A phyloseq-class object.

References

http://qiime.org/

“QIIME allows analysis of high-throughput community sequencing data.” J Gregory Caporaso, Justin Kuczynski, Jesse Stombaugh, Kyle Bittinger, Frederic D Bushman, Elizabeth K Costello, Noah Fierer, Antonio Gonzalez Pena, Julia K Goodrich, Jeffrey I Gordon, Gavin A Huttley, Scott T Kelley, Dan Knights, Jeremy E Koenig, Ruth E Ley, Catherine A Lozupone, Daniel McDonald, Brian D Muegge, Meg Pirrung, Jens Reeder, Joel R Sevinsky, Peter J Turnbaugh, William A Walters, Jeremy Widmann, Tanya Yatsunenko, Jesse Zaneveld and Rob Knight; Nature Methods, 2010; doi:10.1038/nmeth.f.303

See Also

phyloseq

merge_phyloseq

read_tree

read_tree_greengenes

XStringSet-io

Examples

otufile <- system.file("extdata", "GP_otu_table_rand_short.txt.gz", package="phyloseq")
 mapfile <- system.file("extdata", "master_map.txt", package="phyloseq")
 trefile <- system.file("extdata", "GP_tree_rand_short.newick.gz", package="phyloseq")
 import_qiime(otufile, mapfile, trefile)

Description

Now a legacy-format, older versions of QIIME produced an OTU file that typically contains both OTU-abundance and taxonomic identity information in a tab-delimted table. If your file ends with the extension .biom, or if you happen to know that it is a biom-format file, or if you used default settings in a version of QIIME of 1.7 or greater, then YOU SHOULD USE THE BIOM-IMPORT FUNCTION instead, import_biom.

Usage

import_qiime_otu_tax(file, parseFunction = parse_taxonomy_qiime,
  verbose = TRUE, parallel = FALSE)

Arguments

Details

This function uses chunking to perform both the reading and parsing in blocks of optional size, thus constrain the peak memory usage. feature should make this importer accessible to machines with modest memory, but with the caveat that the full numeric matrix must be a manageable size at the end, too. In principle, the final tables will be large, but much more efficiently represented than the character-stored numbers. If total memory for storing the numeric matrix becomes problematic, a switch to a sparse matrix representation of the abundance – which is typically well-suited to this data – might provide a solution.

Value

A list of two matrices. $otutab contains the OTU Table as a numeric matrix, while $taxtab contains a character matrix of the taxonomy assignments.

See Also

import

merge_phyloseq

phyloseq

import_qiime

read_tree

read_tree_greengenes

import_env_file

Examples

otufile <- system.file("extdata", "GP_otu_table_rand_short.txt.gz", package="phyloseq")
 import_qiime_otu_tax(otufile)

Description

QIIME produces several files that can be analyzed in the phyloseq-package, This includes the map-file, which is an important input to QIIME that can also indicate sample covariates. It is converted naturally to the sample_data component data type in phyloseq-package, based on the R data.frame.

Usage

import_qiime_sample_data(mapfilename)

Arguments

Details

See import_qiime for more information about QIIME. It is also the suggested function for importing QIIME-produced data files.

Value

A sample_data object.

See Also

import

merge_phyloseq

phyloseq

import_qiime

import_qiime_otu_tax

import_env_file

Examples

mapfile <- system.file("extdata", "master_map.txt", package = "phyloseq")
 import_qiime_sample_data(mapfile)

Description

The RDP cluster pipeline (specifically, the output of the complete linkage clustering step) has no formal documentation for the ".clust" file or its apparent sequence naming convention.

Usage

import_RDP_cluster(RDP_cluster_file)

Arguments

Details

http://pyro.cme.msu.edu/index.jsp

The cluster file itself contains the names of all sequences contained in input alignment. If the upstream barcode and aligment processing steps are also done with the RDP pipeline, then the sequence names follow a predictable naming convention wherein each sequence is named by its sample and sequence ID, separated by a "_" as delimiter:

"sampleName_sequenceIDnumber"

This import function assumes that the sequence names in the cluster file follow this convention, and that the sample name does not contain any "_". It is unlikely to work if this is not the case. It is likely to work if you used the upstream steps in the RDP pipeline to process your raw (barcoded, untrimmed) fasta/fastq data.

This function first loops through the ".clust" file and collects all of the sample names that appear. It secondly loops through each OTU ("cluster"; each row of the cluster file) and sums the number of sequences (reads) from each sample. The resulting abundance table of OTU-by-sample is trivially coerced to an otu_table object, and returned.

Value

An otu_table object parsed from the ".clust" file.

References

http://pyro.cme.msu.edu/index.jsp

Description

Recently updated tools on RDP Pyro site make it easier to import Pyrosequencing output into R. The modified tool “Cluster To R Formatter” can take a cluster file (generated from RDP Clustering tools) to create a community data matrix file for distance cutoff range you are interested in. The resulting output file is a tab-delimited file containing the number of sequences for each sample for each OTU. The OTU header naming convention is "OTU_" followed by the OTU number in the cluster file. It pads “0”s to make the OTU header easy to sort. The OTU numbers are not necessarily in order.

Usage

import_RDP_otu(otufile)

Arguments

Value

A otu_table-class object.

See Also

An alternative “cluster” file importer for RDP results: import_RDP_cluster

The main RDP-pyrosequencing website http://pyro.cme.msu.edu/index.jsp

Examples

otufile <- system.file("extdata", "rformat_dist_0.03.txt.gz", package="phyloseq")
### the gzipped file is automatically recognized, and read using R-connections
ex_otu  <- import_RDP_otu(otufile)
class(ex_otu)
ntaxa(ex_otu)
nsamples(ex_otu)
sample_sums(ex_otu)
head(t(ex_otu))

Import UPARSE file format

Description

UPARSE is an algorithm for OTU-clustering implemented within usearch. At last check, the UPARSE algortihm was accessed via the -cluster_otu option flag. For details about installing and running usearch, please refer to the usearch website. For details about the output format, please refer to the uparse format definition.

Usage

import_uparse(upFile, omitChimeras = TRUE, countTable = TRUE,
  OTUtable = TRUE, verbose = TRUE)

Arguments

Details

Because UPARSE is an external (non-R) application, there is no direct way to continuously check that these suggested arguments and file formats will remain in their current state. If there is a problem, please verify your version of usearch, create a small reproducible example of the problem, and post it as an issue on the phyloseq issues tracker.

See Also

import_usearch_uc

Examples

###

Description

UPARSE is an algorithm for OTU-clustering implemented within usearch. At last check, the UPARSE algortihm was accessed via the -cluster_otu option flag. For details about installing and running usearch, please refer to the usearch website. For details about the output format, please refer to the uc format definition. This importer is intended to read a particular table format output that is generated by usearch, its so-called “cluster format”, a file format that is often given the .uc extension in usearch documentation.

Usage

import_usearch_uc(ucfile, colRead = 9, colOTU = 10,
  readDelimiter = "_", verbose = TRUE)

Arguments

Details

Because usearch is an external (non-R) application, there is no direct way to continuously check that these suggested arguments and file formats will remain in their current state. If there is a problem, please verify your version of usearch, create a small reproducible example of the problem, and post it as an issue on the phyloseq issues tracker. The version of usearch upon which this import function was created is 7.0.109. Hopefully later versions of usearch maintain this function and format, but the phyloseq team has no way to guarantee this, and so any feedback about this will help maintain future functionality. For instance, it is currently assumed that the 9th and 10th columns of the .uc table hold the read-label and OTU ID, respectively; and it is also assumed that the delimiter between sample-name and read in the read-name entries is a single "_". If this is not true, you may have to update these parameters, or even modify the current implementation of this function.

Also note that there is now a UPARSE-specific output file format, uparseout, and it might make more sense to create and import that file for use in phyloseq. If so, you'll want to import using the import_uparse() function.

See Also

import

import_biom

import_qiime

Examples

usearchfile <- system.file("extdata", "usearch.uc", package="phyloseq")
import_usearch_uc(usearchfile)

Description

A specialized function for creating a network representation of microbiomes, sample-wise or taxa-wise, based on a user-defined ecological distance and (potentially arbitrary) threshold. The graph is ultimately represented using the igraph-package.

Usage

make_network(physeq, type="samples", distance="jaccard", max.dist = 0.4, 
    keep.isolates=FALSE, ...)

Arguments

Value

A igraph-class object.

See Also

plot_network

Examples

# # Example plots with Enterotype Dataset
data(enterotype)
ig <- make_network(enterotype, max.dist=0.3)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.3, label=NULL)
# 
ig1 <- make_network(enterotype, max.dist=0.2)
plot_network(ig1, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.3, label=NULL)
# 
# # Three methods of choosing/providing distance/distance-method
# Provide method name available to distance() function
ig <- make_network(enterotype, max.dist=0.3, distance="jaccard")
# Provide distance object, already computed
jaccdist <- distance(enterotype, "jaccard")
ih <- make_network(enterotype, max.dist=0.3, distance=jaccdist)
# Provide "custom" function.
ii <- make_network(enterotype, max.dist=0.3, distance=function(x){vegan::vegdist(x, "jaccard")})
# The have equal results:		
all.equal(ig, ih)
all.equal(ig, ii)
# 
# Try out making a trivial "network" of the 3-sample esophagus data,
# with weighted-UniFrac as distance
data(esophagus)
ij <- make_network(esophagus, "samples", "unifrac", weighted=TRUE)

Description

Takes a comma-separated list of phyloseq objects as arguments, and returns the most-comprehensive single phyloseq object possible.

Usage

merge_phyloseq(...)

Arguments

Details

Higher-order objects can be created if arguments are appropriate component data types of different classes, and this should mirror the behavior of the phyloseq method, which is the suggested method if the goal is simply to create a higher-order phyloseq object from different data types (1 of each class) describing the same experiment.

By contrast, this method is intended for situations in which one wants to combine multiple higher-order objects, or multiple core component data objects (e.g. more than one otu_table) that should be combined into one object.

Merges are performed by first separating higher-order objects into a list of their component objects; then, merging any component objects of the same class into one object according to the behavior desribed in merge_phyloseq_pair; and finally, building back up a merged-object according to the constructor behavior of the phyloseq method. If the arguments contain only a single component type – several otu_table objects, for example – then a single merged object of that component type is returned.

Value

Merges are performed by first separating higher-order objects into a list of their component objects; then, merging any component objects of the same class into one object according to the behavior desribed in merge_phyloseq_pair; and finally, re-building a merged-object according to the constructor behavior of the phyloseq method. If the arguments contain only a single component type – several otu_table objects, for example – then a single merged object of the relevant component type is returned.

Merges between 2 or more tree objects are ultimately done using consensus from the ape package. This has the potential to limit somewhat the final data object, because trees don't merge with other trees in the same granular manner as data tables, and ultimately the species/taxa in higher-order phyloseq objects will be clipped to what is contained in the tree. If this an issue, the tree component should be ommitted from the argument list.

Examples

#
## # Make a random complex object
## OTU1 <- otu_table(matrix(sample(0:5,250,TRUE),25,10), taxa_are_rows=TRUE)
## tax1 <- tax_table(matrix("abc", 30, 8))
## map1 <- data.frame( matrix(sample(0:3,250,TRUE),25,10), 
##   matrix(sample(c("a","b","c"),150,TRUE), 25, 6) ) 
## map1 <- sample_data(map1)
## exam1 <- phyloseq(OTU1, map1, tax1)
## x <- exam1
## x <- phyloseq(exam1)
## y <- tax_table(exam1)
## merge_phyloseq(x, y)
## merge_phyloseq(y, y, y, y)

Description

Internal S4 methods to combine pairs of objects of classes specified in the phyloseq package. These objects must be component data of the same type (class). This is mainly an internal method, provided to illustrate how merging is performed by the more general merge_phyloseq function.

Usage

merge_phyloseq_pair(x, y)

## S4 method for signature 'otu_table,otu_table'
merge_phyloseq_pair(x, y)

## S4 method for signature 'taxonomyTable,taxonomyTable'
merge_phyloseq_pair(x, y)

## S4 method for signature 'sample_data,sample_data'
merge_phyloseq_pair(x, y)

## S4 method for signature 'phylo,phylo'
merge_phyloseq_pair(x, y)

## S4 method for signature 'XStringSet,XStringSet'
merge_phyloseq_pair(x, y)

Arguments

Details

The merge_phyloseq function is recommended in general.

Special note: non-identical trees are merged using consensus.

Value

A single component data object that matches x and y arguments. The returned object will contain the union of the species and/or samples of each. If there is redundant information between a pair of arguments of the same class, the values in x are used by default. Abundance values are summed for otu_table objects for those elements that describe the same species and sample in x and y.

See Also

merge_phyloseq merge_taxa

Examples

#
## # merge two simulated otu_table objects.
## x  <- otu_table(matrix(sample(0:5,200,TRUE),20,10), taxa_are_rows=TRUE)
## y  <- otu_table(matrix(sample(0:5,300,TRUE),30,10), taxa_are_rows=FALSE)
## xy <- merge_phyloseq_pair(x, y)
## yx <- merge_phyloseq_pair(y, x)
## # merge two simulated tax_table objects
## x <- tax_table(matrix("abc", 20, 6))
## y <- tax_table(matrix("def", 30, 8))
## xy <- merge_phyloseq_pair(x, y)
## # merge two simulated sample_data objects
## x <- data.frame( matrix(sample(0:3,250,TRUE),25,10), 
##   matrix(sample(c("a","b","c"),150,TRUE),25,6) )
## x <- sample_data(x)
## y <- data.frame( matrix(sample(4:6,200,TRUE),20,10), 
##   matrix(sample(c("d","e","f"),120,TRUE),20,8) )
## y <- sample_data(y)
## merge_phyloseq_pair(x, y)
## data.frame(merge_phyloseq_pair(x, y))
## data.frame(merge_phyloseq_pair(y, x))

Description

The purpose of this method is to merge/agglomerate the sample indices of a phyloseq object according to a categorical variable contained in a sample_data or a provided factor.

Usage

merge_samples(x, group, fun=mean)

## S4 method for signature 'sample_data'
merge_samples(x, group, fun = mean)

## S4 method for signature 'otu_table'
merge_samples(x, group)

## S4 method for signature 'phyloseq'
merge_samples(x, group, fun = mean)

Arguments

Details

NOTE: (phylo) trees and taxonomyTable-class are not modified by this function, but returned in the output object as-is.

Value

A phyloseq object that has had its sample indices merged according to the factor indicated by the group argument. The output class matches x.

See Also

merge_taxa, codemerge_phyloseq

Examples

#
data(GlobalPatterns)
GP = GlobalPatterns
mergedGP = merge_samples(GlobalPatterns, "SampleType")
SD = merge_samples(sample_data(GlobalPatterns), "SampleType")
print(SD)
print(mergedGP)
sample_names(GlobalPatterns)
sample_names(mergedGP)
identical(SD, sample_data(mergedGP))
# The OTU abundances of merged samples are summed
# Let's investigate this ourselves looking at just the top10 most abundance OTUs...
OTUnames10 = names(sort(taxa_sums(GP), TRUE)[1:10])
GP10  = prune_taxa(OTUnames10,  GP)
mGP10 = prune_taxa(OTUnames10, mergedGP)
ocean_samples = sample_names(subset(sample_data(GP), SampleType=="Ocean"))
print(ocean_samples)
otu_table(GP10)[, ocean_samples]
rowSums(otu_table(GP10)[, ocean_samples])
otu_table(mGP10)["Ocean", ]

Description

Takes as input an object that describes species/taxa (e.g. phyloseq-class, otu_table-class, phylo-class, taxonomyTable-class), as well as a vector of species that should be merged. It is intended to be able to operate at a low-level such that related methods, such as tip_glom and tax_glom can both reliably call merge_taxa for their respective purposes.

Usage

merge_taxa(x, eqtaxa, archetype=1)

## S4 method for signature 'phyloseq'
merge_taxa(x, eqtaxa,
  archetype = eqtaxa[which.max(taxa_sums(x)[eqtaxa])])

## S4 method for signature 'sample_data'
merge_taxa(x, eqtaxa, archetype = 1L)

## S4 method for signature 'otu_table'
merge_taxa(x, eqtaxa,
  archetype = eqtaxa[which.max(taxa_sums(x)[eqtaxa])])

## S4 method for signature 'phylo'
merge_taxa(x, eqtaxa, archetype = 1L)

## S4 method for signature 'XStringSet'
merge_taxa(x, eqtaxa, archetype = 1L)

## S4 method for signature 'taxonomyTable'
merge_taxa(x, eqtaxa, archetype = 1L)

Arguments

Value

The object, x, in its original class, but with the specified species merged into one entry in all relevant components.

See Also

tip_glom, tax_glom, merge_phyloseq, merge_samples

Examples

#
data(esophagus)
tree <- phy_tree(esophagus)
otu  <- otu_table(esophagus)
otutree0 <- phyloseq(otu, tree)
# plot_tree(otutree0)
otutree1 <- merge_taxa(otutree0, 1:8, 2)
# plot_tree(esophagus, ladderize="left")

Description

Originally, this function was for accessing microbiome datasets from the microbio.me/qiime public repository from within R. As you can see by clicking on the above link, the QIIME-DB sever is down indefinitely. However, this function will remain supported here in case the FTP server goes back up, and also for phyloseq users that have downloaded one or more data packages prior to the server going down.

Usage

microbio_me_qiime(zipftp, ext = ".zip",
  parsef = parse_taxonomy_greengenes, ...)

Arguments

Value

A phyloseq-class object if possible, a component if only a component could be imported, or NULL if nothing could be imported after unzipping the file. Keep in mind there is a specific naming-convention that is expected based on the current state of the microbio.me/qiime servers. Several helpful messages are catted to standard out to help let you know the ongoing status of the current download and import process.

See Also

See download.file and url for details about URL formats – including local file addresses – that might work here.

import_biom

import_qiime

Examples

# This should return TRUE on your system if you have internet turned on
# and a standard R installation. Indicates whether this is likely to
# work on your system for a URL or local file, respectively.
capabilities("http/ftp"); capabilities("fifo")
# A working example with a local example file included in phyloseq
zipfile = "study_816_split_library_seqs_and_mapping.zip"
zipfile = system.file("extdata", zipfile, package="phyloseq")
tarfile = "study_816_split_library_seqs_and_mapping.tar.gz"
tarfile = system.file("extdata", tarfile, package="phyloseq")
tarps = microbio_me_qiime(tarfile)
zipps = microbio_me_qiime(zipfile) 
identical(tarps, zipps)
tarps; zipps
plot_heatmap(tarps)
# An example that used to work, before the QIIME-DB server was turned off by its host.
# # Smokers dataset
# smokezip = "ftp://thebeast.colorado.edu/pub/QIIME_DB_Public_Studies/study_524_split_library_seqs_and_mapping.zip"
# smokers1 = microbio_me_qiime(smokezip)
# # Alternatively, just use the study number
# smokers2 = microbio_me_qiime(524)
# identical(smokers1, smokers2)

Description

Please note that it is up to you to perform any necessary normalizing / standardizing transformations prior to these tests. See for instance transform_sample_counts.

Usage

mt(physeq, classlabel, minPmaxT = "minP", method = "fdr", ...)

## S4 method for signature 'phyloseq,ANY'
mt(physeq, classlabel, minPmaxT = "minP",
  method = "fdr", ...)

## S4 method for signature 'otu_table,integer'
mt(physeq, classlabel, minPmaxT = "minP",
  method = "fdr", ...)

## S4 method for signature 'otu_table,numeric'
mt(physeq, classlabel, minPmaxT = "minP",
  method = "fdr", ...)

## S4 method for signature 'otu_table,logical'
mt(physeq, classlabel, minPmaxT = "minP",
  method = "fdr", ...)

## S4 method for signature 'otu_table,character'
mt(physeq, classlabel, minPmaxT = "minP",
  method = "fdr", ...)

## S4 method for signature 'otu_table,factor'
mt(physeq, classlabel, minPmaxT = "minP",
  method = "fdr", ...)

Arguments

Value

A dataframe with components specified in the documentation for mt.maxT or mt.minP, respectively.

See Also

mt.maxT

mt.minP

p.adjust

Examples

## # Simple example, testing genera that sig correlate with Enterotypes
data(enterotype)
# Filter samples that don't have Enterotype
x <- subset_samples(enterotype, !is.na(Enterotype))
# (the taxa are at the genera level in this dataset)
res = mt(x, "Enterotype", method="fdr", test="f", B=300)
head(res, 10)
## # Not surprisingly, Prevotella and Bacteroides top the list.
## # Different test, multiple-adjusted t-test, whether samples are ent-2 or not.
## mt(x, get_variable(x, "Enterotype")==2)

Description

Unlike, nodeplotdefault and nodeplotboot, this function does not return a function, but instead is provided directly to the nodelabf argument of plot_tree to ensure that node labels are not added to the graphic. Please note that you do not need to create or obtain the arguments to this function. Instead, you can provide this function directly to plot_tree and it will know what to do with it. Namely, use it to avoid plotting any node labels.

Usage

nodeplotblank(p, nodelabdf)

Arguments

Value

The same input object, p, provided as input. Unmodified.

See Also

nodeplotdefault

nodeplotboot

plot_tree

Examples

data("esophagus")
plot_tree(esophagus)
plot_tree(esophagus, nodelabf=nodeplotblank)

Description

Is not a labeling function itself, but returns one. The returned function is specialized for labeling bootstrap values. Note that the function that is returned has two completely different arguments from the four listed here: the plot object already built by earlier steps in plot_tree, and the data.frame that contains the relevant plotting data for the nodes (especially x, y, label), respectively. See nodeplotdefault for a simpler example. The main purpose of this and nodeplotdefault is to provide a useful default function generator for arbitrary and bootstrap node labels, respectively, and also to act as examples of functions that can successfully interact with plot_tree to add node labels to the graphic.

Usage

nodeplotboot(highthresh=95L, lowcthresh=50L, size=2L, hjust=-0.2)

Arguments

Value

A function that can add a bootstrap-values layer to the tree graphic. The values are represented in two ways; either as black filled circles indicating very high-confidence nodes, or the bootstrap value itself printed in small text next to the node on the tree.

See Also

nodeplotdefault

nodeplotblank

plot_tree

Examples

nodeplotboot()
nodeplotboot(3, -0.4)

Description

Is not a labeling function itself, but returns one. The returned function is capable of adding whatever label is on a node. Note that the function that is returned has two completely different arguments to those listed here: the plot object already built by earlier steps in plot_tree, and the data.frame that contains the relevant plotting data for the nodes (especially x, y, label), respectively. See nodeplotboot for a more sophisticated example. The main purpose of this and nodeplotboot is to provide a useful default function generator for arbitrary and bootstrap node labels, respectively, and also to act as examples of functions that will successfully interact with plot_tree to add node labels to the graphic.

Usage

nodeplotdefault(size=2L, hjust=-0.2)

Arguments

Value

A function that can add a node-label layer to a graphic.

See Also

nodeplotboot

nodeplotblank

plot_tree

Examples

nodeplotdefault()
nodeplotdefault(3, -0.4)

Description

Get the number of samples.

Usage

nsamples(physeq)

## S4 method for signature 'ANY'
nsamples(physeq)

## S4 method for signature 'phyloseq'
nsamples(physeq)

## S4 method for signature 'otu_table'
nsamples(physeq)

## S4 method for signature 'sample_data'
nsamples(physeq)

Arguments

Value

An integer indicating the total number of samples.

See Also

taxa_names, sample_names, ntaxa

Examples

#
data("esophagus")
tree <- phy_tree(esophagus)
OTU1 <- otu_table(esophagus)
nsamples(OTU1)
physeq1 <- phyloseq(OTU1, tree)
nsamples(physeq1)

Description

Get the number of taxa/species.

Usage

ntaxa(physeq)

## S4 method for signature 'ANY'
ntaxa(physeq)

## S4 method for signature 'phyloseq'
ntaxa(physeq)

## S4 method for signature 'otu_table'
ntaxa(physeq)

## S4 method for signature 'taxonomyTable'
ntaxa(physeq)

## S4 method for signature 'phylo'
ntaxa(physeq)

## S4 method for signature 'XStringSet'
ntaxa(physeq)

Arguments

Value

An integer indicating the number of taxa / species.

See Also

taxa_names

Examples

data("esophagus")
ntaxa(esophagus)
phy_tree(esophagus)
ntaxa(phy_tree(esophagus))

Description

This function wraps several commonly-used ordination methods. The type of ordination depends upon the argument to method. Try ordinate("help") or ordinate("list") for the currently supported method options.

Usage

ordinate(physeq, method = "DCA", distance = "bray", formula = NULL,
  ...)

Arguments

Value

An ordination object. The specific class of the returned object depends upon the ordination method, as well as the function/package that is called internally to perform it. As a general rule, any of the ordination classes returned by this function will be recognized by downstream tools in the phyloseq package, for example the ordination plotting function, plot_ordination.

See Also

The plot_ordination Tutorial

Related component ordination functions described within phyloseq:

DPCoA

Described/provided by other packages:

cca/rda, decorana, metaMDS, pcoa, capscale

NMDS and MDS/PCoA both operate on distance matrices, typically based on some pairwise comparison of the microbiomes in an experiment/project. There are a number of common methods to use to calculate these pairwise distances, and the most convenient function (from a phyloseq point of view) for calculating these distance matrices is the

distance

function. It can be thought of as a distance / dissimilarity-index companion function for ordinate, and indeed the distance options provided to ordinate are often simply passed on to distance.

A good quick summary of ordination is provided in the introductory vignette for vegan:

vegan introductory vignette

The following R task views are also useful for understanding the available tools in R:

Analysis of Ecological and Environmental Data

Multivariate Statistics

Examples

# See http://joey711.github.io/phyloseq/plot_ordination-examples
# for many more examples.
# plot_ordination(GP, ordinate(GP, "DCA"), "samples", color="SampleType")

Description

This is the suggested method for both constructing and accessing Operational Taxonomic Unit (OTU) abundance (otu_table-class) objects. When the first argument is a matrix, otu_table() will attempt to create and return an otu_table-class object, which further depends on whether or not taxa_are_rows is provided as an additional argument. Alternatively, if the first argument is an experiment-level (phyloseq-class) object, then the corresponding otu_table is returned.

Usage

otu_table(object, taxa_are_rows, errorIfNULL=TRUE)

## S4 method for signature 'phyloseq'
otu_table(object, errorIfNULL = TRUE)

## S4 method for signature 'otu_table'
otu_table(object, errorIfNULL = TRUE)

## S4 method for signature 'matrix'
otu_table(object, taxa_are_rows)

## S4 method for signature 'data.frame'
otu_table(object, taxa_are_rows)

## S4 method for signature 'ANY'
otu_table(object, errorIfNULL = TRUE)

Arguments

Value

An otu_table-class object.

See Also

phy_tree, sample_data, tax_table phyloseq, merge_phyloseq

Examples

#
# data(GlobalPatterns)
# otu_table(GlobalPatterns)

Description

Because orientation of these tables can vary by method, the orientation is defined explicitly in the taxa_are_rows slot (a logical). The otu_table class inherits the matrix class to store abundance values. Various standard subset and assignment nomenclature has been extended to apply to the otu_table class, including square-bracket, t, etc.

Details

taxa_are_rows

A single logical specifying the orientation of the abundance table.

.Data

This slot is inherited from the matrix class.

Description

Assign a new OTU Table to x

Usage

otu_table(x) <- value

## S4 replacement method for signature 'phyloseq,otu_table'
otu_table(x) <- value

## S4 replacement method for signature 'otu_table,otu_table'
otu_table(x) <- value

## S4 replacement method for signature 'phyloseq,phyloseq'
otu_table(x) <- value

Arguments

Examples

# data(GlobalPatterns)
# # An example of pruning to just the first 100 taxa in GlobalPatterns.
# ex2a <- prune_taxa(taxa_names(GlobalPatterns)[1:100], GlobalPatterns)
# # The following 3 lines produces an ex2b that is equal to ex2a
# ex2b <- GlobalPatterns
# OTU <- otu_table(GlobalPatterns)[1:100, ]
# otu_table(ex2b) <- OTU
# identical(ex2a, ex2b)
# print(ex2b)
# # Relace otu_table by implying the component in context.
# ex2c <- GlobalPatterns
# otu_table(ex2c) <- ex2b
# identical(ex2a, ex2c)

Description

These are provided as both example and default functions for parsing a character vector of taxonomic rank information for a single taxa. As default functions, these are intended for cases where the data adheres to the naming convention used by greengenes (http://greengenes.lbl.gov/cgi-bin/nph-index.cgi) or where the convention is unknown, respectively. To work, these functions – and any similar custom function you may want to create and use – must take as input a single character vector of taxonomic ranks for a single OTU, and return a named character vector that has been modified appropriately (according to known naming conventions, desired length limits, etc. The length (number of elements) of the output named vector does not need to be equal to the input, which is useful for the cases where the source data files have extra meaningless elements that should probably be removed, like the ubiquitous “Root” element often found in greengenes/QIIME taxonomy labels. In the case of parse_taxonomy_default, no naming convention is assumed and so dummy rank names are added to the vector. More usefully if your taxonomy data is based on greengenes, the parse_taxonomy_greengenes function clips the first 3 characters that identify the rank, and uses these to name the corresponding element according to the appropriate taxonomic rank name used by greengenes (e.g. "p__" at the beginning of an element means that element is the name of the phylum to which this OTU belongs). Most importantly, the expectations for these functions described above make them compatible to use during data import, specifcally the import_biom function, but it is a flexible structure that will be implemented soon for all phyloseq import functions that deal with taxonomy (e.g. import_qiime).

Usage

parse_taxonomy_default(char.vec)

parse_taxonomy_greengenes(char.vec)

parse_taxonomy_qiime(char.vec)

Arguments

Value

A character vector in which each element is a different taxonomic rank of the same OTU, and each element name is the name of the rank level. For example, an element might be "Firmicutes" and named "phylum". These parsed, named versions of the taxonomic vector should reflect embedded information, naming conventions, desired length limits, etc; or in the case of parse_taxonomy_default, not modified at all and given dummy rank names to each element.

See Also

import_biom import_qiime

Examples

taxvec1 = c("Root", "k__Bacteria", "p__Firmicutes", "c__Bacilli", "o__Bacillales", "f__Staphylococcaceae")
 parse_taxonomy_default(taxvec1)
 parse_taxonomy_greengenes(taxvec1)
 taxvec2 = c("Root;k__Bacteria;p__Firmicutes;c__Bacilli;o__Bacillales;f__Staphylococcaceae")
 parse_taxonomy_qiime(taxvec2)

Retrieve phylogenetic tree (phylo-class) from object.

Description

This is the suggested method for accessing the phylogenetic tree, (phylo-class) from a phyloseq-class. Like other accessors (see See Also, below), the default behavior of this method is to stop with an error if physeq is a phyloseq-class but does not contain a phylogenetic tree (the component data you are trying to access in this case).

Usage

phy_tree(physeq, errorIfNULL=TRUE)

## S4 method for signature 'ANY'
phy_tree(physeq, errorIfNULL = TRUE)

## S4 method for signature 'phylo'
phy_tree(physeq)

Arguments

Details

Note that the tip labels should be named to match the taxa_names of the other objects to which it is going to be paired. The phyloseq constructor automatically checks for exact agreement in the set of species described by the phlyogenetic tree and the other components (taxonomyTable, otu_table), and trims as-needed. Thus, the tip.labels in a phylo object must be named to match the results of taxa_names of the other objects to which it will ultimately be paired.

Value

The phylo-class object contained within physeq; or NULL if physeq does not have a tree. This method stops with an error in the latter NULL case be default, which can be over-ridden by changing the value of errorIfNULL to FALSE.

See Also

otu_table, sample_data, tax_table refseq, phyloseq, merge_phyloseq

Examples

data(GlobalPatterns)
 phy_tree(GlobalPatterns)

Description

Assign a (new) phylogenetic tree to x

Usage

phy_tree(x) <- value

## S4 replacement method for signature 'phyloseq,phylo'
phy_tree(x) <- value

## S4 replacement method for signature 'phyloseq,phyloseq'
phy_tree(x) <- value

Arguments

Examples

#
data("esophagus")
# An example of pruning to just the first 20 taxa in esophagus
ex2a <- prune_taxa(taxa_names(esophagus)[1:20], esophagus)
# The following 3 lines produces an ex2b that is equal to ex2a
ex2b <- ex2a
phy_tree(ex2b) <- phy_tree(esophagus)
identical(ex2a, ex2b)

Description

See the ape package for details about this type of representation of a phylogenetic tree. It is used throughout the ape package.

See Also

phylo, setOldClass

Description

phyloseq() is a constructor method, This is the main method suggested for constructing an experiment-level (phyloseq-class) object from its component data (component data classes: otu_table-class, sample_data-class, taxonomyTable-class, phylo-class).

Usage

phyloseq(...)

Arguments

Value

The class of the returned object depends on the argument class(es). For an experiment-level object, two or more component data objects must be provided. Otherwise, if a single component-class is provided, it is simply returned as-is. The order of arguments does not matter.

See Also

merge_phyloseq

Examples

data(esophagus)
x1 = phyloseq(otu_table(esophagus), phy_tree(esophagus))
identical(x1, esophagus)
# # data(GlobalPatterns)
# # GP <- GlobalPatterns
# # phyloseq(sample_data(GP), otu_table(GP))
# # phyloseq(otu_table(GP), phy_tree(GP))
# # phyloseq(tax_table(GP), otu_table(GP))
# # phyloseq(phy_tree(GP), otu_table(GP), sample_data(GP))
# # phyloseq(otu_table(GP), tax_table(GP), sample_data(GP))
# # phyloseq(otu_table(GP), phy_tree(GP), tax_table(GP), sample_data(GP))

Description

No testing is performed by this function. The phyloseq data is converted to the relevant DESeqDataSet object, which can then be tested in the negative binomial generalized linear model framework of the DESeq function in DESeq2 package. See the phyloseq-extensions tutorials for more details.

Usage

phyloseq_to_deseq2(physeq, design, ...)

Arguments

Value

A DESeqDataSet object.

See Also

vignette("phyloseq-mixture-models")

The phyloseq-extensions tutorials.

DESeq

results

DESeqDataSetFromMatrix

Examples

# Check out the vignette phyloseq-mixture-models for more details.
 # vignette("phyloseq-mixture-models")
 data(soilrep)
 phyloseq_to_deseq2(soilrep, ~warmed)

Description

No testing is performed by this function. The phyloseq data is converted to the relevant MRexperiment-class object, which can then be tested in the zero-inflated mixture model framework (e.g. fitZig) in the metagenomeSeq package. See the phyloseq-extensions tutorials for more details.

Usage

phyloseq_to_metagenomeSeq(physeq, ...)

Arguments

Value

A MRexperiment-class object.

See Also

fitTimeSeries fitLogNormal fitZig MRtable MRfulltable

Examples

# Check out the vignette metagenomeSeq for more details.
 # vignette("metagenomeSeq")
 data(soilrep)
 phyloseq_to_metagenomeSeq(soilrep)

Description

Contains all currently-supported component data classes: otu_table-class, sample_data-class, taxonomyTable-class ("tax_table" slot), phylo-class ("phy_tree" slot), and the XStringSet-class ("refseq" slot). There are several advantages to storing your phylogenetic sequencing experiment as an instance of the phyloseq class, not the least of which is that it is easy to return to the data later and feel confident that the different data types “belong” to one another. Furthermore, the phyloseq constructor ensures that the different data components have compatible indices (e.g. OTUs and samples), and performs the necessary trimming automatically when you create your “experiment-level” object. Downstream analyses are aware of which data classes they require – and where to find them – often making your phyloseq-class object the only data argument required for analysis and plotting functions (although there are many options and parameter arguments available to you).

Details

In the case of missing component data, the slots are set to NULL. As soon as a phyloseq-class object is to be updated with new component data (previously missing/NULL or not), the indices of all components are re-checked for compatibility and trimmed if necessary. This is to ensure by design that components describe the same taxa/samples, and also that these trimming/validity checks do not need to be repeated in downstream analyses.

slots:

otu_table

a single object of class otu_table.

sam_data

a single object of class sample_data.

tax_table

a single object of class taxonomyTable.

phy_tree

a single object of the phylo-class, from the ape package.

refseq

a biological sequence set object of a class that inherits from the XStringSet-class, from the Biostrings package.

See Also

The constructor, phyloseq, the merger merge_phyloseq, and also the component constructor/accessors otu_table, sample_data, tax_table, phy_tree, and refseq.

Description

These will be migrated to "defunct" status in the next release, and removed completely in the release after that. These functions are provided for compatibility with older version of the phyloseq package. They may eventually be completely removed.

Usage

deprecated_phyloseq_function(x, value, ...)

Arguments

Details

Description

There are many useful examples of phyloseq barplot graphics in the phyloseq online tutorials. This function wraps ggplot2 plotting, and returns a ggplot2 graphic object that can be saved or further modified with additional layers, options, etc. The main purpose of this function is to quickly and easily create informative summary graphics of the differences in taxa abundance between samples in an experiment.

Usage

plot_bar(physeq, x="Sample", y="Abundance", fill=NULL,
 title=NULL, facet_grid=NULL)

Arguments

Value

A ggplot2 graphic object – rendered in the graphical device as the default print/show method.

See Also

phyloseq online tutorials.

psmelt

ggplot

qplot

Examples

data("GlobalPatterns")
gp.ch = subset_taxa(GlobalPatterns, Phylum == "Chlamydiae")
plot_bar(gp.ch)
plot_bar(gp.ch, fill="Genus")
plot_bar(gp.ch, x="SampleType", fill="Genus")
plot_bar(gp.ch, "SampleType", fill="Genus", facet_grid=~Family)
# See additional examples in the plot_bar online tutorial. Link above.

Description

Create a ggplot summary of gap statistic results

Usage

plot_clusgap(clusgap, title = "Gap Statistic results")

Arguments

Value

A ggplot plot object. The rendered graphic should be a plot of the gap statistic score versus values for k, the number of clusters.

See Also

gapstat_ord

clusGap

ggplot

Examples

# Load and process data
data("soilrep")
soilr = rarefy_even_depth(soilrep, rngseed=888)
print(soilr)
sample_variables(soilr)
# Ordination
sord  = ordinate(soilr, "DCA")
# Gap Statistic
gs = gapstat_ord(sord, axes=1:4, verbose=FALSE)
# Evaluate results with plots, etc.
plot_scree(sord)
plot_ordination(soilr, sord,  color="Treatment")
plot_clusgap(gs)
print(gs, method="Tibs2001SEmax")
# Non-ordination example, use cluster::clusGap function directly
library("cluster")
pam1 = function(x, k){list(cluster = pam(x, k, cluster.only=TRUE))}
gs.pam.RU = clusGap(ruspini, FUN = pam1, K.max = 8, B = 60)
gs.pam.RU
plot(gs.pam.RU, main = "Gap statistic for the 'ruspini' data")
mtext("k = 4 is best .. and  k = 5  pretty close")
plot_clusgap(gs.pam.RU)

Description

There are many useful examples of phyloseq heatmap graphics in the phyloseq online tutorials. In a 2010 article in BMC Genomics, Rajaram and Oono show describe an approach to creating a heatmap using ordination methods to organize the rows and columns instead of (hierarchical) cluster analysis. In many cases the ordination-based ordering does a much better job than h-clustering. An immediately useful example of their approach is provided in the NeatMap package for R. The NeatMap package can be used directly on the abundance table (otu_table-class) of phylogenetic-sequencing data, but the NMDS or PCA ordination options that it supports are not based on ecological distances. To fill this void, phyloseq provides the plot_heatmap() function as an ecology-oriented variant of the NeatMap approach to organizing a heatmap and build it using ggplot2 graphics tools. The distance and method arguments are the same as for the plot_ordination function, and support large number of distances and ordination methods, respectively, with a strong leaning toward ecology. This function also provides the options to re-label the OTU and sample axis-ticks with a taxonomic name and/or sample variable, respectively, in the hope that this might hasten your interpretation of the patterns (See the sample.label and taxa.label documentation, below). Note that this function makes no attempt to overlay hierarchical clustering trees on the axes, as hierarchical clustering is not used to organize the plot. Also note that each re-ordered axis repeats at the edge, and so apparent clusters at the far right/left or top/bottom of the heat-map may actually be the same. For now, the placement of this edge can be considered arbitrary, so beware of this artifact of this graphical representation. If you benefit from this phyloseq-specific implementation of the NeatMap approach, please cite both our packages/articles.

Usage

plot_heatmap(physeq, method = "NMDS", distance = "bray",
  sample.label = NULL, taxa.label = NULL, low = "#000033",
  high = "#66CCFF", na.value = "black", trans = log_trans(4),
  max.label = 250, title = NULL, sample.order = NULL,
  taxa.order = NULL, first.sample = NULL, first.taxa = NULL, ...)

Arguments

Details

This approach borrows heavily from the heatmap1 function in the NeatMap package. Highly recommended, and we are grateful for their package and ideas, which we have adapted for our specific purposes here, but did not use an explicit dependency. At the time of the first version of this implementation, the NeatMap package depends on the rgl-package, which is not needed in phyloseq, at present. Although likely a transient issue, the rgl-package has some known installation issues that have further influenced to avoid making NeatMap a formal dependency (Although we love both NeatMap and rgl!).

Value

A heatmap plot, in the form of a ggplot2 plot object, which can be further saved and modified.

References

Because this function relies so heavily in principle, and in code, on some of the functionality in NeatMap, please site their article if you use this function in your work.

Rajaram, S., & Oono, Y. (2010). NeatMap–non-clustering heat map alternatives in R. BMC Bioinformatics, 11, 45.

Please see further examples in the phyloseq online tutorials.

Examples

data("GlobalPatterns")
gpac <- subset_taxa(GlobalPatterns, Phylum=="Crenarchaeota")
# FYI, the base-R function uses a non-ecological ordering scheme,
# but does add potentially useful hclust dendrogram to the sides...
gpac <- subset_taxa(GlobalPatterns, Phylum=="Crenarchaeota")
# Remove the nearly-empty samples (e.g. 10 reads or less)
gpac = prune_samples(sample_sums(gpac) > 50, gpac)
# Arbitrary order if method set to NULL
plot_heatmap(gpac, method=NULL, sample.label="SampleType", taxa.label="Family")
# Use ordination
plot_heatmap(gpac, sample.label="SampleType", taxa.label="Family")
# Use ordination for OTUs, but not sample-order
plot_heatmap(gpac, sample.label="SampleType", taxa.label="Family", sample.order="SampleType")
# Specifying both orders omits any attempt to use ordination. The following should be the same.
p0 = plot_heatmap(gpac, sample.label="SampleType", taxa.label="Family", taxa.order="Phylum", sample.order="SampleType")
p1 = plot_heatmap(gpac, method=NULL, sample.label="SampleType", taxa.label="Family", taxa.order="Phylum", sample.order="SampleType")
#expect_equivalent(p0, p1)
# Example: Order matters. Random ordering of OTU indices is difficult to interpret, even with structured sample order
rando = sample(taxa_names(gpac), size=ntaxa(gpac), replace=FALSE)
plot_heatmap(gpac, method=NULL, sample.label="SampleType", taxa.label="Family", taxa.order=rando, sample.order="SampleType")
# # Select the edges of each axis. 
# First, arbitrary edge, ordering
plot_heatmap(gpac, method=NULL)
# Second, biological-ordering (instead of default ordination-ordering), but arbitrary edge
plot_heatmap(gpac, taxa.order="Family", sample.order="SampleType")
# Third, biological ordering, selected edges
plot_heatmap(gpac, taxa.order="Family", sample.order="SampleType", first.taxa="546313", first.sample="NP2")
# Fourth, add meaningful labels
plot_heatmap(gpac, sample.label="SampleType", taxa.label="Family", taxa.order="Family", sample.order="SampleType", first.taxa="546313", first.sample="NP2")

Description

There are many useful examples of phyloseq network graphics in the phyloseq online tutorials. A custom plotting function for displaying networks using advanced ggplot2 formatting. Note that this function is a performance and interface revision to plot_network, which requires an igraph object as its first argument. This new function is more in-line with other plot_* functions in the phyloseq-package, in that its first/main argument is a phyloseq-class instance. Edges in the network are created if the distance between nodes is below a (potentially arbitrary) threshold, and special care should be given to considering the choice of this threshold. However, network line thickness and opacity is scaled according to the similarity of vertices (either samples or taxa), helping to temper, somewhat, the effect of the threshold. Also note that the choice of network layout algorithm can have a large effect on the impression and interpretability of the network graphic, and you may want to familiarize yourself with some of these options (see the laymeth argument).

Usage

plot_net(physeq, distance = "bray", type = "samples", maxdist = 0.7,
  laymeth = "fruchterman.reingold", color = NULL, shape = NULL,
  rescale = FALSE, point_size = 5, point_alpha = 1,
  point_label = NULL, hjust = 1.35, title = NULL)

Arguments

Value

A ggplot2 network plot. Will render to default graphic device automatically as print side effect. Can also be saved, further manipulated, or rendered to a vector or raster file using ggsave.

See Also

Original network plotting functions:

make_network

plot_network

Examples

data(enterotype)
plot_net(enterotype, color="SeqTech", maxdist = 0.3)
plot_net(enterotype, color="SeqTech", maxdist = 0.3, laymeth = "auto")
plot_net(enterotype, color="SeqTech", maxdist = 0.3, laymeth = "svd")
plot_net(enterotype, color="SeqTech", maxdist = 0.3, laymeth = "circle")
plot_net(enterotype, color="SeqTech", shape="Enterotype", maxdist = 0.3, laymeth = "circle")

Description

There are many useful examples of phyloseq network graphics in the phyloseq online tutorials. A custom plotting function for displaying networks using advanced ggplot2 formatting. The network itself should be represented using the igraph package. For the phyloseq-package it is suggested that the network object (argument g) be created using the make_network function, and based upon sample-wise or taxa-wise microbiome ecological distances calculated from a phylogenetic sequencing experiment (phyloseq-class). In this case, edges in the network are created if the distance between nodes is below a potentially arbitrary threshold, and special care should be given to considering the choice of this threshold.

Usage

plot_network(g, physeq=NULL, type="samples", 
	color=NULL, shape=NULL, point_size=4, alpha=1,
	label="value", hjust = 1.35, 
	line_weight=0.5, line_color=color, line_alpha=0.4,
	layout.method=layout.fruchterman.reingold, title=NULL)

Arguments

Value

A ggplot2 plot representing the network, with optional mapping of variable(s) to point color or shape.

References

This code was adapted from a repo original hosted on GitHub by Scott Chamberlain: https://github.com/SChamberlain/gggraph

The code most directly used/modified was first posted here: http://www.r-bloggers.com/basic-ggplot2-network-graphs/

See Also

make_network

Examples

data(enterotype)
ig <- make_network(enterotype, max.dist=0.3)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.3, label=NULL)
# Change distance parameter
ig <- make_network(enterotype, max.dist=0.2)
plot_network(ig, enterotype, color="SeqTech", shape="Enterotype", line_weight=0.3, label=NULL)

Description

There are many useful examples of phyloseq ordination graphics in the phyloseq online tutorials. Convenience wrapper for plotting ordination results as a ggplot2-graphic, including additional annotation in the form of shading, shape, and/or labels of sample variables.

Usage

plot_ordination(physeq, ordination, type = "samples", axes = 1:2,
  color = NULL, shape = NULL, label = NULL, title = NULL,
  justDF = FALSE)

Arguments

Value

A ggplot plot object, graphically summarizing the ordination result for the specified axes.

See Also

Many more examples are included in the phyloseq online tutorials.

Also see the general wrapping function:

plot_phyloseq

Examples

# See other examples at
# http://joey711.github.io/phyloseq/plot_ordination-examples
data(GlobalPatterns)
GP = prune_taxa(names(sort(taxa_sums(GlobalPatterns), TRUE)[1:50]), GlobalPatterns)
gp_bray_pcoa = ordinate(GP, "CCA", "bray")
plot_ordination(GP, gp_bray_pcoa, "samples", color="SampleType")

Description

There are many useful examples of phyloseq graphics functions in the phyloseq online tutorials. The specific plot type is chosen according to available non-empty slots. This is mainly for syntactic convenience and quick-plotting. See links below for some examples of available graphics tools available in the phyloseq-package.

Usage

plot_phyloseq(physeq, ...)

## S4 method for signature 'phyloseq'
plot_phyloseq(physeq, ...)

Arguments

Value

A plot is created. The nature and class of the plot depends on the physeq argument, specifically, which component data classes are present.

See Also

phyloseq frontpage tutorials.

plot_ordination plot_heatmap plot_tree plot_network plot_bar plot_richness

Examples

data(esophagus)
plot_phyloseq(esophagus)

Description

There are many useful examples of alpha-diversity graphics in the phyloseq online tutorials. This function estimates a number of alpha-diversity metrics using the estimate_richness function, and returns a ggplot plotting object. The plot generated by this function will include every sample in physeq, but they can be further grouped on the horizontal axis through the argument to x, and shaded according to the argument to color (see below). You must use untrimmed, non-normalized count data for meaningful results, as many of these estimates are highly dependent on the number of singletons. You can always trim the data later on if needed, just not before using this function.

Usage

plot_richness(physeq, x = "samples", color = NULL, shape = NULL,
  title = NULL, scales = "free_y", nrow = 1, shsi = NULL,
  measures = NULL, sortby = NULL)

Arguments

Details

NOTE: Because this plotting function incorporates the output from estimate_richness, the variable names of that output should not be used as x or color (even if it works, the resulting plot might be kindof strange, and not the intended behavior of this function). The following are the names you will want to avoid using in x or color:

c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher").

Value

A ggplot plot object summarizing the richness estimates, and their standard error.

See Also

estimate_richness

estimateR

diversity

There are many more interesting examples at the phyloseq online tutorials.

Examples

## There are many more interesting examples at the phyloseq online tutorials.
## http://joey711.github.io/phyloseq/plot_richness-examples
data("soilrep")
plot_richness(soilrep, measures=c("InvSimpson", "Fisher"))
plot_richness(soilrep, "Treatment", "warmed", measures=c("Chao1", "ACE", "InvSimpson"), nrow=3)
data("GlobalPatterns")
plot_richness(GlobalPatterns, x="SampleType", measures=c("InvSimpson"))
plot_richness(GlobalPatterns, x="SampleType", measures=c("Chao1", "ACE", "InvSimpson"), nrow=3)
plot_richness(GlobalPatterns, x="SampleType", measures=c("Chao1", "ACE", "InvSimpson"), nrow=3, sortby = "Chao1")