1. Thanks for quickly introducing the issues I am reporting.
2. Here is another thing with the plot_composition functon.

When I am using without x.label optiion everything works fine:

plot_composition(phyl1,taxonomic.level = "Phylum",relative.abundance = T,sort.by = "lesion")

but with this option the labels are gone:

plot_composition(phyl1,taxonomic.level = "Phylum",relative.abundance = T,sort.by = "lesion",x.label = "lesion")

Remove log10 option from hotplot and instead handle that in the ggplot scale_x_log10 for clarity

Hi, I'm using the following to produce a core heat map:

library(microbiome)
pseq.core <- core(relative.abundance, detection = 0.001, prevalence = 95/100)
pseq.core
tax_table(pseq.core) #to get names for comparison (need to go through individually)
#total core abundance in each sample
core.abundance <- sample_sums(pseq.core)
core.abundance
write.table(as.data.frame(core.abundance), file="core.abundance.txt")

#Set different detection levels and prevalence
prevalences <- seq(.5, 1, .5) #0.5 = 95% prevalence
detections <- 10^seq(log10(1e-3), log10(.2), length = 10) #(1e-3) = 0.001% abundance; change "-3" to -2 to increase to 0.01%
p <- plot_core(pseq.core, plot.type = "heatmap", colours = rev(brewer.pal(5, "Spectral")),
min.prevalence = 0.01, prevalences = prevalences, detections = detections) +
xlab("Detection Threshold (Relative Abundance (%))")
print(p)

but the graph prints with the SV names (i.e. SV1, SV2, SV3 etc) rather than the taxonomic information (i.e. genus). Is it possible to coerce these names into the figure?
Also, is there a way to agglomerate the data at a higher taxonomic level (e.g. family or order) and then look at the core microbiome of that?

Many thanks in advance

in the syntax,

boxplot_abundance <- function (pseq, x, y, line = NULL, color = NULL, log10 = TRUE, title = NULL)

if I run the function boxplot_abundance function independently (without the package), y does not take a Family/Genus "f__Prevotellaceae" name as suggested in your vignette. Rather takes an OTUID for y. Is that normal?

Add usage examples of phylogenetic tree aware methods (tree based discriminant analysis, agPCA, phylofactor for phylogeny based factorization, structSSI for hierarchical multiple testing correction etc) in the online tutorial to demonstrate compatibility with other tools and to give reproducible usage examples.

As with plot_composition. plot_core gives OTU id instead of taxonomic identities. Tested on DynamicsIBD dataset. The rownames in otu_table when out of seq based approach has OTU ids and hence the functions return OTU IDs. This may need

library(microbiome)
data("DynamicsIBD")
pseq <- DynamicsIBD
pseq
colnames(tax_table(pseq)) <- c("Domain","Phylum", "Class", "Order", "Family", "Genus", "Species")
colnames(tax_table(pseq))

pseq.1 <- subset_samples(pseq, ibd_subtype == "HC")
pseq.1
pseq.rel <- microbiome::transform(pseq.1, "compositional")
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-3), log10(.2), length = 10)

library(RColorBrewer)
healthycore <- plot_core(pseq.rel, plot.type = "heatmap", colours = rev(brewer.pal(5, "Spectral")),
               prevalences = prevalences, detections = detections, min.prevalence = .5) +
  xlab("Detection Threshold (Relative Abundance (%))")
print(healthycore)

Hi,
I was trying the potential analysis and got the following errors:

1. diversity <- exp(microbiome::diversity(pseq)$shannon)

Here, I think microbiome::diversity needs to be corrected to microbiome::diversities

2. By running res <- movpotential_ews(diversity, age), I got...

Error in maxima[[j]] : subscript out of bounds
In addition: Warning messages:
1: In movpotential_ews(diversity, age) :
The data contains NAs, removing the associated samples from X and param input arguments.
2: In min(maxima) : no non-missing arguments to min; returning Inf
3: In max(maxima) : no non-missing arguments to max; returning -Inf

I found NAs in age and tried to run the command again after removing them, but didn't work out.

Thank you in advance!

@antagomir the plot_composition function has changed and the arguments in tutorial no longer exists. I would prefer the previous function with options for taxonomic.level, transform in this function as colleagues at MIB found it more useful.

Cheers,
Sudarshan

'summarize_phyloseq' object to give like 'biom' summary.

1. Total number of reads,
2. Average number of reads,
3. Minimum library size,
4. Maximum library size,
5. Percent classification of OTUs to various taxonomic levels,
6. Normality test for library size,
7. Number of singletons.

The rda_bagged function is missing some info on the outputs (@return field in roxygen) The documentation could be revised before Bioc release. The function was renamed and simplified, and only the main function exported in order to streamline the package for upcoming Bioc submission. The I.thrsh argument was removed as I did not see it being used anywhere; the ncore argument was removed as I bumped into problems in modifying the functions and it did not seem critical feature for the first phase. The earlier versions are available and this can be extended later with some of that code.

See also the RDA documentation page, which can be modified from the corresponding Rmd file.

For functions read_biom2phyloseq, read_csv2phyloseq, read_mothur2phyloseq etc: add runnable examples to Roxygen by addin small examples files in inst/extdata and using system.file("inst/extdata/..", package = "microbiome") to read with the function.

Hi,

I used to have an installed version of mcirobiome in RStudio which stopped working. I could never invoke the library or the functions in it. So I removed the folder from the win-library directory. But I am now unable to install it again as I keep getting the below error.

install_github("microbiome/microbiome")
Downloading GitHub repo microbiome/microbiome@master
from URL https://api.github.com/repos/microbiome/microbiome/zipball/master
Installation failed: cannot open file 'C:/Users/yugandhar.reddybs/AppData/Local/Temp/Rtmp2JVzOP/devtools1844756e4631/microbiome-microbiome-a708087/inst/extras/ghpages/inst/extras/microbiome.Rcheck/00_pkg_src/microbiome/vignettes/Betadiversity_files/figure-markdown_strict/heterogeneity-example2bbb-1.png': No such file or directory

my system information is as below:

version
_
platform x86_64-w64-mingw32
arch x86_64
os mingw32
system x86_64, mingw32
status
major 3
minor 3.1
year 2016
month 06
day 21
svn rev 70800
language R
version.string R version 3.3.1 (2016-06-21)
nickname Bug in Your Hair

How can I get over this problem?

Thanks in advance,
Yugandhar

Feature request: make this package available via bioconda as has been done with phyloseq.

Due to recent changes in how r-base is compiled for bioconda or conda-forge, R package installations using install.packages() and biocLite() with R installed via conda will now often fail because R cannot find the correct compilers. This happens because the compilers must be in the conda environment, but install.packages() and biocLite() won't know if they are in the environment or not.

Due to this issue, if R is installed via conda (in order to easily incorporate R + Jupyter), R packages like microbiome often cannot be installed, because the package isn't available via conda and biocLite() will likely fail (unless the proper compilers are installed into the conda environment). Thus, anyone using Jupyter for reproducible bioinformatics is left out if the R package isn't available via conda. Anyone can contribute a package to bioconda, but it's probably best if the package developer is in charge of it.

According to some feedback we have received, the phyloseq package does not provide full capabilities for mothur (and qiime?) imports, and updates to the package may be slow. We were asked to expand the read_phyloseq function so that it can also import simple Mothur OTU + Taxonomy tables in addition to simple CSV:s. That format is also different from the standard Mothur import of the phyloseq package and would hence provide a useful complementary tool. Also other simple table formats could be added as an option.

Following this article add evenness metrics that perform well according to a series of statistical requirements, including lacking a strong bias toward very large or very small abundances, independence of richness (S), and scaling between zero (no evenness) and one (perfect evenness). These included Smith and Wilson’s indices (Evar and EQ), Simpson’s evenness (E1/D), Bulla’s index (O), and Camargo’s index (E′) [See Magurran AE, McGill BJ , eds (2011) Biological Diversity: Frontiers in Measurement and Assessment (Oxford Univ Press, Oxford), Vol 12 as well as http://www.jstor.org/stable/3545749].

Hi,
I just wanted to highlight that there should be an option to sort the samples according to more than one factor.

Lets assume that I have two variables "disease" and "lesion". I would like my samples to be sorted now according to both of these variables meaning first I want the lesional samples to be grouped and there can be two groups Lesional or non-lesional, then I want the non-lesional to be grouped again into two subcategories: "disease+" and "disease-".

As it is imposible for a "disease-" sample to fall into the lesional category it will work perfectly.

Hope this clarifies what I mean.

Best,
W

It would handy to have an indicator in phyloseq class definition that would indicate whether the data is compositional, or absolute counts (or perhaps some other transformation). This may affect how certain downstream analyses are carried out.

As in the title

Would be good to complement the write_phyloseq function with a read_phyloseq that knows how to read such data back in R. What you think ?

The read_phyloseq, read_biom2phyloseq etc useq datasets from QIITA and DynamicsIBD. We had to remove these data sets from the package due to Bioc constraints for package size (4Mb). I have commented out the respective parts in the documentation (Rd files as well as docs/Rmd files) but these need to be updated. One solution is to make a separate example data repository online together with the tutorial. The larger data sets can then be fetched from there on-the-fly, This will require access to internet connection however.

Hi Antagomir,

Thanks for developing the package,

I would like to be able to plot boxplot_abudance on paired fecal and biopsy samples (paired as in they are obtained from the same participant). What happens if that particular OTU occurs in fecal but does not occur in biopsy, right now it's pairing multiple samples in the biopsy group to just one sample in the fecal group (which is not the exact corresponding pair)
Can you help me interpret the plot better? I would like that each biopsy sample is paired with its corresponding fecal counterpart at 0 even if that OTU is absent in the plot.
Do multiple lines mean that there are multiple overlapping points at that same 'signal' value? Because it happens at non-zero signal values as well.

Also, is the lower bound of Y axis set to 0? I tried to set the limits in ggplot2 it interfered with the scale10 transformation. If it can show 0 that will be helpful

Finally, in the syntax,
boxplot_abundance <- function (pseq, x, y, line = NULL, color = NULL, log10 = TRUE, title = NULL)

if I run the function boxplot_abundance function independently (i.e without the package, sorry I have set a version control since I have multiple dissertation analyses with various dependencies and cannot update all packages for successful microbiome installation currently), y does not take a Family/Genus "f__Prevotellaceae" name as suggested in your vignette. Rather takes an OTUID for y. Is that normal?

Thanks a lot for your time again

Hi- The was an error in the 'Error in travis build' for the updates i did as suggested previously. Could you please help me debug this.
https://travis-ci.org/microbiome/microbiome/builds/210866256

Cheers,
Sudarshan

I would like to be able to plot boxplot_abudance on paired fecal and biopsy samples (paired as in they are obtained from the same participant). What happens if that particular OTU occurs in fecal but does not occur in biopsy, right now it's pairing multiple samples in the biopsy group to just one sample in the fecal group (which is not the exact corresponding pair)
Can you help me interpret the plot better? I would like that each biopsy sample is paired with its corresponding fecal counterpart at 0 even if that OTU is absent in the plot.
Do multiple lines mean that there are multiple overlapping points at that same 'signal' value? Because it happens at non-zero signal values as well.

Simplify the transform function to make more modification easier + implementation faster. This is coupled with the abundances function.

R CMD check gives:

• checking R code for possible problems ... NOTE
  read_biom2phyloseq: no visible binding for global variable
  ‘parse_taxonomy_default’
  Undefined global functions or variables:
  parse_taxonomy_default

The axis labels in bRDA plot for a phyloseq object are sometimes NA, sometimes extremely small (1e-5..), and sometimes thousands of percents. Not yet tested whether this is an issue with matrix input as well, or limited to phyloseq. Here one example:

library(microbiome)
data(peerj32)
plot_rda_bagged(rda_bagged(peerj32$phyloseq, "group", bs.iter=10))

When you create new tutorial pages, please add a link to vignettes/index.Rmd - this helps to organize the tutorials thematically.

Please clarify if this is intentional and propose a possible workaround.

In order to easily visually compare samples, I would suggest to plot always the OTUs in the same "level" of a plot.

The way plot_composition() works right now is to plot the OTUs based on their abundance: the most abundant at the bottom of the bar and the least abundant at the top.

I think there should be a possibility to override this. So that I always would get OTU nr 1 on the bottom and OTU 100 on top. based on their abundance in the whole database (not just restricted to one sample)

When I plot with phyloseq i get this:

PLEASE LOOK AT THE TWO HIGHEST ABUNDANCES and the CORRESPONDING SAMPLE NAMES

and when I plot with plot_composition() I get this:

Notice that the sample names are totally incorrect!

Hi, I'm trying to make a heat map of the core microbiome, using the following code:

library(microbiome)
pseq <- final.4
meta <- meta(pseq)
taxonomy <- tax_table(pseq)
pseq.rel <- microbiome::transform(pseq, "compositional")
colnames(tax_table(pseq.rel)) <- c("Domain","Phylum", "Class", "Order", "Family", "Genus")
colnames(tax_table(pseq.rel))
pseq = merge_phyloseq(pseq, tax_table(final.4))

#Simple anlaysis
core.taxa.standard <- core_members(pseq.rel, detection = 0, prevalence = 50/100)
core.taxa.standard
core.taxa.standard = merge_phyloseq(core.taxa.standard, tax_table(pseq))
core.taxa.standard

#Creates phyloseq object for core microbiome
pseq.core <- core(pseq.rel, detection = 0, prevalence = .5)
pseq.core
pseq.core = merge_phyloseq(pseq.core, tax_table(pseq))
pseq.core

#Use relative abundance
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-3), log10(.2), length = 10)
p <- plot_core(pseq.rel, plot.type = "heatmap", colours = rev(brewer.pal(5, "Spectral")),
prevalences = prevalences, detections = detections) +
xlab("Detection Threshold (Relative Abundance (%))")
print(p)

However, the resulting figure just gives SV's on the y axis, rather than the actual genus name. Do you know how to coerce the proper taxa information back into this object? Some of the code above might be redundant - it was my attempt to fix this problem!

Many thanks

would be useful functions sometimes.. to be implemented.

Format the code so that line widths are <80 characters (Bioc style requirement).

Also do for other functions maintained.

Merge the package with microbiomeSeq

https://academic.oup.com/bioinformatics/article/33/8/1179/2907823/Scater-pre-processing-quality-control

This article addresses pre-processing-quality-control which is crucial for microbiome data also. Maybe we can try to incorporate this in the next version update?

boxplot_abundance Also, is the lower bound of Y axis set to 0? I tried to set the limits in ggplot2 it interfered with the scale10 transformation. If it can show 0 that will be helpful

I am using data uploaded in phyloseq for the R package microbiome. When creating a community composition plot, I receive the OTU ID as my legend rather than the bacterial name. Is there a way to transform the phyloseq-formatted data, so I can utilize the microbiome package scripts?

Two requests:

1. This could utilize directly the more generic phyloseq functions parse_taxonomy_qiime / _default / _greengenes (check the help of those functions for examples) that know how to perform a similar parsing than is being done in the beginning of this script. Since phyloseq is a central dependency anyway it would make sense to utilize the readily available solution from there. This could simplify the code within format_phyloseq, help to standardize the procedures, and help to avoid bugs by not duplicating solutions.

2. Then it is unclear whether the format_phyloseq function is needed in the first place. Perhaps we could obtain the same by just writing 1-2 lines of code like 'tax_table(x) <- parse_taxonomy_qiime(tax_table(x))' and this could be shown in the tutorial examples. Better not write new functions that just duplicate what is available otherwise. But: I did not examine this properly, so if there is essential functionality on top of the phyloseq parsers then that can justify the format_phyloseq function. Even then, the phyloseq parsers should be used where possible to simplify the code and avoid duplication.

Could you have a look when the time allows ?

SummarizedExperiments is an alternative class structure for multi-assay data objects. This could be very useful in many cases. Start by implementing this for the peerj32 multiassay example data set and exploring the capabilities before considering for the rest of the package.

I cannot use the plot_composition() function because it trows an error:

Error in sample_data(data.frame(x)[i, j, drop = FALSE]) : 
error in evaluating the argument 'object' in selecting a method for function 'sample_data': Error in   `[.data.frame`(data.frame(x), i, j, drop = FALSE) : 
undefined columns selected

Everything works fine if I use the plot_bar() function from the phyloseq package.

> sessionInfo ()
R version 3.2.2 (2015-08-14)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.3 LTS

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=pl_PL.UTF-8            LC_COLLATE=en_US.UTF-8    
[5] LC_MONETARY=pl_PL.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=pl_PL.UTF-8       LC_NAME=C                 
[9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=pl_PL.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] microbiome_0.99.62  RPA_1.26.0          affy_1.48.0         Biobase_2.30.0          BiocGenerics_0.16.0 plyr_1.8.3         
[7] dplyr_0.4.3         ggplot2_1.0.1       tidyr_0.3.1         phyloseq_1.14.0   

We should consider rewriting phyloseq::import_biom to circumwent the problems and to remove the dependency. Same applies to many other phyloseq functions that are currently being imported.

The codecov integration has to be fixed.

(i) Description of purpose: an introduction that lets readers know what the vignette is for; (ii) Resources: short descriptions of the data and r packages utilized in the vignette, to give readers context of the analyses; (iii) Workflow: the bulk of the vignette, normally an analysis starting from data import to inspecting or visualizing the results; (iv) References/Acknowledgements: includes all works cited as well as the authors and contributors of the vignette; (v) R session information: the operating system, version of r and all r packages present when the vignette was rendered, including their versions, installation date and source. The latter is considered a best practice in the r community. Wickham H (2015) R Packages. O'Reilly Media, Sebastopol, California.

Error received when trying to run global() on phyloseq object:
Error in colSums(xcomp[cm, ]) : 'x' must be an array of at least two dimensions

I think I traced down the problem. It seems that colSums() fails on line 25 of core_abundance.R when cm has a length of exactly one.

Here's the problematic section:

if (ncol(xcomp) > 1) {
    colSums(xcomp[cm, ])
} else {
    sum(xcomp[cm, ])
}

I noticed that you're checking the number of columns in xcomp. Should you also make sure that the length of cm is greater than one before calling colSums()?

@antagomir @microsud error in plot_regression.
plot_regression(diversity ~ Age, meta(pseq)) this gives error

Error: All select() inputs must resolve to integer column positions.
The following do not:
*  DV

diversity column already in pseq object in example data atlas1006. but with my data how do i add the diversity column in metadata.

This is how i have added the column in metadata file and it resolved the issue

diversity <- estimate_richness(pseq, measures = "shannon")
metadata <- meta(pseq)
metadata$Shannon <- diversity$Shannon
p <- plot_regression(Shannon ~ Age, metadata)
print(p)

Following many requests it seems we should provide Chao1, ACE etc indicators as part of the harmonized "global" function although they are already in phyloseq as well.

I was thinking to give a more descriptive name for save_tables function .For instance, write_phyloseq_tables or export_phyloseq - any comments ?

Hello everyone, I'm trying to use microbiome tutorial for rda analysis .. the commands should be like :

rda.result <- vegan::rda(t(otu) ~ factor(metadata$time),na.action = na.fail, scale = TRUE)

but if i want to test rda for multiple variables contained within metadata how can I perform this analysis? thanks in advance
An.

.. is named misleadingly, and the presumed shifting functionality has been disabled on the way - investigate and fix.

It was possible to modify the composition_plot() with additional ggplot2 functions.

Now whenever I would add a
+scale_fill_manual()

or any othe ggplot2 function it trows an error:

non-numeric argument to binary operator