目录
- Rmatch script style1年前
- imagesAdded additional functions for plots and normalizations.5年前
- instfinished up final details2年前
- manadded requested changes2年前
- vignettesimplemented various requested changes2年前
- .Rbuildignorev0.18年前
- .gitignoreUpdates and golay correction added5年前
- DESCRIPTIONfinal adjustments2年前
- LICENSEAdded MIT8年前
- NAMESPACEadded requested changes2年前
- README.mdfinished up final details2年前
邀请码
版权所有：中国计算机学会技术支持：开源发展技术委员会
京ICP备13000930号-9
京公网安备 11010802032778号
Tutorial: Integrating QIIME2 and R for data visualization and analysis using qiime2R (v0.99.6)
Background
The qiime artifact is a method for storing the input and outputs for QIIME2 along with associated metadata and provenance information about how the object was formed. This method of storing objects has a number of obvious advantages; however, on the surface it does not lend itself to easy import to R for the R-minded data scientist. In reality, the .qza file is a compressed directory with an intuitive structure.
While it is possible to export data on a one-by-one basis from the qiime artifacts using qiime’s built in suite of export features this is problematic and runs antithetical to the purpose of the artifact format for the following reasons:
This package is trying to simplify the process of getting the artifact into R without discarding any of the associated data through a simple
read_qzafunction. The artifact is unpacked in to a temporary directory and the raw data and associated metadata are read into a named list (see below). Data are typically returned as either a data.frame, phylo object (trees), or DNAStringSets (nucleic acid sequences).Functions
read_qza()- Function for reading artifacts (.qza).qza_to_phyloseq()- Imports multiple artifacts to produce a phyloseq object.read_q2metadata()- Reads qiime2 metadata file (containing q2-types definition line)write_q2manifest()- Writes a read manifest file to import data into qiime2theme_q2r()- A ggplot2 theme for for clean figures.print_provenance()- A function to display provenance information.is_q2metadata()- A function to check if a file is a qiime2 metadata file.parse_taxonomy()- A function to parse taxonomy strings and return a table where each column is a taxonomic class.parse_ordination()- A function to parse the internal ordination format.read_q2biom()- A function for reading QIIME2 biom files in format v2.1make_clr()- Transform feature table using centered log2 ratio.make_proportion()- Transform feature table to proportion (sum to 1).make_percent()- Transform feature to percent (sum to 100).interactive_table()- Create an interactive table in Rstudio viewer or rmarkdown html.summarize_taxa()- Create a list of tables with abundances sumed to each taxonomic level.taxa_barplot()- Create a stacked barplot using ggplot2.taxa_heatmap()- Create a heatmap of taxonomic abundances using gplot2.corner()- Show top corner of a large table-like obejct.min_nonzero()- Find the smallest non-zero, non-NA in a numeric vector.mean_sd()- Return mean and standard deviation for plotting.subsample_table()- Subsample a table with or without replacement.filter_features()- Remove low abundance features by number of counts and number of samples they appear in.Installing qiime2R
qiime2R is currently available via github which can easily be installed in R via the following command:
Reading Artifacts (.qza)
This example is using data derived from the moving pictures tutorial. The main function we will need is
read_qza():We will start be reading in a table of sequence variants (SVs):
When the artifact is imported, there are a number of pieces of information included. To see them we can use the names command:
To access the actual data stored within the object, access the data as below:
In the above file each row denotes a sequence variant where in the actual text is the hash of the complete sequence. See here for an example tool for generating hashes.
We can also look at the unique identifier for this object:
We can see the type of artifact:
We can also get a complete list of the files within the artifact and their sizes:
We can also print the providence; however, it is probably easier to use the q2view tool for a graphical aid in its interpretation.
Reading Metadata
If you are using a qiime2 metadata file (outlined here), you can use the supplied function
read_q2metadata()as below. If using a standard tsv or csv file, you can useread.table(),readr::read_tsv(), orreadr::csv().Reading Taxonomy
Note, when taxonomy is imported, a single string is returned along with a confidence score. For many analysis we will want to break up this string and for that purpose the
parse_taxonomy()function is provided:Creating a Phyloseq Object
A wrapper function called
qza_to_phyloseq()is provided which links multipleread_qza()calls together to create a phyloseq object for subsequent analysis as per the phyloseq tutorials. An example usage is shown below:Creating a TreeSummarizedExperiment (TSE) Object
A wrapper function called
qza_to_tse()is provided which links multipleread_qza()calls together to create a TSE object for subsequent analysis. You may find some useful tutorials on: OMA tutorials. An example usage is shown below:Example Visualizations
In this next section, methods for publication-ready figures will be demonstrated through extensive use of the tidyverse, and particularly ggplot2 and its many extensions. If starting with R, learning these packages early on will have a bit of a learning curve, but pay off massively in the long term. Each section contains a self-contained code snippets capable of replicating the figures and links to the data used. A lot can be learned from modifying parameters of the graphing commands to customize and play around with the visualizations. Note, the pipe operator (%>%) is being used frequently here and passes the output of one line to the next. Learn more about piping in R here.
All the data we will use is from the “Moving Pictures” tutorial and/or related content which can be found here.
Alpha Diversity Over Time
Data:
This dataset provides a good example of why you should always look at your data. There are an extra 3 samples in the metadata which do not have an assigned Shannon diversity value as seen below:
If this was your dataset, you should look into why this was and try to see why these samples have been lost. We will continue after sorting them out:
Note how we have now added a column for the shannon diversity to our table which will make everything easy to plot. Given this dataset is longitudinal we can plot it as such. Note: the code I am using here would handle replicate data and plot error bars; however, this case does not have replicates. Also, a quick note is that R does generally not like
-s in column names. To use these names we can surround them with ticks.Imagine we wanted to instead ask if there was an effect of antibiotics on diversity. We can change the type of plot:
Maybe instead the question is if the two individuals have different diversity?
Note: It is important to remember that these a graphical tools, but to make a claim, you should also generate some form of statistical support taking the nature of the data into consideration (ie that it uses repeated sampling over time).
Plotting PCoA
Here I am going to demo a PCoA with extra metadata mapped to it. I indicate the body site, antibiotic usage, and the shannon diversity by mapping to different aesthetic parameters.
Data:
Plotting a Heatmap
Data:
Heat maps provide a very good overview of the composition of samples but tend to not perform well for complex samples or visualizing small differences. qiime2R includes a function for making heatmaps called
taxa_heatmap().Making a taxonomic barplot
Data:
Heat maps provide a very good overview of the composition of samples but tend to not perform well for complex samples or visualizing small differences. qiime2R includes a function for making heatmaps called
taxa_heatmap().Differential Abundance Analysis
There are many ways to find features that are deferentially abundant across conditions. In this case we are going to use the output of q2-aldex as described here.
Volcano Plot
Plotting Per-Feature Abundances
If we only wanted to plot one significantly different feature, say 4b5eeb300368260019c1fbc7a3c718fc, we can plot only its abundances. The abundances are not exported by ALDEx2 so instead we will recalculate the CLR-abundances using an approximation of ALDEx; however, in R this could be easily calculated using the package itself.
Plotting a Phylogenetic Tree
Finally, we can visualize the results on a phylogenetic tree with the help of ggtree.
Note: please report problems or suggestions on the qiime2 forum or the github issue tracker.