Data Import

Getting your GDC count table, map file and tree into a phyloseq object, and how to work with it.

phyloseq provides a convenient way to organise amplicon sequencing data in R. Rather than working with separate count tables, taxonomy tables, metadata, and trees, it stores them together in a single object that many analysis packages understand, so that the sample names and OTU/zOTU IDs across all of them remain synchronized. Once your data is a phyloseq object, most downstream analyses, diversity metrics, ordinations, subsetting, plotting, agglomerating by taxonomic rank, operate directly on it without you having to manually re-align matrices by hand every time.

In practice, this means a phyloseq object behaves like a single container. Instead of passing separate count tables, taxonomy tables, and metadata into every analysis function individually, you generally pass just ps, and the function extracts whatever components it actually needs from inside it.

This page is specifically about importing GDC output, structured the way our pipeline produces it (see Data Prep Output). If your data comes from QIIME2, DADA2, or mothur instead, the same general phyloseq concepts apply, but the specific import steps here won't match your files.

What You Need

See the Data Prep Output page for the full explanation of these files. In short:

The annotated count table (*_Count_Sintax.txt, or *_Count_<reference name>.txt), required
The map (metadata) file, required
The tree file (*.tre), optional but recommended if you'll use UniFrac or other phylogeny-aware methods
The FASTA file (*.fa.gz), optional

Each file maps onto one phyloseq component:

Count table   →  otu_table()
Taxonomy      →  tax_table()
Map file      →  sample_data()
Tree file     →  phy_tree()   (via read_tree())
FASTA file    →  refseq()     (via readDNAStringSet())
                       ↓
         phyloseq()  combines them all

phyloseq() doesn't just bundle these into a list. It builds an S4 object, R's more formal, structured class system, with each component stored in its own defined slot. That's what makes the whole accessor system work the way it does: otu_table(ps), tax_table(ps), and the rest aren't just convenient shorthand, they're S4 generic functions that know how to reach into the right slot and, in some cases, check that what's there is actually valid (for example, that sample names in the count table and the metadata line up) before handing it back to you. It's also why you generally interact with a phyloseq object through these accessor functions rather than by poking at its internal structure directly.

Standard Workflow

If you generated your files with the GDC pipeline, this is all you need. Everything after this section is for special cases, unusual file layouts, or extra control, safe to skip until you actually need it.

Although the name comes from QIIME1, import_qiime() works well directly against the specific format our pipeline produces: it correctly links the count table, taxonomy, and metadata in our _Count_<reference>.txt and map files, and it's what we actually use day to day. Treat this as a pipeline-specific recommendation for GDC output, not a general phyloseq recommendation, the function's name will otherwise suggest it's only for QIIME.

Our annotated count table already carries taxonomy for every (z)OTU inline, so a single call gets you everything at once, count table, taxonomy, map file, tree, and reference sequences:

### Load Data ----
otu  <- "<projectNumber>_<runID>_<locus>_ZOTU_Count_<reference>.txt"
map  <- "<projectNumber>_<runID>_<locus>_MapFile.txt"
tree <- "<projectNumber>_<runID>_<locus>_ZOTU_MSA.tre"
ref  <- "<projectNumber>_<runID>_<locus>_ZOTU.fa"

ps <- import_qiime(otufilename = otu, mapfilename = map, treefilename = tree, refseqfilename = ref)
ps  # e.g. 144 samples / 6773 zOTU

## Explore Data ----

<reference> is whatever reference or classifier produced the annotation, for example Sintax, Silva138, or UNITE10 (see the Data Prep Output page's note on filenames). Drop tree = and/or ref = from the call if you don't have those files.

Worth a Quick Taxonomy Check

import_qiime() reporting the right number of samples and zOTUs confirms the count table came in correctly, but it's worth a separate glance at tax_table(ps) specifically, to make sure taxonomy parsed into proper rank columns rather than one long unparsed string. If it looks wrong for your specific file, the manual construction method further down gives you full control over how taxonomy gets split.

Save the Object, Don't Repeat the Import

Importing and validating your data is normally a one-time cost. Save the object once, and every later session can pick straight back up without redoing the import or re-running your validation checks:

saveRDS(ps, file = "phyloseq_object.rds")
ps <- readRDS("phyloseq_object.rds")

Check That It Imported Correctly

Before doing anything else, check the basics. This is the point where import problems are easiest to catch, before they've propagated into every downstream result. The most common failures aren't wrong dimensions, they're mismatched IDs: hidden whitespace in a sample name, or metadata rows that don't actually line up with your count table.

ps                       # summary: dimensions of each component
ntaxa(ps)                # does this match your expected number of OTUs/zOTUs?
nsamples(ps)             # does this match your number of samples?
sample_sums(ps)[1:5]     # do these look like plausible read depths?
rank_names(ps)           # are your taxonomic ranks actually present, and not all NA?
sample_names(ps)[1:5]    # do these look like your actual sample IDs, no stray whitespace?
taxa_names(ps)[1:5]      # do these look like your actual OTU/zOTU IDs?

If any of these look off, work backwards: check sample_names(ps) and taxa_names(ps) against your original files for mismatched IDs.

Exploring the phyloseq Object

# Phyloseq object
summary(ps)
# Summarize phyloseq object
microbiome::summarize_phyloseq(ps)

Accessors

Components of a phyloseq object can be extracted using accessor functions. These are the ones you'll actually reach for day to day:

otu_table(ps)         # count table
tax_table(ps)         # taxonomy table
sample_data(ps)       # metadata
ntaxa(ps)             # number of OTUs
nsamples(ps)          # number of samples
taxa_names(ps)        # OTU labels
sample_names(ps)      # sample IDs
sample_sums(ps)       # counts per sample, i.e. sequencing depth

Additional Accessors

Used less often, but still available directly:

taxa_sums(ps)         # count per OTU, summed across samples
rank_names(ps)        # taxonomic rank names available
sample_variables(ps)  # variables defined in the map file
get_variable(ps)      # metadata

get_taxa() and get_sample() Are Named Confusingly

This one trips people up, and it's phyloseq's own naming, not something to work around: get_sample(physeq, i) takes a taxon name and returns its abundance across every sample, while get_taxa(physeq, i) takes a sample name and returns every taxon's abundance in that sample. The function called "get_taxa" is the one you feed a sample name to, and vice versa. Worth double-checking the documentation (?get_taxa) the first few times you use these, since the names alone will lead you the wrong way.

get_taxa(ps, i = "SB1")     # count table: all OTU counts for one sample
get_sample(ps, i = "OTU1")  # count table: one OTU's counts across all samples

For Experienced Users: Accessing Internal Slots Directly

access() returns the component object specified by slot, or NULL if that slot doesn't exist on the object. It's a lower-level accessor than the convenience functions above, more often used internally by other functions than typed directly.

access(ps, "otu_table")
access(ps, "tax_table")
access(ps, "phy_tree")

Access Examples

# Library size / sequencing depth for specific samples
sample_sums(ps)[c("SB1", "SB2")]

# Total abundance of specific OTU(s)
taxa_sums(ps)[c("OTU1", "OTU15")]

# Counts per sample for a specific OTU
get_sample(ps, i = "OTU1")

# Counts per OTU for specific samples
get_taxa(ps, i = c("SB2", "SC2"))

# Values for specific metadata variables
get_variable(ps, varName = c("pH", "Temp"))

Modify via Accessor

# Add a derived variable to the sample metadata
sample_data(ps)$NewVariable <- ifelse(sample_data(ps)$Temp >= 21, "High", "Low")
sample_data(ps)
sample_variables(ps)

# Edit individual values directly
otu_table(ps)["OTU15", "SC2"] <- 0
tax_table(ps)["OTU7", "Genus"] <- "Clostridiales"

Editing Values Directly Is a Last Resort

Editing counts or taxonomy by hand, as in the last two lines above, bypasses whatever process originally generated those values and won't be reflected anywhere else, your report files, your raw data, or anyone else's copy of the analysis. The same principle we mentioned on the Get Ready page applies here too: if a value is wrong, fix it at the source and rebuild the phyloseq object, rather than patching the object itself. Direct editing is occasionally the right tool for a one-off, deliberate, well-documented correction, but if you find yourself doing it repeatedly, that's a sign the fix belongs further upstream.

Exporting to Plain Files

write_phyloseq() isn't part of phyloseq itself, it comes from the microbiome package, and writes the OTU table, taxonomy, and metadata out as separate, plain CSV files, not the phyloseq object itself, unlike saveRDS() above. Useful for sharing a snapshot with a collaborator who doesn't work in R, or for archiving alongside a publication, since CSV is about as universally readable as it gets, but it won't reload as a phyloseq object the way saveRDS()/readRDS() will.

library(microbiome)
write_phyloseq(ps, "OTU")
write_phyloseq(ps, "TAXONOMY")
write_phyloseq(ps, "METADATA")

Advanced / Special Cases

Everything below is for situations the standard workflow above doesn't cover: a second taxonomy, a differently laid-out file, or wanting full manual control over each component. Skip this section unless you have a specific reason to need it.

Adding Taxonomy as a Separate Step

Useful when you want to attach a taxonomy that didn't come bundled with the count table, for example a second, independently run classification against a different reference. Start from the count table without annotation instead:

otu_plain <- "<projectNumber>_<runID>_<locus>_ZOTU_Count.txt"
ps <- import_qiime(otufilename = otu_plain, mapfilename = map, treefilename = tree, refseqfilename = ref)

Then build a taxonomy table separately, from a .tax file or any other source, and attach it. It needs one row per OTU/zOTU matching taxa_names(ps), and one column per rank:

tax_df <- read.delim("<projectNumber>_<runID>_<locus>_ZOTU.tax", row.names = 1)
tax_table(ps) <- tax_table(as.matrix(tax_df))

Not Yet Verified Against Real Data

We haven't tested hands-on whether this assignment matches rows to taxa_names(ps) by name automatically, or requires them to already be in the same order. Treat this as a starting point rather than a verified recipe, and check tax_table(ps) carefully afterward before trusting it.

Attaching More Than One Taxonomy, or Tree

It's possible to attach more than one taxonomy (for example, classifications from two different reference databases) or more than one tree to the same phyloseq object, but the details take more space than fits here, and we'd rather not hand you an untested code example for something this easy to get subtly wrong. We'll cover this properly on a dedicated page later. If you need this now, ask us.

Build a phyloseq Object Manually

If your file layout doesn't match what import_qiime() expects, for instance a differently structured taxonomy file, you can build the object component by component instead and combine them with the phyloseq() constructor, following the same mapping shown near the top of this page.

library(phyloseq)

## 1. Count table, without annotation
otu_file <- "<projectNumber>_<runID>_<locus>_ZOTU_Count.txt"
otu_mat  <- as.matrix(read.delim(otu_file, row.names = 1, check.names = FALSE))
OTU      <- otu_table(otu_mat, taxa_are_rows = TRUE)

## Taxonomy, from its own file. See "Handling Taxonomy" below, this is the
## one step you'll need to adapt to your file's actual layout.
tax_file <- "<projectNumber>_<runID>_<locus>_ZOTU.tax"
tax_df   <- read.delim(tax_file, row.names = 1)
TAX      <- tax_table(as.matrix(tax_df))

## 2. Map (metadata) file
map_file <- "<projectNumber>_<runID>_<locus>_MapFile.txt"
map      <- read.delim(map_file, row.names = 1)
MAP      <- sample_data(map)

## 3. Tree file (optional)
tree_file <- "<projectNumber>_<runID>_<locus>_ZOTU_MSA.tre"
TREE      <- read_tree(tree_file)   # phyloseq's wrapper around ape::read.tree

## 4. Reference sequences (optional, see caution below before including this)
library(Biostrings)
ref_file <- "<projectNumber>_<runID>_<locus>_ZOTU.fa"
REFSEQ   <- readDNAStringSet(ref_file)

## Assemble everything into one phyloseq object.
## Omit TREE and/or REFSEQ here if you don't have them.
d <- phyloseq(OTU, TAX, MAP, TREE, REFSEQ)
d

Think Before Attaching Reference Sequences to the Object

Loading the FASTA file into the phyloseq object works fine for a small project, but for a large one, thousands of zOTUs, long amplicons, or both, it can make the object noticeably heavier: slower to save and load, slower to subset and manipulate, and a much bigger .rds file to pass around or back up.

Pros of including it: sequences travel with the rest of your data as one self-contained object; convenient for building trees, running BLAST, or checking specific sequences directly from your R session without switching to another file.

Cons of including it: memory and disk overhead that scales with project size; every operation on the phyloseq object, even ones that have nothing to do with sequences, carries that extra weight along with it.

What to do instead, if it becomes a problem: skip refseq() when building the object, and keep the FASTA file as a separate, external file instead. taxa_names(ps) gives you the zOTU IDs, which you can use to look up specific sequences in the FASTA file only when you actually need them (e.g. with Biostrings::readDNAStringSet() on a subset, or command-line tools like samtools faidx / seqkit grep), rather than carrying every sequence around in memory for the whole analysis.

Why taxa_are_rows = TRUE?

phyloseq needs to know which orientation your count table uses, since it can't always tell rows from columns just by looking. Our count tables put OTUs/zOTUs as rows and samples as columns (see the OTU/Count Table page), which is why we set taxa_are_rows = TRUE. If you ever import a table the other way round and get a confusing dimension error, this is the first thing to check.

Handling Taxonomy

The tax_file step above depends on exactly how your .tax file lays out taxonomy: a single delimited SINTAX string in one column (most likely, something like d:Bacteria,p:Firmicutes,c:..., possibly with confidence values attached), or already-split columns per rank. Open the file and check before running this. If it's a single delimited string, you'll need an extra parsing step (e.g. tidyr::separate(), or a small custom function to strip the rank prefixes and any confidence values) to turn it into a proper rank-by-column matrix before tax_table() treats it the way you expect. Ask us if you're not sure what your specific file looks like.