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

# Human Microbiome Project

Example importing into [phyloseq](http://joey711.github.com/phyloseq/) the files produced by [Qiime](http://qiime.org/) run on [the HMPv35 dataset](http://hmpdacc.org/micro_analysis/microbiome_analyses.php), which is avilable from [HMP-DACC](http://hmpdacc.org/). Even more specifically this data is from [the HMP QIIME Community Profiling data page](http://hmpdacc.org/HMQCP/). This is important to note because there are also "parallel versions" of this data [generated by mothur](http://hmpdacc.org/HMMCP/) which could also be imported by phyloseq; and there are also a growing number of data links and displays added to the main HMP-DACC page that might otherwise obfuscate where this data came from. However, the precise links are shown in the code below. In this way we are demonstrating both the `import_qiime` function in phyloseq, as well as the import of web-facing data in general using R.

Note that this is a document produced by [the markdown package](http://cran.r-project.org/web/packages/markdown/index.html) from an [R Markdown file](http://rstudio.org/docs/authoring/using_markdown). Markdown is a simple formatting syntax for authoring web pages, and in this case the structure includes reproducible code.

## Load [the phyloseq package](http://joey711.github.com/phyloseq/)
```{r load_phyloseq, message=FALSE, warning=FALSE}
library("phyloseq"); packageVersion("phyloseq")
```
Because this demo includes an example of manually importing a Fasta-formatted DNA sequence file using a function from [the Biostrings package](), we will need to load that as well (it is already a phyloseq dependency, so nothing exotic here)
```{r load_Biostrings, message=FALSE, warning=FALSE}
library("Biostrings"); packageVersion("Biostrings")
```

As data sources could be updated (or their links removed), and the phyloseq package is itself under active development, I am posting here the date/time stamp that this demo was built, as extra compatibility metadata to supplement the package version numbers shown above

```{r timestamp}
date()
```


## Download the files to a temporary location 
Note that R's file handling is sophisticated enough that you do not need to unpack/unzip the `.gz`- or `.bz2`-compressed files. They are understood automatically.

### Download the OTU abundance file
Create a temporary local file where there large otu-tax file will be downloaded and read-from.

```{r download-otutax}
temp_otutax = tempfile()
otu_url     = "http://downloads.hmpdacc.org/data/HMQCP/otu_table_psn_v35.txt.gz"
download.file(otu_url, temp_otutax)
```
### Download the Sample Map File

```{r download-sample-map}
temp_map = tempfile()
map_url  = "http://downloads.hmpdacc.org/data/HMQCP/v35_map_uniquebyPSN.txt.bz2"
download.file(map_url, temp_map)
temp_map = bzfile(temp_map)
```

### Download the phylogenetic tree
```{r download-tree}
temp_tree = tempfile()
tree_url = "http://downloads.hmpdacc.org/data/HMQCP/rep_set_v35.tre.gz"
download.file(tree_url, temp_tree)
```

It turns out that the tree has weird quotes added around tip names, for example. The following will read the `.gz` tree directly (no extra uncompression step needed).
```{r import-tree-fix}
tree = read_tree(temp_tree) 
head(taxa_names(tree))
```

These extra quotes ( `"'"` ) must be removed/replaced. This can be accomplished using the following `gsub` command. Even for `r ntaxa(tree)` OTUs this is very fast.
```{r}
tree$tip.label <- gsub("\'", "", tree$tip.label, fixed=TRUE)
head(tree$tip.label)
```

### Download the representative sequences
Unfortunately, the `readDNAStringSet` function for importing the representative sequences -- which we will call `refseq` (for "reference sequences") -- does not accept standard R file connections for reading files, possibly because it is performing the parsing at a very low level? In any case, unlike the other files we are importing, this must be unzipped first to an uncompressed `.fna` file, before we read and parse.
```{r download-refseq}
zipped_refseq   = tempfile()
refseq_url = "http://downloads.hmpdacc.org/data/HMQCP/rep_set_v35.fna.gz"
download.file(refseq_url, zipped_refseq)
```

I use the `gzfile` function to open an R connection to the `.gz` compressed representative sequence file that we just downloaded. Then, use `readLines` and `writeLines` to essentially "pipe" this uncompressed to our second temporary file, called `unzipped_rsfile`.
```{r decompress-refseq}
unzipped_rsfile = tempfile()
refseqgzcon = gzfile(zipped_refseq)
writeLines(readLines(refseqgzcon), unzipped_rsfile)
```

Now that I have an uncompressed version of the file saved in as temporary (but known) file, I can use the `readDNAStringSet` function without error and store the initial import result as the object, `refseq`.
```{r decompress}
refseqs = readDNAStringSet(unzipped_rsfile)
refseqs
head(names(refseqs))
```

We need to parse the "names" of each sequence entry in `refseq` such that they match the format of the OTU names in the other files. In this case, that means removing everything from the first space on. There are multiple ways to do this using R, and in the following code we use regular-expression pattern matching and substitution available in the `gsub` function.
```{r}
names(refseqs) = gsub("(^[[:print:]]{7,13})([[:space:]]{1,})([[:print:]]{0,}$)", "\\1", names(refseqs))
head(names(refseqs))
```


## Import the Data and Create phyloseq Object.
I have included a timing function call so that you can see how long this large file took to import using my laptop. Also, the number of "chunks" is shown by the number of "dots" printed after the function call. Note that the phylogenetic tree was already imported, slightly modified in the previous code chunk, and we want this to provide this `tree` as the tree argument, rather than the raw tree file, which won't exactly match the `taxa_names` in the OTU-Tax file.

We have actually already imported the tree and reference sequences as `tree` and `refseqs`, respectively; and we are merely passing them along to be included in the object that is output by `import_qiime`. The hard work and the bulk of the processing time for `import_qiime` will be parsing the OTU-abundance/taxonomy text file in RAM-manageable chunks, prior to combining with the other data components and storing as the phyloseq object, `HMPv35`.

Note that you will multiple instances of the following warning:
```{r}
## Warning: No greengenes prefixes were found.  Consider using
## parse_taxonomy_default() instead if true for all OTUs.  Dummy ranks may be
## included among taxonomic ranks now.
```
These arise from OTUs that have only been assigned a "root" taxonomic classification -- which is meaningless -- and no other. In this dataset, the "root" has no greengenes prefix but all the other taxonomic classifications do. When the default greengenes parsing function (`parse_taxonomy_default`) finds these OTUs, it complains that there was not actually a greengenes prefix in any of the taxonomy elements. That's okay, your data (and this data) is fine. It is just a warning, and all the meaningless "root" classification elements are assigned to an equally meaningless dummy taxonomic rank called "Rank1". See for yourself below in the output of `rank_names`.

```{r import-HMPv35, warning=FALSE}
system.time(HMPv35 <- import_qiime(temp_otutax, temp_map, tree, refseqs))
```

## Investigate Basic Features of HMPv35
For more advanced tools for exploring the data, see [additional phyloseq documentation](http://joey711.github.com/phyloseq/)

```{r explore-HMPv35}
ntaxa(HMPv35)
nsamples(HMPv35)
sample_variables(HMPv35)
rank_names(HMPv35)
head(tax_table(HMPv35))
```

## Save HMPv35 to an RData File
Save `HMPv35` to an `".RData"`" file so that you don't have to run this import ever again. In the following example, I'm saving it to the current directory of my process. Your preferred file path for saving R data may (and probably should) vary.

```{r}
save(HMPv35, file="HMPv35.RData", compress="bzip2")
```

## Alternatively, Download the Processed File from Us
We re-run this import code each time we update/rebuild these demo vignettes. This [already-imported `HMPv35` dataset](HMPv35.RData) is available by cloning or downloading the phyloseq-demo repository. Typically this `RData` binary file will load into R much, much faster than importing it again from the public HMP data files.

