ranacapa: Explore output from Anacapa in R

How to install:

0. Install the phyloseq and devtools packages into R if you don't already have them:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("multtest")
BiocManager::install("phyloseq")

if (!requireNamespace("devtools", quietly = TRUE))
      install.packages('devtools')

1. Install ranacapa:

library(devtools)

devtools::install_github("gauravsk/ranacapa")

2. The package installation may take a few minutes.

How to run the interactive shiny app:

1. Launch the shiny app by running the following line in R:

library(ranacapa)
ranacapa::runRanacapaApp()

How to run, option 2 (Only relevant for command line Anacapa users. Scripted output of all possible graphs and tables, not interactive):

1. Download the ranacapa_automated.R file from the Anacapa repo: https://github.com/limey-bean/Anacapa/blob/master/Anacapa_db/scripts/downstream-analyses/ranacapa_automated.R

2. In the terminal the script runs as:

Rscript ranacapa_automated.R /path/to/input_biom.txt /path/to/input_metadata.txt /path/to/desired_output_directory rarefaction_depth rarefaction_replicates

---

Structure of the package

Last updated 22 April 2018

Input files

1. As of right now the package is designed to visualize the output from the Anacapa pipeline. At minimum, the Anacapa output file has one column that contains the full taxonomic path assigned to the read, and one column for each of the samples sequenced. A file may look something like this:

site_1 site_2 site_3 site_4 site_5 site_6 site_7 site_8 site_9 sum.taxonomy                                                         
      1      0      0      0      0      0      0      0      0 Annelida;Clitellata;Haplotaxida;Megascolecidae;Amynthas;Amynthas sze…
      0      0      0      0      0      0      0      0      0 Nemertea;Palaeonemertea;NA;Cephalothricidae;Cephalothrix;Cephalothri…
      0      0      0      0      0      0      0      0      0 ""                                                                   
      0      0      0      0      0      0      0      0      0 Nematoda;Chromadorea;Monhysterida;Monhysteridae;NA;Monhysteridae sp.…
      0      1      0      0      0      0      0      0      0 NA;Oomycetes;Pythiales;Pythiaceae;Pythium;Pythium rostratifingens    

In some cases, the output file also contains a column called [xxx]_seq_number that contains a unique sequence identifier for each row. Here's what the file would look like if the sequence number was written by Anacapa:

CO1_seq_number    site_1 site_2 site_3 site_4 site_5 site_6 site_7 site_8 site_9 sum.taxonomy                                       
forward_CO1_9136       1      0      0      0      0      0      0      0      0 Annelida;Clitellata;Haplotaxida;Megascolecidae;Amy…
forward_CO1_13513      0      0      0      0      0      0      0      0      0 Nemertea;Palaeonemertea;NA;Cephalothricidae;Cephal…
forward_CO1_14071      0      0      0      0      0      0      0      0      0 ""                                                 
forward_CO1_13891      0      0      0      0      0      0      0      0      0 Nematoda;Chromadorea;Monhysterida;Monhysteridae;NA…
forward_CO1_3833       0      1      0      0      0      0      0      0      0 NA;Oomycetes;Pythiales;Pythiaceae;Pythium;Pythium …

Anacapa outputs this information because the user can use this sequence ID to go back to a FASTA file and find the sequence; it is of no use right now to the ranacapa package. There are functions that exist within ranacapa that scrub out columns with seq_number in their name.

Eventually, ranacapa will also be able to accept qiime-formatted biom tables.

NOTE: As of now there is inconsistency in the terminology for this particular file type. I some times call it the "anacapa out" file and frequently use the variable name ana_out in functions to deal with it; sometimes, I might refer to it as the "biom table" or "biom file". One goal is to make the terminology more consistent.

2. The other input file that ranacapa cares about is called the "mapping file", which contains the metadata for all of the samples. For example, we might include the season in which the sample was collected, whether it was collected from under an Oak or in the grasslands, etc., in this file. The mapping file has (at least) as many rows of information as there are samples in the study, and one column per "metadata" variable. The first column of the mapping file must match the names of the samples (i.e. all non-taxonomy column names of the biom file. A valid mapping file looks like:

sample_name sample_id season texture depth main_plant
site_1         186B1       s      sand      10. other     
site_2         16A2        w      sand      20. other     
site_3         201A1       s      clay      10. oak       
site_4         201C2       s      sand      30. sage      
site_5         181C1       s      sand      30. oak       
site_6         PPM6mo      p      clay      30. other     
site_7         181B2       s      clay      10. sage      
site_8         182C1       s      clay      20. other     
site_9         15C2        w      silt      30. other     

rancapa under the hood

Apart from the shiny app, which does all the visualizations, most of ranacapa functions to convert between files that the users upload (i.e. the Anacapa output file and the associated Metadata file) and the phyloseq objects that ultimately get used for a lot of the analyses. This converting between formats happens in just a few core functions:

1. group_anacapa_by_taxonomy(): sometimes, in anacapa output, the same taxonomic path can appear in multiple rows- for example, if two different sequences had the same taxon identification. This function groups such rows together, and adds up the values in the sample columns. It's important because physeq objects can only have one row per identified taxonomy.

2. convert_anacapa_to_phyloseq() This function takes the output file and the mapping file as its inputs, and spits out a phyloseq class object that has both the taxonomic and the metadata information in the same object. To do this, the function does the folllowing:

• A. Convert the table into a phyloseq otu_table object
• B. Then, it extracts the rownames of the otu_table created in step 1 (the rownames contain the full taxonomic path), and splits it up on each semicolon (;) to make a new matrix with 6 columns, one for each of the taxonomic ranks. The matrix is then converted into a tax_table object.
• C. The otu_table and tax_table from steps A and B are united into the phyloseq style object.
• D. The information from the mapping file is converted to class sample_data
• E. The phyloseq object in step C is merged with the sample_data from step D with the phyloseq function merge_phyloseq .

3. vegan_otu() This function just takes a phyloseq style object, e.g. the one made by convert_anacapa_to_phyloseq, and extracts the otu table in the format of a community matrix that can be read by vegan functions.

4. A series of validator functions for the input files, currently only in the development branch.