Skip to content

Commit

Permalink
Docs
Browse files Browse the repository at this point in the history
  • Loading branch information
davidealbanese committed Jun 4, 2019
1 parent 05e7449 commit 9462b49
Showing 1 changed file with 69 additions and 22 deletions.
91 changes: 69 additions & 22 deletions doc/source/phyloseq.rst
Original file line number Diff line number Diff line change
@@ -1,10 +1,9 @@
An introduction to the downstream analysis with R and phyloseq
==============================================================

In this tutorial we describe what is needed to make some common analysis on the
data processed with micca using `R <https://www.r-project.org/>`_ and `phyloseq
<https://joey711.github.io/phyloseq/>`_. In particular, we will discuss the
following topics:
In this tutorial we describe a `R <https://www.r-project.org/>`_ pipeline for
the downstream analysis starting from the output of micca. In particular, we
will discuss the following topics:

- rarefaction;
- taxonomy and relative abundances;
Expand All @@ -16,14 +15,14 @@ following topics:

* This tutorial requires :doc:`pairedend_97` to be done.

* The tutorial is tested with R 3.5.3, phyloseq 1.26.1, ggplot2 3.1.0, vegan
* The tutorial is tested on R 3.5.3, phyloseq 1.26.1, ggplot2 3.1.0, vegan
2.5-4 and DESeq2 1.22.2.

Import data and preparation
---------------------------

We can import the micca processed data (the BIOM file, the phylogenetic tree and
the representative sequences) into the R environment using the ``import_biom()``
Import the micca processed data (the BIOM file, the phylogenetic tree and the
representative sequences) into the R environment using the ``import_biom()``
function available in `phyloseq <https://joey711.github.io/phyloseq/>`_ library.

.. code-block:: R
Expand All @@ -49,7 +48,52 @@ table (which contains the OTU counts for each sample), the sample data matrix
taxonomy for each OTU), the phylogenetic tree, and the OTU representative
sequences.

At this point, we can plot the rarefaction curves using vegan:
Print the metadata using the phyloseq function ``sample_data()``:

.. code-block:: R
> sample_data(ps)
Sample Data: [34 samples by 4 sample variables]:
Season Depth Month Year
B0214D1-PL1-D1 Winter 1 2 14
B0214D2-PL1-E1 Winter 10 2 14
B0214D3-PL1-F1 Winter 20 2 14
B0314D1-PL1-G1 Spring 1 3 14
B0314D2-PL1-H1 Spring 10 3 14
B0314D3-PL1-A2 Spring 20 3 14
B0414D1-PL1-B2 Spring 1 4 14
B0414D2-PL1-C2 Spring 10 4 14
B0414D3-PL1-D2 Spring 20 4 14
B0514D1-PL1-E2 Spring 1 5 14
B0514D2-PL1-F2 Spring 10 5 14
B0514D3-PL1-G2 Spring 20 5 14
B0614D1-PL1-H2 Summer 1 6 14
B0614D2-PL1-A3 Summer 10 6 14
B0714D2-PL1-B3 Summer 10 7 14
B0714D3-PL1-C3 Summer 20 7 14
B0814D1-PL1-D3 Summer 1 8 14
B0814D2-PL1-E3 Summer 10 8 14
B0814D3-PL1-F3 Summer 20 8 14
B0914D1-PL1-G3 Fall 1 9 14
B0914D2-PL1-H3 Fall 10 9 14
B0914D3-PL1-A4 Fall 20 9 14
B1014D1-PL1-B4 Fall 1 10 14
B1014D2-PL1-C4 Fall 10 10 14
B1014D3-PL1-D4 Fall 20 10 14
B1114D1-PL1-E4 Fall 1 11 14
B1114D2-PL1-F4 Fall 10 11 14
B1114D3-PL1-G4 Fall 20 11 14
B1214D1-PL1-H4 Winter 1 12 14
B1214D2-PL1-A5 Winter 10 12 14
B1214D3-PL1-B5 Winter 20 12 14
Bar0114D1-PL1-A1 Winter 1 1 14
Bar0114D2-PL1-B1 Winter 10 1 14
Bar0114D3-PL1-C1 Winter 20 1 14
The sample data contains 4 features for each sample: the season of sampling,
the sampling depth, the month and the year of sampling .

Plot the rarefaction curves using vegan function ``rarecurve()``:

.. code-block:: R
Expand All @@ -59,9 +103,13 @@ At this point, we can plot the rarefaction curves using vegan:
:align: center
:scale: 95%

Now we can rarefy the samples. Rarefaction is used to simulate even number of
reads per sample. In this example, the rarefaction depth chosen is 90% of the
minimum sample depth in the dataset (459 reads per sample).
``otu_table()`` is a phyloseq function which extract the OTU table from the
phyloseq object.

Rarefy the samples without replacement. Rarefaction is used to simulate even
number of reads per sample. In this example, the rarefaction depth chosen is the
90% of the minimum sample depth in the dataset (in this case 459 reads per
sample).

.. code-block:: R
Expand All @@ -75,7 +123,6 @@ minimum sample depth in the dataset (459 reads per sample).

* Remember to set the random seed (``rngseed``) for repeatable experiments.


.. admonition:: Exercise

Plot the samples depths before and after the rarefaction using the
Expand All @@ -85,9 +132,9 @@ minimum sample depth in the dataset (459 reads per sample).
Plot abundances
---------------

Using the rarefied dataset, make a stacked barplot of the abundances and color
each OTU (i.e. each bar) according its classified phylum (in this case
``Rank2``):
Using the rarefied dataset, make a stacked barplot of the abundances (read
counts) and color each OTU (i.e. each bar) according its classified phylum (in
this case ``Rank2``):

.. code-block:: R
Expand All @@ -111,8 +158,7 @@ according to the season:

Alternatively, we can merge the OTUs at the phylum level and build a new phyloseq
object. Given a taxonomic rank (in this case the phylum), the phyloseq function
``tax_glom`` merges the OTU with the same taxonomy, summing the relative
abundance values:
``tax_glom`` merges the OTUs with the same taxonomy, summing the abundances:

.. code-block:: R
Expand All @@ -126,24 +172,25 @@ abundance values:
refseq() DNAStringSet: [ 35 reference sequences ]
The option ``NArm`` set to ``FALSE`` forces the function to keep the
unclassified OTUs at the phylum level. In this case we obtain 35 phyla.
Now we can make a cleaner bar plot:
unclassified OTUs at the phylum level. Now we can make a cleaner bar plot:

.. code-block:: R
> plot_bar(ps.phylum, fill="Rank2") + facet_wrap(~Season, scales= "free_x", nrow=1)
.. image:: /images/garda_barplot3.png
:align: center
:scale: 75%

.. admonition:: Exercise

Make a stacked barplot at class level (``Rank3``).

Alpha diversity
---------------

Now we can plot the number of observed OTUs in each month, coloring the values
according to the sampling depth:
Plot the number of OTUs at each month and coloring the values according to
the sampling depth (``Depths``):

.. code-block:: R
Expand Down

0 comments on commit 9462b49

Please sign in to comment.