`r opts_chunk$set(cache=TRUE, fig.width=10)`
<link href="http://joey711.github.com/phyloseq/markdown.css" rel="stylesheet"></link>


DIY: public restroom bacteria
========================================================
[Originally hosted here](https://github.com/joey711/phyloseq-demo/Restroom-Biogeography)

### Goals of this demonstration
- Highlight an interesting use of R (analyzing bacteria in restrooms)
- Demonstrate for R users how to access publicly-available microbiome census data. In this case counts of bacteria from a survey of public restroom sufaces in 2011.
- Scholarly reproduction / reanalysis of the findings/techniques in the original article, noting places where reproduction is difficult or impossible and where there may be weaknesses in the statistical approach.
- Demonstrate the use of [phyloseq](http://joey711.github.io/phyloseq/) and R for conducting reproducible statistical analysis of microbiome census data.
- Demonstrate some of the powerful features of [the ggplot2 package](http://cran.r-project.org/web/packages/ggplot2/index.html) in creating graphics about this data.

The **[original article by Flores et al.](http://dx.plos.org/10.1371/journal.pone.0028132)** was published in 2011 under an open access license by the journal [PLoS ONE](http://www.plosone.org/). As of the time of this writing, it has been viewed almost 22,000 times.
Flores GE, Bates ST, Knights D, Lauber CL, Stombaugh J, et al. (2011) **Microbial Biogeography of Public Restroom Surfaces**. PLoS ONE 6(11): e28132.

**The data** is hosted at [microbio.me/qiime](http://www.microbio.me/qiime/), with Study ID **1335**, Project Name `Flores_restroom_surface_biogeography`. The zipfile of the data is exposed at [this FTP address](ftp://thebeast.colorado.edu/pub/QIIME_DB_Public_Studies/study_1335_split_library_seqs_and_mapping.zip) -- note that you may have to create an account and login to access the `.zip` file. Otherwise I would include downloading and unzipping as part of the R source code in this demo.

The time that this demo was built: `r date()`

## Study design and findings

Here is an illustration of the study design (their Figure 3), namely, **where in public restrooms they sampled**; as well as a few overall observations that they included in their caption.

<img src="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0028132.g003/largerimage" width="750px" />

###### **Cartoon illustrations of the relative abundance of discriminating taxa on public restroom surfaces.** Light blue indicates low abundance while dark blue indicates high abundance of taxa. (Panel A) Although skin-associated taxa (Propionibacteriaceae, Corynebacteriaceae, Staphylococcaceae and Streptococcaceae) were abundant on all surfaces, they were relatively more abundant on surfaces routinely touched with hands. (Panel B) Gut-associated taxa (Clostridiales, Clostridiales group XI, Ruminococcaceae, Lachnospiraceae, Prevotellaceae and Bacteroidaceae) were most abundant on toilet surfaces. (Panel C) Although soil-associated taxa (Rhodobacteraceae, Rhizobiales, Microbacteriaceae and Nocardioidaceae) were in low abundance on all restroom surfaces, they were relatively more abundant on the floor of the restrooms we surveyed. Figure not drawn to scale.


You can also check out [this image at the original article](http://www.ncbi.nlm.nih.gov/core/lw/2.0/html/tileshop_pmc/tileshop_pmc_inline.html?title=Click%20on%20image%20to%20zoom&p=PMC3&id=3223236_pone.0028132.g003.jpg).


## Load phyloseq and other packages

```{r load-packages, message=FALSE, warning=FALSE}
library("phyloseq"); packageVersion("phyloseq")
```

I also want to make some nice graphics using the ggplot2 package, so I will also load that and adjust its default theme

```{r ggplot-load-theme}
library("ggplot2"); packageVersion("ggplot2")
theme_set(theme_bw())
```


## Import data
The microbio.me/qiime site has provided a zip file with two key data files. Namely, a file with the abundance data and some metadata, and a tab delimited text file with sample meta data. This following code chunk imports both files, and then combines them in one integrated phyloseq object that we will use in the remaining (re)analysis.

First you need to download the file directly from the link provided above, or use the zipfile we have mirrored in our phyloseq-demo GitHub repository.

Here is how you would download it from the FTP address to a random temporary file defined by your operating system. We will call this file `zipfile`.
```{r download-zipfile, eval=FALSE}
zipftp = "ftp://thebeast.colorado.edu/pub/QIIME_DB_Public_Studies/study_1335_split_library_seqs_and_mapping.zip"
zipfile = tempfile("RestroomBiogeography")
download.file(zipftp, zipfile)
```

Since I have reposted the zipped data file on the phyloseq-demo repository on GitHub -- in the same directory as the original source file for this post -- the `zipfile` in this example can also be defined locally, as shown here:
```{r zipfile-ghpages}
zipfile = "study_1335_split_library_seqs_and_mapping.zip"
```

In either case, it is now possible with R to unzip to a randomly labeled temporary diretory created specific to your operating system, which we will name with the symbol `import_dir`. My understanding is that eventually temporary files and directories are deleted at system shutdown or startup. Probably system dependent exactly when they get deleted.

Here is how to create the temporary directory for the unpacked file(s) from the `.zip` file.
```{r unzip-file}
import_dir <- tempdir()
unzip(zipfile, exdir=import_dir)
```

Import from the `.biom` file.
```{r import-data-biom}
biomfile = paste0(import_dir, "/study_1335_closed_reference_otu_table.biom")
biom = import_biom(biomfile, parseFunction=parse_taxonomy_greengenes)
```

Import from the sample data file.
```{r import-data-sampledata}
sdfile = paste0(import_dir, "/study_1335_mapping_file.txt")
sample_metadata = import_qiime_sample_data(sdfile)
```

Combine (merge) the two data objects.
```{r import-data-merge}
restroom = merge_phyloseq(biom, sample_metadata)
```

Now we have available a new combined data object, called `restroom`, that contains all the data we should need for this tutorial. We can query different features of `restroom` using the phyloseq API in the form of accessor functions/methods. We can also just print it to standard out and we get some summary information.

```{r accessor-examples}
restroom
rank_names(restroom)
sample_variables(restroom)
levels(sample_data(restroom)$BUILDING)
levels(sample_data(restroom)$SURFACE)
```

Okay, let's move on to preprocessing the data as needed, and reproducing the figures in the article.

## Preprocessing
The phyloseq package includes a [tutorial on preprocessing microiome census data](http://joey711.github.io/phyloseq/preprocess.html) with several broad and flexible methods. Here we want to use one of those preprocessing approaches as well as the one used by Flores et al. in this article. Unfortunately, they used random subsampling of their data with on reported seed or random number generation method, so it is impossible to exactly recapitulate their randomly subsampled data. However, we can do our best using a function included in phyloseq for exactly this purpose, `rarefy_even_depth`. 

Are there any OTUs included in this dataset that have no counted reads at all prior to preprocessing? If so, how many?

```{r remove-zeros}
any(taxa_sums(restroom) == 0)
sum(taxa_sums(restroom) == 0)
```
So of the `r ntaxa(restroom)` OTUs in the dataset, there were `r sum(taxa_sums(restroom) == 0)` OTUs that were apparently not observed even once in any of the samples (`r 100*round(sum(taxa_sums(restroom) == 0)/ntaxa(restroom), 1)`%). This sounds very odd, but more likely this is evidence of some data "massaging" that removed the abundance values of those taxa, but their entries in the `.biom` file are inexplicably included.

Here I remove those "unobserved" OTUs. Note how I first save the originally-imported data as `restroom0` in case I want to go back to it after modifying/preprocessing `restroom`.

```{r remove-zero-OTUs}
restroom0 = restroom
restroom = prune_taxa(taxa_sums(restroom) > 0, restroom)
```

Any empty samples? Apparently not, though the code for pruning "empty" samples is also shown. And procedes quickly since there is nothing in `restroom` to modify.

```{r remove-zero-samples}
any(sample_sums(restroom) == 0)
restroom = prune_samples(sample_sums(restroom) > 0, restroom)
```

What about the total reads per sample, and what does the distribution look like?
```{r, fig.height=2.5, warning=FALSE}
readsumsdf = data.frame(
  nreads = sort(taxa_sums(restroom), TRUE),
  sorted = 1:ntaxa(restroom),
  type="OTUs"
)
readsumsdf = rbind(readsumsdf, data.frame(
    nreads = sort(sample_sums(restroom), TRUE),
    sorted = 1:nsamples(restroom),
    type="Samples"
  )
)
title = "Total number of reads"
p = ggplot(readsumsdf, aes(x=sorted, y=nreads)) + geom_bar(stat="identity")
p + ggtitle(title) + scale_y_log10() + facet_wrap(~type, 1, scales="free")
```

This looks pretty typical for the distribution of reads from an amplicon-based microbiome census, if not even surprisingly evenly distributed across most samples... I've seen much, much worse. In any case, one of the main preprocessing steps we will do is to transform the abundances of each sample in a way that will make the represented "sequencing effort", AKA the sum of the counts, the same across all samples. As mentioned earlier, I will redo a random subsampling ("rarefying") to 500 reads as the original authors did; and, separately, I will also simply transform the abundance vectors of each sample by multiplying the OTU proportion within each sample by the arbitrary value, `500` used in the random subsampling. We can later compare the how the two datasets perform.

First, I will remove the samples with a total number of reads below the threshold set by the original uauthors of `500`. 

```{r remove-500}
restroom = prune_samples(sample_sums(restroom) > 500, restroom)
```

This resulted in the removal of `r nsamples(restroom0) - nsamples(restroom)` sample(s) from the dataset.

### Transform to equal sample sum

**Rarefy**
We will call this `restroomR`. The way to accomplish this in phyloseq is to **first define a random number seed** so that others can reproduce your subsetted data, and then use the `rarefy_even_depth` function. Here I use the "ePub" number assigned to this article by PLoS ONE as the seed. Seemed appropriate, but totally arbitrary.

```{r rarefy}
set.seed(28132)
restroomR = rarefy_even_depth(restroom, sample.size=500)
```

And now an alternative to random subsampling, a simple proportional transformation, calling it `restroomP`.

```{r transform-proportional}
restroomP = transform_sample_counts(restroom, function(x) 500*x/sum(x) )
```

For sanity-check, let's replot the sample sums of each of these new data objects, to convince ourselves that all of the samples now sum to 500.
```{r sample-sums-transformed, fig.height=3}
par(mfrow=c(1, 2))
title = "Sum of reads for each sample, restroomR"
plot(sort(sample_sums(restroomR), TRUE), type="h", main=title, ylab="reads",
  ylim=c(0, 1000))
title = "Sum of reads for each sample, restroomP"
plot(sort(sample_sums(restroomP), TRUE), type="h", main=title, ylab="reads", 
  ylim=c(0, 1000))
```


---

## Figure 1
### Original Figure

<img src="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0028132.g001/largerimage" width="750px" />

###### Original Figure caption, quoted here: "Figure 1. Taxonomic composition of bacterial communities associated with public restroom surfaces.(A) Average composition of bacterial communities associated with restroom surfaces and potential source environments. (B) Taxonomic differences were observed between some surfaces in male and female restrooms. Only the 19 most abundant taxa are shown. For a more detailed taxonomic breakdown by gender including some of the variation see [Supplemental Table S2](http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0028132#pone.0028132.s002)."

### Figure 1 - Reproduced
I will use the phyloseq plotting function `plot_bar` to reproduce Figure 1. I have taken liberties to improve the representation of the data to make certain features more comparable.

First, they have mixed taxonomic ranks in this plot between family and order, but most are family so I will use that. Furthermore, they showed only the most-abundant 19 OTUs, so we will add a variant taxonomic rank to the data, called `Family19`, that will only have non-`r NA` values if the OTU is among those top 19.
```{r}
top19otus = names(sort(taxa_sums(restroomR), TRUE)[1:19])
taxtab19  = cbind(tax_table(restroomR), family19=NA)
taxtab19[top19otus, "family19"] <- as(tax_table(restroomR)[top19otus, "Family"], "character")
tax_table(restroomR) <- tax_table(taxtab19)
```

Now create the abundance bar plot of `restroomR` using the `plot_bar` function. 
```{r Figure1-redo-01-A}
title = "Figure 1 Part A (remake), attempt 1"
plot_bar(restroomR, "SURFACE", fill="family19", title=title) + coord_flip()
```

Hey that was easy! But the abundance values get over 500. How come? Oh, each of these categories includes more than one sample, so we expect them to exceed 500. In fact, you can see that the totals are multiples of 500. This is easy to fix by re-doing our earlier transformation on `restroomR` after mergeing by category. The phyloseq package includes a function specifically for doing this kind of categorical merge, called `merge_samples`.

```{r surface-merge-figure-01}
restroomRm = merge_samples(restroomR, "SURFACE")
```
Repair the merged values associated with each surface after merge.
```{r surface-merge-repair}
sample_data(restroomRm)$SURFACE <- levels(sample_data(restroomR)$SURFACE)
```

Transform to percentages of total available.
```{r transform-surface-merge}
restroomRm = transform_sample_counts(restroomRm, function(x) 100*x/sum(x))
```

Now we can plot Figure 1A with each category on equivalent footing. 
```{r Figure1-redo-01-A-2}
title = "Figure 1 Part A (remake), attempt 2"
plot_bar(restroomRm, "SURFACE", fill="family19", title=title) + coord_flip() + ylab("Percentage of Sequences")
```

To better match Figure 1 Panel A from the original article, I can remove the gray rectangles that represent OTUs that were not among the most abundant 19. 

```{r prune-non-19}
restroomRm19 = prune_taxa(top19otus, restroomRm)
title = "Figure 1 Part A (remake), attempt 3"
plot_bar(restroomRm19, "SURFACE", fill="family19", title=title) + coord_flip() + ylab("Percentage of Sequences") + ylim(0, 100)
```
###### Re-plot of Figure 1A from original article using phyloseq. Each rectangle is a different OTU. Since these are multi-sample categories, rather than simply samples, the length of each bar represents the mean fraction abundance of that OTU among the normalized samples in that category. Multiple rectangles with the same color and category represent relative abundances of different OTUs in the same taxonomic family. Gray boxes are OTUs that were among the top 19 most abundant, but did not have a taxonomic family classification in the taxonomy table.


---

### Figure 1, Panel B -- Gender
What about separating by gender? 

First, we will merge samples by a dummy variable representing both `GENDER` and `SURFACES`, repair flattened factors, and transform to proportional abundances.
```{r gender-surface-merge}
sample_data(restroomR)$gensurf <- paste0(sample_data(restroomR)$GENDER, sample_data(restroomR)$SURFACE) 
restroomRgsm = merge_samples(restroomR, "gensurf")
# repair factors
sample_data(restroomRgsm)$SURFACE <- levels(sample_data(restroomR)$SURFACE)[get_variable(restroomRgsm, "SURFACE")]
sample_data(restroomRgsm)$GENDER <- levels(sample_data(restroomR)$GENDER)[get_variable(restroomRgsm, "GENDER")]
# transform to proportions
restroomRgsmp = transform_sample_counts(restroomRgsm, function(x) 100*x/sum(x))
```

Now plot with panelling by surface for an easier gender-wise comparison than in the original figure.

```{r Figure-01-partB-v01}
restroomRgsm19 = prune_taxa(top19otus, restroomRgsmp)
title = "Figure 1 Part B, Restroom Surface and Gender"
p = plot_bar(restroomRgsm19, "GENDER", fill="family19", title=title) + coord_flip() + labs(colour="family")
p + facet_wrap(~SURFACE)
```




---

## Figure 2 - Dimensional Reduction
The original authors use [Multidimensional Scaling]() (also called Principle Coordinates Analysis) to decompose a pairwise distance matrix between all the microbial samples. The original distances are calculated by what microbial ecologists call the ["unweighted UniFrac" distance](). This approach is useful for graphically displaying high-level trends across all the samples in the dataset, usually by overlaying additional information about each sample, for instance, where it came from.

Unfortunately, calculating the UniFrac distance between each sample requires having an evolutionary (phylogenetic) tree of the bacterial species in the dataset, and this wasn't provided. We will instead attempt to create an equally-illuminating scatterplot using a different distance matrix.

### Figure 2 - Original Article

<img src="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0028132.g002/largerimage" width="700px" />
**Figure 2. Relationship between bacterial communities associated with ten public restroom surfaces.**
Communities were clustered using PCoA of the unweighted UniFrac distance matrix. Each point represents a single sample. Note that the floor (triangles) and toilet (asterisks) surfaces form clusters distinct from surfaces touched with hands.

### Figure 2 - Reproduced

The multidimensional scaling of a UniFrac distance matrix is one instance of a larger category of what we call *ordination methods*. These methods attempt to decompose the variability of higher-dimensional data into a smaller number of orthogonal (perpendicular) axes that turn out to be pretty useful for plotting and certain clustering methods.

```{r phyloseq-ordinate-figure02}
title = "Figure 2 remake, PCoA of Bray-Curtis distance"
restroom_ord = ordinate(restroomR, "PCoA", "bray")
p = plot_ordination(restroomR, restroom_ord, color="SURFACE", shape="GENDER")
p = p + geom_point(size=6, alpha=0.7) + ggtitle(title)
p
p + facet_wrap(~SURFACE)
p + facet_wrap(~GENDER)#, nrow=3)
```


# Table 1
<img src="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0028132.t001/largerimage" width="700px" />
###### Table 1 original article. ANOSIM test. Results of pairwise comparisons for unweighted UniFrac distances of bacterial communities associated with various surfaces of public restrooms on the University of Colorado campus using the ANOSIM test in Primer v6.

In addition to `Primer v6`, the ANOSIM method is implemented in [the vegan package](http://cran.r-project.org/web/packages/vegan/index.html) in R, and described nicely in [the anosim documentation](http://www.inside-r.org/packages/cran/vegan/docs/anosim). It is designed to operate directly on a dissimilarity (distance) matrix, so it is unclear why the authors first decomposed their unweighted UniFrac distance matrix with Principle Coordinate Analysis, and then provided those coordinates (presumably from just the first two axes?) to the ANOSIM implementation in *Primer v6*. 

**Important note about multiple testing**
The `anosim` function performs a non-parametric test of the significance of the sample-grouping you provide against a permutation-based null distribution, generated by randomly permuting the sample labels many times (999 permutations is the default, used here). The table shown in the original article appears to be the result of the `r choose(10, 2)` two-surface (pairwise) `anosim` tests, excluding water. The precise details are not explained, and we cannot repeat them here because **(1)** we cannot exactly reproduce the randomly subsampled ("rarefied") abundance data from the somewhat raw abundance matrix that is hosted on `microbio.me`, and more importantly **(2)** the dataset does not include the phylogenetic tree that is necessary for calculating the unweighted UniFrac distance matrix. Critically, the authors made **no mention of any correction for having conducted `r choose(10, 2)` simultaneous tests**. They did take a slightly stricter-than-usual threshold for significance of `P <= 0.01`, but we don't know without a formal correction if the tests with `P` close to `0.01` are actually significant (even at the more common `0.05` level), and the P-values themselves were not reported, just the statistic with and without asterisk. 

Here is how to re-perform `anosim` using R and the vegan package:

```{r vegan}
library("vegan"); packageVersion("vegan")
```

First look at results for surface types.
```{r anosim-surface}
surface_group = get_variable(restroomR, "SURFACE")
surface_ano = anosim(distance(restroomR, "bray"), surface_group)
surface_ano$signif
surface_ano$statistic
```

Now test bathroom genders. First remove samples for which there was no gender to make this a pairwise test.
```{r anosim-gender}
restroomRg = subset_samples(restroomR, !GENDER == "None")
gender_group = get_variable(restroomRg, "GENDER")
gender_ano = anosim(distance(restroomRg, "bray"), gender_group)
gender_ano$signif
gender_ano$statistic
```

Here we settled for a Bray-Curtis distance matrix for our own version of just two `anosim` calculations, one of which would not be considered significant by most standards (gender). 


### Add a demonstration of a multiple-hypothesis correction?

I considered adding a demonstration of how to add a multiple correction for the approach the authors took, but I'm not convinced (yet) that this was the optimal statistical approach to arrive at the knowledge that they were seeking. For one thing, it quickly loses power as the number of categories increases because of the larger number of corrections, not to mention it becomes difficult to calculate. Practicality aside, even the documentation for `anosim` in R suggests users use `adonis` (a permutation MANOVA, also in the vegan package) as "a more robust alternative".

```{r adonis}
df = as(sample_data(restroomR), "data.frame")
d = distance(restroomR, "bray")
restroomadonis = adonis(d ~ SURFACE + GENDER + BUILDING + FLOOR, df)
restroomadonis
plot(restroomadonis$aov.tab)
```


---

# Figure 4 -- Microbial Source Tracking

## Figure 4 -- Original Figure
<img src="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0028132.g004/largerimage" width="700px" />


## Soil source

The soil source data was taken from [this 2009 article by Lauber et al. in AEM](http://aem.asm.org/content/75/15/5111.full)

## SourceTracker

The authors describe using an algorithm called "SourceTracker", previously described in [Knights et al. (2011) Nature Methods](http://www.nature.com/nmeth/journal/v8/n9/full/nmeth.1650.html), and provided as a [commented R script published at sourceforge/sourcetracker](http://sourceforge.net/projects/sourcetracker/).

In perhaps another post, or a continuation of this one, I will include their code and apply it to the available data here. One additional challenge is the availability of the refernce data, which was not repackaged or republished with this dataset. Hopefully it is all still available.

# Feedback welcome!

Please feel free to post comments, suggestions. Also, ideas for other demos are strongly encouraged. A good place for that is probably [the issue tracker for the phyloseq-demo repository](https://github.com/joey711/phyloseq-demo/issues) where this document is versioned and where many other phyloseq-related demos are hosted.

# Thanks to original authors
Thanks to original authors for an interesting study and for making (most of) their data so publicly available. I just subjected it to some scholarly scrutiny, which is a good thing (and I hope they agree); and only possible because of what they made available. If you like the idea of being able to re-analyze this and other studies, then clamor for **better reproducible research**, which should hopefully the original *and* intermediate data *as well as* analysis source code ... and strive to do the same in your own publications. We will all do better, faster work because of it :)

