d

Author: WisamSaleem

.Rhistory

@@ -1,45 +1,3 @@
-labs(y = "Diversity", x = "", title = paste0("Kruskal test p=", round(pv, 2)))
-df <- meta(phy)
-df <- bind_cols(df, A)
-df$index <- df[[index]]
-pv <- kruskal.test(data = df, index ~ factor(bmi_group))$p.value
-p1 <- ggplot(df, aes(x = bmi_group, y = index, fill = bmi_group)) + geom_boxplot() +
-labs(y = "Diversity", x = "", title = paste0("Kruskal test p=", round(pv, 2)))
-print(p1)
-p2 <- ggplot(df, aes(fill = bmi_group, x = index)) +
-geom_density(alpha = 0.5) + labs(x = "Diversity", y = "Density")
-print(p2)
-p1 <- ggplot(df, aes(x = bmi_group, y = index, fill = bmi_group)) + geom_boxplot() +theme(axis.line.x = element_text(angle = 60, hjust = 1)) +
-labs(y = "Diversity", x = "", title = paste0("Kruskal test p=", round(pv, 2)))
-df <- meta(phy)
-df <- bind_cols(df, A)
-df$index <- df[[index]]
-pv <- kruskal.test(data = df, index ~ factor(bmi_group))$p.value
-p1 <- ggplot(df, aes(x = bmi_group, y = index, fill = bmi_group)) + geom_boxplot() + theme(axis.text.x = element_text(angle = 60, hjust = 1)) +
-labs(y = "Diversity", x = "", title = paste0("Kruskal test p=", round(pv, 2)))
-print(p1)
-p2 <- ggplot(df, aes(fill = bmi_group, x = index)) +
-geom_density(alpha = 0.5) + labs(x = "Diversity", y = "Density")
-print(p2)
-library(dplyr)
-source("functions.R")
-library(dplyr)
-#source("functions.R")
-phy.clr <- phy %>% microbiome::transform("compositional")
-genus.stress.clr <- phy.clr %>% abundances() %>% t
-permutations <- 999
-controls <- c()
-independent_vars <- "Group" # Groupings to compare
-alpha_level <- 0.05
-library(vegan)
-library(magrittr)
-library(tibble)
-permanova_results <- lapply(independent_vars, function(var) do_permanova(pseq = phy.clr,
-group = var,
-permutations = permutations)) %>%
-set_names(independent_vars) %>%
-do.call(rbind, .) %>%
-m_neat(colnames = c("variable", "R2", "p_value")) %>%
 mutate(padj_value = p.adjust(p_value, "BH")) %>%
 arrange(variable)
 library(dplyr)
@@ -510,3 +468,45 @@ rmarkdown::render('continuous_variable.rmd','md_document')
 getwd()
 setwd('D:/Seafile/My Library/Important/Final_microbiome/tutorials')
 rmarkdown::render('Tutorial.rmd','all')
+# Load libraries
+library(microbiome)
+library(ggplot2)
+library(magrittr)
+library(dplyr)
+# Probiotics intervention example data
+data(dietswap)
+# Only check the core taxa to speed up examples
+pseq <- core(dietswap, detection = 50, prevalence = 80/100)
+library(phyloseq)
+library(reshape2)
+library(DESeq2)
+library(knitr)
+library(magrittr)
+# Running the DESeq2 analysis
+ds2 <- phyloseq_to_deseq2(pseq, ~ nationality)
+dds <- DESeq(ds2)
+res <- results(dds)
+df <- as.data.frame(res)
+df$taxon <- rownames(df)
+df <- df %>% arrange(log2FoldChange, padj)
+library(knitr)
+print(head(kable((df))))
+# Identify top taxa based on standard ANOVA
+source(system.file("extdata/check_anova.R", package = "microbiome"))
+ano <- check_anova(pseq, "nationality");
+ano$log2FC <- log2(ano$ave.AFR) - log2(ano$ave.AAM)
+taxa.anova <- as.character(subset(ano, padj < 0.01 & abs(log2FC) > log2(2))$taxa)
+# lowPick the top taxa based on DESEq2
+taxa.deseq <- subset(res, padj < 0.01 & abs(log2FoldChange) > log2(2))$taxon
+# Check overlap
+# Most DESEq2 taxa are confirmed with ANOVA
+library(gplots)
+#venn( list(ANOVA = taxa.anova,DESeq2 = taxa.deseq) )
+# Also the est p-values are well correlated (higher not so)
+mf<-data.frame(df$padj, ano$padj)
+p<-ggplot(mf, aes(x = log10(df$padj), y = log10(ano$padj))) + xlab('DESeq2 adjusted p-value') + ylab('ANOVA adjusted p-value')
+p<- p + geom_point()
+#plot(xlim=log10(df$padj), log10(ano$padj), type = 'l', xlim = "DESeq2 adjusted p-value",  ylim("ANOVA adjusted p-value"))
+venn( list(ANOVA = taxa.anova,DESeq2 = taxa.deseq))
+print(p)
+rmarkdown::render('deseq2.Rmd','all')

Atlas.html

@@ -135,7 +135,7 @@ p <- plot_core(pseq, plot.type = "heatmap",
        		 prevalences = prevalences,
        		 detections = detections,
 		 colours = rev(brewer.pal(5, "Spectral")),
-		 min.prevalence = .2, horizontal = TRUE)
+		 min.prevalence = .2, horizontal = TRUE) 
 print(p)
 ```
 
@@ -323,7 +323,7 @@ detections <- 10^seq(log10(1e-5), log10(.2), length = 10)
 gray <- gray(seq(0,1,length=5))
 p <- plot_core(pseq.rel, plot.type = "heatmap", colours = gray,
     prevalences = prevalences, detections = detections, min.prevalence = .5) +
-    xlab("Detection Threshold (Relative Abundance (%))")
+    labs(x = "Detection Threshold (Relative Abundance (%))")
 print(p)    
 
 
@@ -336,59 +336,3 @@ library(viridis)
 print(p + scale_fill_viridis())
 ```
 
-As it can be seen, we see only OTu IDs and this may not be useful to interpret the data. We need to repreoccess this figure to include taxonomic information. We can do this as follows:  
-
-```{r,core-example4, fig.width=14, fig.heigth=18, fig.show='hold', out.width = '200px', warning=FALSE, eval=FALSE}
-library(RColorBrewer)
-library(knitr)
-# Core with absolute counts and vertical view:
-# and minimum population prevalence (given as percentage)
-detections <- 10^seq(log10(1), log10(max(abundances(pseq.2))/10), length = 10)
-
-healthycore <- plot_core(pseq.2, plot.type = "heatmap", 
-       		 prevalences = prevalences,
-       		 detections = detections,
-		 colours = rev(brewer.pal(5, "Spectral")),
-		 min.prevalence = .90, horizontal = F)
-# get the data used for plotting 
-df <- healthycore$data 
-
-# get the list of OTUs
-list <- df$Taxa 
-
-# check the OTU ids
-# print(list) 
-
-# get the taxonomy data
-tax <- tax_table(pseq.2)
-tax <- as.data.frame(tax)
-
-# add the OTus to last column
-tax$OTU <- rownames(tax)
-
-# select taxonomy of only 
-# those OTUs that are used in the plot
-tax2 <- dplyr::filter(tax, rownames(tax) %in% list) 
-
-# head(tax2)
-
-# We will merege all the column into one except the Doamin as all is bacteria in this case
-tax.unit <- tidyr::unite(tax2, Taxa_level,c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species", "OTU"), sep = "_;", remove = TRUE)
-
-tax.unit$Taxa_level <- gsub(pattern="[a-z]__",replacement="", tax.unit$Taxa_level)
-
-# add this new information into the plot data df
-
-df$Taxa <- tax.unit$Taxa_level
-
-# you can see now we have the taxonomic information
-knitr::kable(head(df))
-
-# replace the data in the plot object
-healthycore$data <- df
-
-plot(healthycore + theme(axis.text.y = element_text(face="italic")))
-```
-
-
-

Betadiversity.html

@@ -46,7 +46,7 @@ This returns a table with selected diversity indicators. Check a separate page o
 
 ps1 <- prune_taxa(taxa_sums(pseq) > 0, pseq)
 
-tab <- alpha(ps1, index = "all")
+tab <- microbiome::alpha(ps1, index = "all")
 kable(head(tab))
 
 ```

Bimodality.html

@@ -158,7 +158,9 @@ data(dietswap)
 x <- dietswap
 # Compositional data 
 x2 <- microbiome::transform(x, "compositional")
+```
 
+```{r data10bb, echo=FALSE, warning=FALSE, message=FALSE, results="hide"}
 # library(zCompositions)
 # This would be recommended for compositional zero-inflated data:
 # KMSS replacement for zero-inflated compositional values
@@ -186,7 +188,7 @@ Assign new fields to metadata
 
 ```{r preprocess4b2, warning=FALSE, message=FALSE}
 # Calculate diversity for samples
-div <- alpha(pseq, index = "shannon")
+div <- microbiome::alpha(pseq, index = "shannon")
 
 # Assign the estimated diversity to sample metadata
 sample_data(pseq)$diversity <- div

Clustering.html

@@ -382,13 +382,13 @@ <h3>Contribute</h3>
 <ul>
 <li><a href="https://github.com/microbiome/microbiome/issues">Issue Tracker</a></li>
 <li><a href="https://github.com/microbiome/microbiome/">Pull requests</a></li>
-<li>Subscribe to the <a href="https://groups.google.com/forum/#!forum/microbiome-devel">mailing list</a> (<a href="mailto:[email protected]" class="email">[email protected]</a>)</li>
+<li>Subscribe to the <a href="https://groups.google.com/forum/#!forum/microbiome-devel">mailing list</a> (<a href="mailto:[email protected]">[email protected]</a>)</li>
 <li><a href="https://github.com/microbiome/microbiome">Star us on the Github page</a></li>
 </ul>
 </div>
 <div id="publications-using-the-microbiome-package" class="section level3">
 <h3>Publications using the microbiome package</h3>
-<p><strong>Kindly cite this work</strong> as follows: &quot;Leo Lahti <a href="https://github.com/microbiome/microbiome/graphs/contributors">et al.</a> (2017). Tools for microbiome analysis in R. Microbiome package version 1.6.0. URL: <a href="https://microbiome.github.io/tutorials">http://microbiome.github.com/microbiome</a>, and the relevant references listed in the manual page of each function. The list of publications is not exhaustive. Let us know if you know of further publications using the microbiome package; we are collecting these on the website.</p>
+<p><strong>Kindly cite this work</strong> as follows: &quot;Leo Lahti <a href="https://github.com/microbiome/microbiome/graphs/contributors">et al.</a> (2017). Tools for microbiome analysis in R. Microbiome package version 1.9.1. URL: <a href="https://microbiome.github.io/tutorials">http://microbiome.github.com/microbiome</a>, and the relevant references listed in the manual page of each function. The list of publications is not exhaustive. Let us know if you know of further publications using the microbiome package; we are collecting these on the website.</p>
 <p><a href="https://academic.oup.com/femsre/article/doi/10.1093/femsre/fuw045/2979411/Intestinal-microbiome-landscaping-insight-in#58802539">Intestinal microbiome landscaping: Insight in community assemblage and implications for microbial modulation strategies</a>. Shetty S, Hugenholtz F, Lahti L, Smidt H, de Vos WM, Danchin A. <em>FEMS Microbiology Reviews</em> fuw045, 2017.</p>
 <p><a href="http://dx.doi.org/10.1016/j.mib.2015.04.004">Metagenomics meets time series analysis: unraveling microbial community dynamics</a> Faust K, Lahti L, Gonze D, de Vos WM, Raes J. <em>Current Opinion in Microbiology</em> 15:56-66 2015.</p>
 <p><a href="http://www.nature.com/ncomms/2014/140708/ncomms5344/full/ncomms5344.html">Tipping elements in the human intestinal ecosystem</a> Lahti L, Salojärvi J, Salonen A, Scheffer M, de Vos WM. <em>Nature Communications</em> 5:4344, 2014.</p>

Comparisons.html

@@ -18,6 +18,10 @@ List of issues to improve the package and tutorials:
 * Lay out overall project roadmap
 * Consider package/tutorial collection instead of a single package
 
+### Missing examples
+
+ * is_compositional
+
 ### Related work
 
  * [metagenomeSeq](https://bioconductor.org/packages/release/bioc/html/metagenomeSeq.html)

Composition.html

@@ -318,6 +318,7 @@ Long tails of rare taxa:
 
 ```{r tail, warning=FALSE, message=print, fig.width=16, fig.height=8, out.width="200px"}
 medians <- apply(abundances(d),1,median)/1e3
+library(reshape2)
 A <- melt(abundances(d))
 A$Var1 <- factor(A$Var1, levels = rev(names(sort(medians))))
 p <- ggplot(A, aes(x = Var1, y = value)) +

CompositionAmplicondata.html

@@ -0,0 +1,56 @@
+As it can be seen, we see only OTu IDs and this may not be useful to interpret the data. We need to reprccess this figure to include taxonomic information. We can do this as follows:  
+
+```{r,core-example4, fig.width=14, fig.heigth=18, fig.show='hold', out.width = '200px', warning=FALSE, eval=FALSE}
+library(RColorBrewer)
+library(knitr)
+# Core with absolute counts and vertical view:
+# and minimum population prevalence (given as percentage)
+detections <- 10^seq(log10(1), log10(max(abundances(pseq.2))/10), length = 10)
+
+healthycore <- plot_core(pseq.2, plot.type = "heatmap", 
+       		 prevalences = prevalences,
+       		 detections = detections,
+		 colours = rev(brewer.pal(5, "Spectral")),
+		 min.prevalence = .9)
+# get the data used for plotting 
+df <- healthycore$data 
+
+# get the list of OTUs
+list <- df$Taxa 
+
+# check the OTU ids
+# print(list) 
+
+# get the taxonomy data
+tax <- tax_table(pseq.2)
+tax <- as.data.frame(tax)
+
+# add the OTus to last column
+tax$OTU <- rownames(tax)
+
+# select taxonomy of only 
+# those OTUs that are used in the plot
+tax2 <- dplyr::filter(tax, rownames(tax) %in% list) 
+
+# head(tax2)
+
+# We will merege all the column into one except the Doamin as all is bacteria in this case
+tax.unit <- tidyr::unite(tax2, Taxa_level,c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species", "OTU"), sep = "_;", remove = TRUE)
+
+tax.unit$Taxa_level <- gsub(pattern="[a-z]__",replacement="", tax.unit$Taxa_level)
+
+# add this new information into the plot data df
+
+df$Taxa <- tax.unit$Taxa_level
+
+# you can see now we have the taxonomic information
+knitr::kable(head(df))
+
+# replace the data in the plot object
+healthycore$data <- df
+
+plot(healthycore + theme(axis.text.y = element_text(face="italic")))
+```
+
+
+