John Bunge jab18@cornell Department of Statistical Science Cornell University

1. EstimatingMicrobial Diversity John Bunge jab18@cornell.edu Department of Statistical ScienceCornell University

2. Thanks to: Amy Willis Fiona Walsh David Mark Welch Colleagues too numerous to mention Bunge, J., Willis, A. and Walsh, F. (2013) Estimating the number of species in microbial diversity studies. Ann. Rev. of Statist. and its Appl. v.1. Forthcoming.

3. Statisticians

4. Bioinformaticists

5. Statistics is not a collection of formulae, nor computer programs, but a conceptual framework, an intellectual stance, a point of view, a theory of knowledge Fundamental idea:distinction between sample and population Classical or frequentist statistics is fundamentally dualistic

6. Plato’s Republic, VII,7 Behold! human beings living in an underground den, which has a mouth open towards the light and reaching all along the den; here they have been from their childhood […]Above and behind them a fire is blazing at a distance, […] you will see, if you look, a low wall built along the way, like the screen which marionette players have in front of them, over which they show the puppets.  […] They see only their own shadows, or the shadows of one another, which the fire throws on the opposite wall of the cave […] To them, I said, the truth would be literally nothing but the shadows of the images.

7. Old Testament Ecclesiastes 1:15 What is crooked cannot be straightened; what is lacking cannot be counted. New Testament Corinthians 13:12 For now we see through a glass, darkly, but then face to face: now I know in part; but then shall I know even as also I am known.

8. The knowledge problem in microbiome studies • Metagenomics is the study of metagenomes, genetic material recovered directly from environmental samples. • Wikipedia • DNA extraction bias notwithstanding, metagenomics is the most unrestricted and comprehensive approach. Our ability to interpret these data is always improving, and we stand on a precipice of unprecedented discovery […] Microbes are not the only group to benefit from these surveys; viruses exist at 10 times the abundance of microbes […]. - Gilbert, 2011 • BUT: Metagenomicsurveys recover only a small fraction of the extant diversity. Nonetheless, many methods treat the observed sample as the population.

9. Machines The fundamental idea of statistics: Distinction betweenPopulation (or universe) and Sample (or data)

10. The sample is a subset of the population Population Universe Reality State of nature Truth parameters Sample Finite, random noise error perturbation shock statistics Statistical inference: Extract maximum information from sample in order to draw conclusions about population Inductive not deductive

11. Question: In a microbial diversity study, What is the population? • Collect 1L seawater @ 500m depth in ocean • From 1L, remove 5ml & exhaustively sequence microbial DNA • Cluster sequences into OTUs • From OTUs, calculate frequency countdata • Compute estimateof total species richness • Question: Richness of what population? Original 1L of water? Surrounding environment? Entire pelagic microbiome? DefinitionThe population is what would be observed if the operative sampling and analysis protocols were carried out to infinite effort.

12. How do we statistically estimate total microbial taxonomic richness?

13. Cluster sequences at some % “identity,” typically 97% {clusters} = {OTUs} OTU = “operational taxonomic unit”

14. Statistical problem: • Estimate total population diversity – number of species, classes, taxa, OTUs – based on frequency count data • Data = # of units observed exactly once in sample (singletons); • # observed exactly twice (doubletons); # observed exactly three times; … .

15. Frequency count data example Microbial ecology Fiona Walsh et al. • Data from soil in apple orchards • Use of antibiotics on bacterial populations in soil ecosystems • Singletons ≈ 2x doubletons – may be 10x! • Goal is to estimate taxonomic richness of community • Change with respect to intervention/covariates/metadata Walsh F, Owens S, Duffy B, Smith DP, Frey JE. 2013. Streptomycin use in apple orchards did not alter the soil bacterial communities

16. Issues: • High diversity • Typical of microbial data • Singletons ~ 2x doubletons • Data acquisition / bioinformatic issues • Spurious singletons? • Correct at what stage? Statistical approach?

17. Statistical inference from frequency count data • Standard Model • C classes/taxa/species in population. Each species independently contributes Poisson-distributed # of representatives to the sample. • Counts ~ zero-truncated mixed Poisson. sample

18. The mixed-Poisson model • Species (taxon) i contributes a Poisson-distributed number Xi of replicates to the sample – i.e., taxon i appears in the sample Xi times. • Units appear independently in the sample • Fundamental problem: heterogeneity, i.e., unequal Poisson means λi • Standard approach: model λi‘s as i.i.d. replicates from some mixing distribution F • Frequency counts fi are then marginally i.i.d. F-mixed Poisson random variables • Zero-truncated since zero counts Xiare unobservable

19. The mixed-Poisson model cont’d • Mixing distribution F, i.e., distribution of sampling intensities λ, is also called species abundance distribution • Probably a misnomer • Mathematical treatment (marginalization) implies that each species contribution to the sample is independent and identically distributed • Both assumptions are certainly wrong • How to account for dependent or differently distributed species counts? Not in standard model.

20. Mixing distributions F • Parametric, low-dimensional parameter vector • None ≡ point mass at λ ≡ all equal species sizes • Gamma (Fisher, 1943) • Lognormal • Inverse Gaussian, generalized inverse Gaussian (Sichel) • Pareto • Log-t • Stable • Finite mixture of exponentials - semiparametric

21. Richness estimation under the Poisson model • Diversity estimate is then • where PF(0) = F-mixed Poisson probability of 0: • is the Horvitz-Thompson estimator (HTE) and is uniformly minimum variance unbiased (UMVU). • Require empirical version of , i.e., require estimate of PF(0) (frequentist version).

22. Richness estimation under the Poisson model, cont’d • Require empirical version of HTE • Estimate θ by ML, using zero-truncated F-mixed Poisson, conditional on # of observed taxa. Final estimator: • SE via Fisher information • CI via (approximation to) profile likelihood

23. CatchAll softwarewww.northeastern.edu/catchall • or: STAMPS! • Developed under NSF grant DEB – 0816638 by JB/LW/SC, in C# & C • Implements • finite mixtures of 0 – 4 exponential components (F) • weighted linear regression procedure • all Chao-type nonparametric procedures • model evaluation/GOF/selection/outlier assessment • Produces estimates, SEs, & CIs • Fast, efficient, platform-independent • Excel graphics (VBA) package • Summary or copious output (text files) Bunge J, Woodard L, Böhning D, Foster JA, Connolly S, Allen HK. 2012b. Estimating population diversity with CatchAll. Bioinformatics 28:1045--47

24. Partial CatchAll summary output for apple orchard data

25. CatchAll fitted models for apple orchard data Τ = 184

26. Data-analytic considerations • Problem of right cutoff point τ • Typically no parametric model will fit complete frequency count dataset • Too many right outliers – highly abundant taxa in sample – with large gaps between counts • Nonparametric methods do even worse with outliers, diverging to ∞ as outliers are included in data • Data-analytic solution: remove large frequency counts for frequencies > some cutoff τ • Chao1: τ = 2 • Chao-type coverage-based nonparametric methods: τ = 10 (arbitrary) • Parametric mixture models: τ selected by goodness-of-fit algorithm • Weighted linear regression model: selected by goodness-of-fit • Further problem: model selection and outlier deletion confounded • Computational solution: compute all methods at every τ • Requires optimized code • Use double selection algorithm to select “best of the best” • Introduces simultaneous inference problem: large number of simultaneous GOF tests. Little theory exists to correct for this.

27. Statistical analysis of standard model: The bigger picture

28. Statistical analysis of standard model – Chao-type nonparametrics • Coverage-based approaches • Coverage = proportion of population represented in sample • Random variable not parameter • Can interpret 1 – PF(0) as surrogate for coverage • Turing’s estimate of PF(0): where n = # of individual units in sample • Good-Turing estimateof diversity: • Chao’s abundance-based coverage estimators (ACE): • Good-Turing + adjustment for heterogeneity Chao, A. & J. Bunge. 2002. Estimating the number of species in a stochastic abundance model. Biometrics 58: 531–539

29. Coverage-based estimators diverge to infinity as large frequency counts are included Hence coverage-based estimators require τ ≤ 10

30. Statistical analysis of standard model: general nonparametrics • Nonparametric maximum likelihood estimation • Leave species abundance distribution F unspecified, i.e., F varies across all possible distributions • Mathematical implications: F is actually non-identifiable • Nevertheless NPMLE is possible in principle. • Computational issues: difficult numerical search, highly complex error estimation. • Software CAMCR Böhning D, Kuhnert R. 2009. CAMCR: Computer-Assisted Mixture model analysis for Capture-Recapture count data. AStA Adv. Stat. Anal. 93:61--71

31. The Bayesian paradigm • Rev. Thomas Bayes • Bayesian statistics: Probabilistic & statistical statements concern degrees of belief • Usually parametric: statements concern values of parameters, e.g., species richness. Nonparametric Bayes is possible but complex. • Procedure: • Investigator first declares existing belief about population value: this is prior distribution • Collect sample data • Update prior, based on data, to obtain posterior, i.e., final state of knowledge or belief about population.

32. The Bayesian paradigm cont’d Bayes’ Theorem: Posterior distribution: Bayesian computation is now fairly well established

33. Bayesian estimation of taxonomic richness based on the standard model • Species abundance distribution F is parametric: F depends on a small number of parameters (typically 2-3), called  • Parameter of interest is total richness C • Procedure: • Establish prior distributions for  and C • Likelihood function is known (based on mixed-Poisson) • Run Bayesian machinery • Obtain posterior distribution, estimate, “credible interval,” etc. • Quince et al. quasi-noninformative priors; Barger et al. formal objective priors. Active research area in statistics. Quince C, Curtis TP, Sloan WT. 2008. The rational exploration of microbial diversity. ISME J. 2:997—1006; Barger K, Bunge J. 2011. Objective Bayesian estimation for the number of species. J. Bayesian Analysis 5:765--86

34. A New Hope • Is it possible to estimate taxonomic richness without • a species abundance distribution • independent species contributions to the sample • identically distributed species contributions to the sample • ? • Yes, using ratios of frequency counts.

35. breakaway: Estimating taxonomic richness based onratios of frequency counts Idea: ratios are ~ linear Project line downward to obtain f0 = # of unobserved species

36. breakaway: Estimating taxonomic richness based onratios of frequency counts, cont’d • Some issues: • Straight-line fit may go negative! • Can be fixed by ad hoc log-transformation (Rocchetti et al.) • Broad generalization: represent ratio of frequency counts as ratio of polynomials • Deep probabilistic justification; corrects negativity Rocchetti I, Bunge J, Böhning D. 2011. Population size estimation based upon ratios of recapture probabilities. Ann. Appl. Stat. 5:1512—33; Willis A. and Bunge J. (2013) in prep.

37. breakaway: Estimating taxonomic richness based onratios of frequency counts, cont’d ################## Smoothed weights ################## The best estimate of total diversity is 1800 with std error 256 The model employed was model_1_1 The function selected was f_{x+1}/f_{x} ~ (beta0+beta1*(x-xbar))/(1+alpha1*(x-xbar)) Coef estimates Coef std errors beta0 1.11078693 0.13241518 beta1 0.05383757 0.02916098 alpha1 0.03002143 0.03840271

38. breakaway: Estimating taxonomic richness based onratios of frequency counts, cont’d • Nonlinear regression • Heteroscedastic (changing variance) • Autocorrelated: f2/f1 is correlated with f3/f2, etc. • Collinear: parameter estimates of α’s and β’s highly correlated unless corrected • Multiple significant numerical challenges • Statistical questions • Model selection – degree of numerator and denominator polynomials • Error estimation • Underlying probability theory: what do these models imply, and what are they implied by?

39. Noise and unreliable low frequency counts Next generation sequencing technology […] has revolutionised the study of microbial diversity as it is now possible to sequence a substantial fraction of the 16S rRNA genes in a community. However, […] because of the large read numbers and the lack of consensus sequences it is vital to distinguish noise from true sequence diversity in this data. Otherwise this leads to inflated estimates of the number of types or operational taxonomic units (OTUs) present. - Quince et al. (2011)

40. I. Fix the data at the source! • Example: PyroNoise and AmpliconNoise • - aim at “separately removing 454 sequencing errors and PCR single base errors.” (Quince 2011) • Direct, non-statistical approach Methods to address unreliable low frequency counts

41. Methods to address unreliable low frequency counts

42. III. Deleting the high-diversity component of a mixture model Methods to address unreliable low frequency counts Bunge J, Böhning D, Allen H, Foster JA. 2012a. Estimating population diversity with unreliable low frequency counts. In Biocomputing 2012: Proceedings of the Pacific Symposium, pp. 203--12. Hackensack, NJ: World Sci. Publ

43. IV. Bayesian approaches • Informative or subjective: investigator specifies non-trivial downweighting or rapidly decreasing prior for higher diversity values • Specific choice of prior? Methods to address unreliable low frequency counts

44. Numerical results from viral phage data:Lower bounds and component deletion

45. Some notes on β-diversity Crucial to distinguish between Statistical inference procedures that (attempt to) account for unobserved as well as observed diversity Procedures (computational, graphical, or qualitative) that treat the observed sample as the population. UniFrac, “ordination” methods, co-inertia. Only the former considered here. Estimation of population parameters, possible hypothesis testing.

46. Statistical inference for comparing taxonomic diversity across populations Simplest version: Estimate richness in each population, with associated standard errors and confidence intervals, & compare (e.g., do CI’s overlap?) Can be done with existing methods: parametric, nonparametric, Bayesian, etc. Exactly ONE known inferential procedure. Lower bound for # of shared taxa: (D12= observed # of shared species, fjk= # of species observed j times in sample 1 and k times in sample 2, aand b = constants) Pan HY, Chao A, Foissner W. 2009. A nonparametric lower bound for the number of species shared by multiple communities. J. Agric. Biol. Environ. Stat. 14:452--68

47. Statistical inference for β-diversity:other scenarios Inference for the Jaccard index, accounting for unobserved species (Chao et al.) Inference for “the probability of a draw from one distribution not being observed in k draws from another distribution.” (Hampton et al.) Statistical work in this area not extensive – very fertile area for research. Chao A, Chazdon RL, Colwell RK, Shen T-J. 2006. Abundance-based similarity indices and their estimation when there are unseen species in samples. Biometrics 62:361—71; Hampton J, Lladser ME. 2012. Estimation of distribution overlap of urn models. PLoS ONE 7:e42368

48. NEVER throw away data when doingstatistical inference “Not even wrong” – Richard Feynman

49. There is no post hoc statistical fix for • Ill-posed research problem • Vaguely defined population • Statistical model not appropriate for • population description • sample generation process • Model must compromise between detailed phenomenological description and parsimony • “To what extent can we idealize the properties of the system and still obtain satisfactory results? The answer to this question can only be given in the end by experiment. Only the comparison of the answers provided by analysis of our model with the results of the experiment will enable us to judge whether the idealization is legitimate.”Andronov (1937) Theory of Oscillators.

50. On the sociology of science • Fact: Universities have statistics departments! • Cornell: www.stat.cornell.edu • At least 131 university stat dept’s in U.S. – random sample of 10: • • University of California, Berkeley, Division of Biostatistics • Princeton University, Program in Statistics and Operations Research • Bowling Green State University, Department of Applied Statistics and Operations Research • University of Illinois, Urbana-Champaign, Department of Statistics • University of South Carolina, Department of Statistics • Columbia School of Public Health, Division of Biostatistics • Medical College of Georgia, Office of Biostatistics and Bioinformatics • Duke University, Institute of Statistics and Decision Sciences • Yale University Department of Statistics • University of Michigan, Department of Biostatistics • Collaboration extremely valuable in both directions (even though academic incentive structure may not immediately reward it) • Be persistent: “Fall down seven times, get up eight”