Documentation updated

doc/source/phyloseq.rst

@@ -8,8 +8,11 @@ An introduction to the downstream analysis with R and phyloseq
 .. note::
 
    This tutorial requires `R <https://www.r-project.org/>`_, `phyloseq
-   <https://joey711.github.io/phyloseq/>`_ and ggplot2 (tested on R v3.4. and
-   phyloseq v1.22.3) to be installed in your system.
+   <https://joey711.github.io/phyloseq/>`_ ggplot2 and vegan (tested on R v3.4
+   and phyloseq v1.22.3) to be installed in your system.
+
+Import data and rarefaction
+---------------------------
 
 We can import the micca processed data (the BIOM file, the phylogenetic tree and
 the representative sequences) into the `R <https://www.r-project.org/>`_
@@ -23,9 +26,9 @@ file.
 
     > library("phyloseq")
     > library("ggplot2")
+    > library(vegan)
     > setwd("denovo_greedy_otus") # set the working directory
-    > ps = import_biom("tables.biom", treefilename="tree_rooted.tree",
-    + refseqfilename="otus.fasta")
+    > ps = import_biom("tables.biom", treefilename="tree_rooted.tree", refseqfilename="otus.fasta")
     > sample_data(ps)$Month <- as.numeric(sample_data(ps)$Month)
     > ps
     phyloseq-class experiment-level object
@@ -35,28 +38,93 @@ file.
     phy_tree()    Phylogenetic Tree: [ 529 tips and 528 internal nodes ]
     refseq()      DNAStringSet:      [ 529 reference sequences ]
 
-At this point, we can compute the number of OTUs as measure of alpha diversity
-after the rarefaction:
+At this point we can plot the rarefaction curves using vegan:
+
+.. code-block:: R
+
+    > rarecurve(t(otu_table(ps)), step=50, cex=0.5)
+
+.. image:: /images/garda_rarecurves.png
+    :align: center
+    :scale: 95%
+
+Now, we can rarefy in order to bring the samples to the same depth:
 
 .. code-block:: R
 
     > # rarefy without replacement
     > ps.rarefied = rarefy_even_depth(ps, rngseed=1, sample.size=0.9*min(sample_sums(ps)), replace=F)
-    > # plot the number of observed OTUs
+   
+Plot abundances
+---------------
+
+.. code-block:: R
+
+    > plot_bar(ps.rarefied, fill="Rank2") + facet_wrap(~Season, scales = "free_x", nrow = 1)
+
+.. image:: /images/garda_bar.png
+    :align: center
+    :scale: 75%
+
+
+Alpha diversity
+---------------
+
+Now we can plot the number of observed OTUs
+
+.. code-block:: R
+
     > plot_richness(ps.rarefied, x="Month", color="Depth", measures=c("Observed"))
 
 .. image:: /images/garda_alpha.png
     :align: center
     :scale: 75%
 
-Finnaly, we can plot the PCoA using the unweighted UniFrac as distance:
+Moreover, we can make a boxplot of the number of OTUs and the Shannon entropy 
+grouping the different months by season:
+
+.. code-block:: R
+
+    > plot_richness(ps.rarefied, x="Season", measures=c("Observed", "Shannon")) + geom_boxplot()
+
+.. image:: /images/garda_alpha2.png
+    :align: center
+    :scale: 75%
+
+Beta diversity
+--------------
+
+Now, we can plot the PCoA using the unweighted UniFrac as distance:
 
 .. code-block:: R
 
     > # PCoA plot using the unweighted UniFrac as distance
-    > ordination = ordinate(ps.rarefied, method="PCoA", distance="unifrac", weighted=F)
+    > wunifrac_dist = distance(ps.rarefied, method="unifrac", weighted=F)
+    > ordination = ordinate(ps.rarefied, method="PCoA", distance=wunifrac_dist)
     > plot_ordination(ps.rarefied, ordination, color="Season") + theme(aspect.ratio=1)
 
 .. image:: /images/garda_beta.png
     :align: center
     :scale: 75%
+
+Now we test whether the seasons differ significantly from each other using the 
+permutational ANOVA (PERMANOVA) analysis:
+
+.. code-block:: R
+
+    > adonis(wunifrac_dist~sample_data(ps.rarefied)$Season)
+    
+    Call:
+    adonis(formula = wunifrac_dist ~ sample_data(ps.rarefied)$Season) 
+
+    Permutation: free
+    Number of permutations: 999
+
+    Terms added sequentially (first to last)
+
+                                    Df SumsOfSqs  MeanSqs F.Model     R2 Pr(>F)    
+    sample_data(ps.rarefied)$Season  3    0.6833 0.227765  4.3451 0.3029  0.001 ***
+    Residuals                       30    1.5726 0.052419         0.6971           
+    Total                           33    2.2559                  1.0000           
+    ---
+    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1