`r knitr::opts_chunk$set(cache=FALSE, fig.width=12)`

# edgeR with phyloseq

Define the extension. Do not need to have anything loaded yet to do this.

```{r}
################################################################################
#' Convert phyloseq OTU count data into DGEList for edgeR package
#' 
#' Further details.
#' 
#' @param physeq (Required).  A \code{\link{phyloseq-class}} or
#'  an \code{\link{otu_table-class}} object. 
#'  The latter is only appropriate if \code{group} argument is also a 
#'  vector or factor with length equal to \code{nsamples(physeq)}.
#'  
#' @param group (Required). A character vector or factor giving the experimental
#'  group/condition for each sample/library. Alternatively, you may provide
#'  the name of a sample variable. This name should be among the output of
#'  \code{sample_variables(physeq)}, in which case
#'  \code{get_variable(physeq, group)} would return either a character vector or factor.
#'  This is passed on to \code{\link[edgeR]{DGEList}},
#'  and you may find further details or examples in its documentation.
#'  
#' @param method (Optional). The label of the edgeR-implemented normalization to use.
#'  See \code{\link[edgeR]{calcNormFactors}} for supported options and details. 
#'  The default option is \code{"RLE"}, which is a scaling factor method 
#'  proposed by Anders and Huber (2010).
#'  At time of writing, the \link[edgeR]{edgeR} package supported 
#'  the following options to the \code{method} argument:
#'  
#'  \code{c("TMM", "RLE", "upperquartile", "none")}.
#'
#' @param ... Additional arguments passed on to \code{\link[edgeR]{DGEList}}
#' 
#' @examples
#' 
phyloseq_to_edgeR = function(physeq, group, method="RLE", ...){
  require("edgeR")
  require("phyloseq")
  # Enforce orientation.
  if( !taxa_are_rows(physeq) ){ physeq <- t(physeq) }
  x = as(otu_table(physeq), "matrix")
  # Add one to protect against overflow, log(0) issues.
  x = x + 1
  # Check `group` argument
  if( identical(all.equal(length(group), 1), TRUE) & nsamples(physeq) > 1 ){
    # Assume that group was a sample variable name (must be categorical)
    group = get_variable(physeq, group)
  }
  # Define gene annotations (`genes`) as tax_table
  taxonomy = tax_table(physeq, errorIfNULL=FALSE)
  if( !is.null(taxonomy) ){
    taxonomy = data.frame(as(taxonomy, "matrix"))
  } 
  # Now turn into a DGEList
  y = DGEList(counts=x, group=group, genes=taxonomy, remove.zeros = TRUE, ...)
  # Calculate the normalization factors
  z = calcNormFactors(y, method=method)
  # Check for division by zero inside `calcNormFactors`
  if( !all(is.finite(z$samples$norm.factors)) ){
    stop("Something wrong with edgeR::calcNormFactors on this data,
         non-finite $norm.factors, consider changing `method` argument")
  }
  # Estimate dispersions
  return(estimateTagwiseDisp(estimateCommonDisp(z)))
}
################################################################################
```

# Citations

If you find this extension or tutorial useful in your work, please cite the following:

### Differential Abundance for Microbiome Data
McMurdie and Holmes (2014)
Waste Not, Want Not: Why Rarefying Microbiome Data is Inadmissible.
[PLoS Computational Biology 10(4):e1003531](http://dx.plos.org/10.1371/journal.pcbi.1003531)

### phyloseq
McMurdie and Holmes (2013)
phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data.
[PLoS ONE 8(4):e61217](http://dx.plos.org/10.1371/journal.pone.0061217)

### edgeR
Robinson MD, McCarthy DJ, Smyth GK (2009) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
[Bioinformatics (Oxford, England) 26: 139–140](http://bioinformatics.oxfordjournals.org/content/26/1/139.short)


# Load example data and try it out

```{r, message=FALSE, warning=FALSE}
date()
library("phyloseq"); packageVersion("phyloseq")
library("edgeR"); packageVersion("edgeR")
```


I use here the same example dataset as in the DESeq2 vignette, the publicly available data from a study on colorectal cancer:

[Genomic analysis identifies association of Fusobacterium with colorectal carcinoma](http://genome.cshlp.org/content/22/2/292.long).
Kostic, A. D., Gevers, D., Pedamallu, C. S., Michaud, M., Duke, F., Earl, A. M., et al. (2012). *Genome research*, 22(2), 292-298. 

This work was published ahead of print in [Genome Research](http://genome.cshlp.org/) alongside a highly-related article from a separate group of researchers (hooray for reproducible observations!): [Fusobacterium nucleatum infection is prevalent in human colorectal carcinoma](http://genome.cshlp.org/content/22/2/299.long). In case you are interested. For the purposes of example, however, we will stick to the data from the former study, with data available at the [microbio.me/qiime](http://www.microbio.me/qiime/) server.

Study ID:  `1457`

Project Name:  `Kostic_colorectal_cancer_fusobacterium`

Study Abstract:  The tumor microenvironment of colorectal carcinoma is a complex community of genomically altered cancer cells, nonneoplastic cells, and a diverse collection of microorganisms. Each of these components may contribute to carcino genesis; however, the role of the microbiota is the least well understood. We have characterized the composition of the microbiota in colorectal carcinoma using whole genome sequences from nine tumor/normal pairs. Fusobacterium sequences were enriched in carcinomas, confirmed by quantitative PCR and 16S rDNA sequence analysis of 95 carcinoma/normal DNA pairs, while the Bacteroidetes and Firmicutes phyla were depleted in tumors. Fusobacteria were also visualized within colorectal tumors using FISH. These findings reveal alterations in the colorectal cancer microbiota; however, the precise role of Fusobacteria in colorectal carcinoma pathogenesis requires further investigation.


```{r}
filepath = system.file("extdata", "study_1457_split_library_seqs_and_mapping.zip", package = "phyloseq")
kostic = microbio_me_qiime(filepath)
# Remove the samples for which the DIAGNOSIS was not included
kosticB = subset_samples(kostic, DIAGNOSIS != "None")
```

### Independent filtering

Independent filtering is done automatically in DESeq2, but not in DESeq. Here we will filter OTUs for which the variance across all samples is very low, and we'll do this before ever passing the data to DESeq. 

```{r var-plot-ind-filter}
kosticB
kosticp = transform_sample_counts(kosticB, function(x){x/sum(x)})
hist(log10(apply(otu_table(kosticp), 1, var)),
     xlab="log10(variance)", breaks=50,
     main="A large fraction of OTUs have very low variance")
varianceThreshold = 1e-5
keepOTUs = names(which(apply(otu_table(kosticp), 1, var) > varianceThreshold))
kosticB = prune_taxa(keepOTUs, kosticB)
kosticB
```

Here we've used an arbitrary but not-unreasonable variance threshold of `r varianceThreshold`. It is important to keep in mind that this filtering is independent of our downstream test. The sample classifications were not used.

Now let's use our newly-defined function to convert the phyloseq data object `kosticB` into an edgeR "DGE" data object, called `dge`.

```{r define-dge}
dge = phyloseq_to_edgeR(kosticB, group="DIAGNOSIS")
# Perform binary test
et = exactTest(dge)
# Extract values from test results
tt = topTags(et, n=nrow(dge$table), adjust.method="BH", sort.by="PValue")
res = tt@.Data[[1]]
alpha = 0.001
sigtab = res[(res$FDR < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(kosticB)[rownames(sigtab), ], "matrix"))
dim(sigtab)
head(sigtab)
```

Here is a bar plot showing the log2-fold-change, showing Genus and Phylum. Uses some ggplot2 commands.

```{r, fig.width=11}
library("ggplot2"); packageVersion("ggplot2")
theme_set(theme_bw())
scale_fill_discrete <- function(palname = "Set1", ...) {
    scale_fill_brewer(palette = palname, ...)
}
sigtabgen = subset(sigtab, !is.na(Genus))
# Phylum order
x = tapply(sigtabgen$logFC, sigtabgen$Phylum, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Phylum = factor(as.character(sigtabgen$Phylum), levels = names(x))
# Genus order
x = tapply(sigtabgen$logFC, sigtabgen$Genus, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Genus = factor(as.character(sigtabgen$Genus), levels = names(x))
ggplot(sigtabgen, aes(x = Genus, y = logFC, color = Phylum)) + geom_point(size=6) + 
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5))
```

As expected from the original study abstract and title, *Fusobacterium* OTUs were most-significantly differentially abundant between the cancerous and healthy samples. If you look closely, two different genera of the Fusobacteria phylum were among the most significantly different, *Leptotrichia* (the winner) as well as *Fusobacterium*.

---

## Paired tests

As mentioned above, the design of this experiment is 95 carcinoma/normal pairs, where each pair comes from the same patient. Although the previous tests are valid, they are conservative in that they do not use this extra information regarding the sample-pairs, and in that sense have forfeited extra power. There is support in edgeR for paired tests, and this is officially described in [one of the edgeR user guides](http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf). It is also demonstrated here in the following.

```{r paired-tests}
Diagnosis = get_variable(kosticB, "DIAGNOSIS")
Patient = get_variable(kosticB, "ANONYMIZED_NAME")
# Notice that we have some patients without one of the pairs.
tapply(Patient, Diagnosis, length)
length(levels(Patient))
any(tapply(Diagnosis, Patient, length) > 2)
sum(tapply(Diagnosis, Patient, length) < 2)
# Keep only patients with both healthy and cancer samples
keepPatients = names(which(tapply(Diagnosis, Patient, function(x){length(unique(x))}) == 2))
kosticBpair = subset_samples(kosticB, ANONYMIZED_NAME %in% keepPatients)
Diagnosis = get_variable(kosticBpair, "DIAGNOSIS")
Patient = get_variable(kosticBpair, "ANONYMIZED_NAME")
# With that accounting out of the way, define the design matrix
design = model.matrix(~ Patient + Diagnosis)
```

Must estimate the dispersions, as always. This is different than in the function shown above.

```{r}
# Add one to protect against overflow, log(0) issues.
x = as(otu_table(kosticBpair), "matrix") + 1L
taxonomy = data.frame(as(tax_table(kosticBpair), "matrix"))
# Now turn into a DGEList
x = DGEList(counts=x, group=Diagnosis, genes=taxonomy, remove.zeros=TRUE)
# Calculate the normalization factors and estimate dispersion
x = calcNormFactors(x, method="RLE")
x = estimateGLMCommonDisp(x, design)
x = estimateGLMTrendedDisp(x, design)
x = estimateGLMTagwiseDisp(x, design)
```

As in the edgeR User's Guide, we proceed to fit a linear model and test for the treatment effect. Note that we can omit the coefficient argument to `glmLRT` because the "treatment effect" (in this case the tissue diagnosis) is the last coeffcient in the model.

```{r paired-test-glmfit}
fit <- glmFit(x, design)
lrt <- glmLRT(fit)
topTags(lrt)
```
This test detects OTUs that are differentially abundant in the tumor colon mucosa relative to the healthy colon mucosa (control), adjusting for baseline differences between the patients. This test can be viewed as a generalization of a paired t-test.

Re-make the plot.

```{r paired-ggplot}
respair = topTags(lrt, n=nrow(x), adjust.method="BH", sort.by="PValue")
respair = respair@.Data[[1]]
alpha = 0.001
sigtab = respair[(respair$FDR < alpha), ]
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(kosticB)[rownames(sigtab), ], "matrix"))
dim(sigtab)
head(sigtab)
sigtabgen = subset(sigtab, !is.na(Genus))
# Phylum order
x = tapply(sigtabgen$logFC, sigtabgen$Phylum, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Phylum = factor(as.character(sigtabgen$Phylum), levels = names(x))
# Genus order
x = tapply(sigtabgen$logFC, sigtabgen$Genus, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Genus = factor(as.character(sigtabgen$Genus), levels = names(x))
ggplot(sigtabgen, aes(x = Genus, y = logFC, color = Phylum)) + geom_point(size=6) + 
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.5)) +
  ggtitle("Log Fold Change of Significant OTUs in a Paired Test")
```

### Other covariates available

As a side note, there are many other interesting patient-sample covariates available in the `sample_data`, which you can access with `sample_data()` and `get_variable()`.

```{r show-variables}
sample_variables(kostic)
```



---

## Other extensions for the phyloseq package:

```{r other-graphics-links, results='asis', echo=FALSE}
htmlfiles = list.files(pattern="\\.html")
htmlstems = gsub(".html", "", htmlfiles)
mdlines = paste("#### [", htmlstems, "](", htmlfiles, ")", sep="")
cat(mdlines, sep="\n\n")
```
