adw96 / breakaway Goto Github PK

true_shannon(c(0, 0.5, 0.5)) gives NaN. Need to filter zero abundances before calculating.

When passing through a phyloseq object for convert(input_data) of class:

[1] "otu_table"
attr(,"package")
[1] "phyloseq"

An error is produced:

Error in convert(input_data) :
  Unsupported input type to function `convert`. You passed in an object of class otu_table
Calls: summary -> <Anonymous> -> colSums -> is.data.frame -> convert

Breakaway seems to fail on some of my samples which have less than 6 non-contiguous frequencies (with force=TRUE). I've included the frequency table in attachment. There is no warning/error message.
temp.txt

Code that I used

rich <- breakaway::breakaway(temp, print = FALSE, plot = FALSE, answers = TRUE, force=TRUE)

objective_bayes_negbin needs mvrnorm from MASS, which is not imported.

Either write a mini-MVN random generator, or import MASS.

Quick question regarding betta() output, would it be correct to assume that betta$blups would be considered the "corrected" richness estimations?

i.e.

devtools::install_github("adw96/breakaway")
library(breakaway)
library(phyloseq)
data("GlobalPatterns")
ba <- breakaway(GlobalPatterns)
ss <- summary(ba) 
x_m <- model.matrix( ss$sestimate ~ ss$x_1 + Sample_Type)
bt <- betta(ss$estimate, ss$error, x_m)

bt$blups ==> richness estimates for GlobalPatterns corrected for Sample_Type and x_1
?

Switch sample_names over to the first column:

sample_names estimate error lower upper name model

Possible issue with not skipping frequency counts that have value zero in the weighted linear regression model implementations. TODO: test and confirm

Error in check_format(output_data) :
Input data to check_format is of length 1 or 0.
Calls: summary ... FUN -> breakaway.default -> convert -> check_format

Throw an error? Normalise without warning? Throw warning?

shannon(c(0.2, 0.8))
[1] 0.5004024
shannon(c(0.02, 0.08))
[1] 0.2802988

@bryandmartin Could you come up with some rules for deciding on the default axis limits for plot.alpha_estimates? I don't think it should need to contain the largest upper limit of the CIs. e.g., I think in this example it would be fine to have the lines extending off the top of the page. Thanks!

When using the breakaway plot function of alpha estimates the following ggplot2 error occurs:

Error: Aesthetics must be either length 1 or the same as the data (26): x, y, colour, shape

shannon_e

In test-alpha_estimates.R:

  
  ## this one
  expect_silent({y <- GlobalPatterns %>% 
    subset_samples(SampleType %in% c("Mock")) %>%
    tax_glom("Phylum") %>%
    breakaway})
  
  expect_is(y, "alpha_estimates")
  expect_is(summary(y), "tbl")
  expect_is(plot(y), "ggplot")
}

The y object has all NA estimates, so it breaks the plotting test. Temporarily replaced breakaway with sample_richness in the test.

GlobalPatterns %>% breakaway_nof1

Error in input_data[, 1] : incorrect number of dimensions

Tried to make it work but @bryandmartin is the expert

Hi,
I am studying the ammonia oxidizing archaeal diversity in long term fertilized paddy fields. Roche 454 was performed targeting the archaeal amoA gene and raw sequences were processed using Mothur.
Following that, I would like to calculate the diversity indices using breakaway package.

I have now the OTU table with me and would like to use it for calculation. I tried a few times to convert the table to frequency count table using the command build_frequency_count_tables(the_table), but was unable to do as I am not expert in R language. Below is a brief description how I processed.

I created a .csv of my OTU table with title OTU_table.csv and imported it in R and then used:

build_frequency_count_tables(the_table = "OTU_table")

but it showed

Error in FUN(newX[, i], ...) :
breakaway no longer supports file paths as inputs

On that note, if I can get some help on how can I proceed and calculate the indices it will be really appreciated. I am attaching the OTU table. It contains the OTU numbers in rows and treatments in columns. I have analyzed 5 treatments with three replicates each.
If you can kindly guide how should be the input format and also what will be the steps for calculation
it will be really helpful.

Thanks in advance.

Regards,

Sandipan

OTU_table .xlsx

Hi
I started to play around with breakaway and when reading the Getting started document and in the section Estimating alpha diversity:

BTW, you can use the above to look sample size variability as well. If you have very different numbers of reads across samples, this can introduce additional variability. Let’s look at the distribution of reads in the current dataset.

The code in the documentation is:

ns <- unlist(lapply(frequencytablelist, function(x) sum(x[,2])))

Maybe I am missing something, but this is just summing up the frequencies; shouldn't it be

ns <- unlist(lapply(frequencytablelist, function(x) sum(x[,1]*x[,2])))

to get the number of reads per sample?

Many thanks for developing such great methods!!!
Antonio

Possibly because of phyloseq? Amy to investigate with Bryan

Hello,
Quick question, is it possible to use betta() with multiple covariates and if so, how?
-N

• taxa_are_rows = FALSE
• nof1
• minibreak
• minibreak_nof1
• make taxglom an object so you only have to do it once

Defunct object call in hill_pict(). Fix ASAP

Issue: if f1 = 0, chao_bunge f0-hat = 0, so CI has length zero

Hello Amy,

I was doing the tutorial seen here with breakaway as well as applying it to own data.
In both cases, the breakaway command gives me the following warning

In cutoff_wrap(my_data, requested = cutoff) :
  ignoring requested cutoff; it's too low

Consequently, when I try to plot the breakaway object, I get the follow error
Error: Aesthetics must be either length 1 or the same as the data (11): x, y, colour, shape

I can still run betta() on the breakaway object but plotting directly seems to not work?
Help??

Best,
Nur

Hey folks,
I'm just gearing up to finally give Breakaway a go (woo) and am looking for some clarification. I understand it's important to not discard singletons for proper diversity estimation, as per the q2-breakaway tutorial:

We recommend using a FeatureTable that has been generated from deblur/vsearch or dada2 in R with pool = TRUE to make sure that singletons have not been completely filtered out

I understand denoising pooled reads is more sensitive to rare ASVs and should be used but it comes at quite the computation cost so I'm just wondering if this has been tested and is in fact necessary in the case of paired-end reads. I ask because running DADA2 on paired-end reads still produces singletons since it removes singletons before merging and post-merge singletons are still possible and are not removed as they are considered true ASVs. Would these be adequate for breakaway's estimation models?
edit: What about 'psuedo-pooling' option of DADA2, might be a good compromise as well?

Also worth mentioning, q2-deblur's default parameters also removes rare features so it might be worth clarifying to manually adjust those parameters in future tutorials?

--p-min-reads INTEGER Retain only features appearing at least
min_reads times across all samples in the
resulting feature table. [default: 10]
--p-min-size INTEGER In each sample, discard all features with an
abundance less than min_size. [default: 2]

Thanks!

Hi Amy,

It is great to see the continued development the breakaway package!

The vignettes are extremely helpful and informative. Extending breakaway to better sync with phyloseq and a post/tutorial about when microbial diversity estimation is and isn't possible (i.e. given issues with calling rare ASV/OTUs) would be very informative (saw your recent post on the dada2 repo).

I write to ask if breakaway::betta can be used (or extended) to account for correlated sampling such as that arising from a two arm longitudinal experimental design? We have several ongoing studies where we have an interest in modeling the trajectory of diversity indices/estimators over time and ideally would like to account for the inexhaustive sampling. Our inferences are largely focused on the fixed effects.

Thank you in advance for any thoughts.

Question regarding betta() function. From what I understand, we use that to "test for heterogeneity of total diversity (observed plus unobserved) across multiple sites" and can correct for various covariates that may impact the diversity.
Is it possible to perform residual analysis with the parameter estimates generated by betta()?

breakaway::objective_bayes_negbin(data = my_table[[i]], plot = F, print = F, answers = T)
Error in data.frame(w = w, n = n, NP.est.C = NP.est.n0 + w, tau = tau, :
object 'hist.points' not found

Did I get rid of something that was not, in fact, superfluous? Investigate!

Hi Amy,
when running through the breakaway tutorial, I get stuck at the function "build_frequency_count_tables"

I get the following error message:
Error in build_frequency_count_tables(otu_data) :
could not find function "build_frequency_count_tables"

I have installed breakaway version 4.6.8.
Thank you and best regards,
Maren

Needs to be more stable

Hello Amy,
I am trying breakaway. I am trying to compare fungal and bacterial observed vs. rarefied and also modeled (see McMurdie and Holmes, 2014 PlosONE) richness, in a set of soil samples. I wasn't able to find richness estimation models for single samples in breakaway I noticed the git-hub page is still a working progress that is why I am asking you directly.
Thanks, Gian

Howdy! Leaving the process I did to get NaNs for breakaway's error and breakaway_nof1's estimate+error in a little script. Thanks so much for taking a look.

breakaway_na_ex.txt

Hi,
I have been playing around with the Breakaway package and I think I am slowly geting the grasp of it. However, I when it comes to the output of the beta() model, I can't get my head around the results$ssq_u value.
I assume is testing the heterogeinity of the groups; can I get a significant value from it? I am not sure how to interprete this.

Thank You,
Paula

Use alpha_estimates class in betta:

Currently to run:

bt <- betta(summary(ba)$estimate,
summary(ba)$error,
DivNet::make_design_matrix(water, "SampleType"))
bt$table

ideally could run
betta(ba, water, "SampleType") to get the same output

Opened in adw96/stamps2018#23

a2 <- apples[-(1:20),]
sum(a2[,2])
chao_bunge(a2)
chao_bunge(a2, cutoff = 30)
chao_bunge(a2, cutoff = 300)

All of these have NA standard errors

It would be great to have the names of the samples embedded in the alpha_estimates object if it is created from a phyloseq object. TODO(Amy)

Hi, I noticed that in wlrm, setting cutoff for maximum frequency count to be a number greater than number of rows in dataset throws an error. Thanks!

Hey awesome breakaway crew,

I m running breakaway on my own data and I ve come an error I can't really figure out.
running breakaway(ps2) returns

Error in alpha_estimate(estimate = diversity, error = diversity_se, estimand = "richness", :
object 'f0' not found
this doesn't happen if I use a preloaded dataset from phyloseq like Gobalpatterns .

this is my Phyloseq object
ps2

phyloseq-class experiment-level object
otu_table() OTU Table: [ 3013 taxa and 292 samples ]
sample_data() Sample Data: [ 292 samples by 43 sample variables ]
tax_table() Taxonomy Table: [ 3013 taxa by 8 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 3013 tips and 2977 internal nodes ]

thanks for the help and greetings from Germany :)
Rima

..and have them output objects of class alpha_estimates

yy <- chao_bunge(GlobalPatterns)
zz <- wlrm_transformed(GlobalPatterns)
aa <- wlrm_untransformed(GlobalPatterns)
bb <- chao1(GlobalPatterns)
cc <- chao1_bc(GlobalPatterns)
dd <- good_turing(GlobalPatterns)

and then
plot(yy)
plot(zz)
...

should all work.

Assigned to @bryandmartin

TODO (Amy)

mvtnorm needs to be removed from the DESCRIPTION file

alphas <- breakaway(GlobalPatterns %>% 
                        subset_samples(SampleType == "Soil"))

throws possibly serious warnings; investigate

to facilitate easy install

breakaway_nof1 should always produce some estimate, and doesn't if no models converge. Perhaps have breakaway()-like pipeline that cycles through WLRM_ut and WLRM_t?

Running Atacama data and the following error is produced. Using breakaway_4.6.6

library(phyloseq)
library(breakaway)
library(RCurl)
x <- getURL("https://raw.githubusercontent.com/paulinetrinh/data/master/otu_table_atacama.txt")
a_table <- read.table(text =x,
                    skip = 0, 
                    header = F, 
                    row.names = NULL)[,-1]
colnames(a_table) <- colnames(read.csv(text = x, nrows=1, skip=1, sep = "\t"))[-1]
dim(a_table) #Should be 401 rows X 61 columns
a_ps <- phyloseq(otu_table(a_table, taxa_are_rows = TRUE)) 
a_breakaway <- breakaway(a_ps)

It would be nice to have a link to the documentation for breakaway from the main GitHub Readme page. If it's there and I missed it, my apologies.

Utilize the structures data(apples) in order to avoid overwriting existing variables.

My apologies to anyone whose apples have been eaten by breakaway!

set.seed(1)
y <- rnbinom(10000, size=0.01, prob=0.01)
y %>% breakaway::make_frequency_count_table %>%
breakaway::chao_bunge

in docs/

y <- cbind(c(1,3,5,6,11,13,17,21), c(13,2,2,3,2,1,1,2))
y
chao_bunge(y)