---
abstract: |
    We propose two unsupervised clustering methods that are designed for human microbiome data.  Existing clustering approaches do not fully address the challenges of microbiome data, which are typically structured as counts with a fixed sum constraint.  In addition to accounting for this structure, we recognize that high-dimensional microbiome datasets often contain uninformative features, or “noise" operational taxonomic units (OTUs), that hinder successful clustering.  To address this challenge, we select features which are useful in differentiating groups during the clustering process.  By taking a Bayesian modeling approach, we are able to learn the number of clusters from the data, rather than fixing it upfront.  We first describe a basic version of the model using Dirichlet multinomial distributions as mixture components which does not require any additional information on the OTUs.  When phylogenetic or taxonomic information is available, however, we rely on Dirichlet tree multinomial distributions, which capture the tree-based topological structure of microbiome data.  We test the performance of our methods through simulation, and illustrate their application first to gut microbiome data of children from different regions of the world, and then to a clinical study exploring differences in the microbiome between long and short term pancreatic cancer survivors.  Our results demonstrate that the proposed methods have performance advantages over commonly used unsupervised clustering algorithms and the additional scientific benefit of identifying informative features. 

    **Contact:** cbpeterson@mdanderson.org

    **Keywords:** Unsupervised clustering, Bayesian mixture model, Mixture of finite mixtures model, Feature selection, Dirichlet tree Multinomial
author:
- 'Yushu Shi [^1]'
- 'Liangliang Zhang [^2]'
- 'Kim-Anh Do'
- 'Robert Jenq [^3]'
- 'Christine Peterson [^4]'
title: Bayesian Approaches for Flexible and Informative Clustering of Microbiome Data
---

Introduction
============

The microbiome plays an important role in human health.  Previous studies have shown that human stool microbiome samples naturally form clusters that are associated with dietary and geographic factors [@MartSteg15].  In addition, clinical studies have shed light on the influence of microbiome on the natural history of diseases and the host immune response to treatments [@GopaSpen18].  In the past decade, next generation sequencing methodologies have enabled researchers to cheaply and comprehensively analyze microbial communities [@GilbBlas18; @KnigCall17].  Downstream analysis of microbiome data requires grouping similar sequences, and the highest resolution bacterial sequence is referred as an operational taxonomic unit (OTU).  Though many high-dimensional data analysis methodologies are applicable to microbiome data [@LoveWolf14], accounting for the relatedness of OTUs remains challenging.  Given the representative sequence of each OTU, researchers can construct a tree structure by combining similar OTUs [@PricDeha10], which can approximate evolutionary relatedness [@DArg18; @GoodDi14].  Alternatively, researchers map these reference sequences based on their similarity to known species, which can be organized into a classical taxonomic tree.  This traditional taxonomy is based on microbiological features but more recently increasingly takes genetic similarity into account. 

Both machine-learning methods and model-based methods have been proposed for clustering of microbiome samples.  Most machine-learning methods require first determining pair-wise distances between samples using metrics such as Bray-Curtis dissimilarity [@BrayCurt57], unweighted UniFrac distance [@LozuKnig05], or weighted UniFrac distance [@LozuHama07], and then applying algorithms, such as PAM [@KaufRous08], K-means [@MacQ67] or hierarchical clustering, with the number of clusters pre-specified by the user.  When the number of clusters is unspecified, one can calculate the silhouette width [@Rous87] for a range of number of clusters.  A higher silhouette width implies more appropriateness of the choice.  Other commonly-used methods for determining the number of clusters, such as the gap statistic, can only be applied to Euclidean distances, and are therefore not suitable for microbiome data. 

As an alternative to machine-learning methods, a model-based unsupervised clustering method was proposed by @HolmHarr12 using Dirichlet multinomial mixtures (DMM).  They proposed determining the number of clusters by calculating model evidence via Laplace approximation over a range of possible values and choosing the maximum.  The DMM method, which does not perform feature selection, can be applied to a few hundred variables, and may therefore be used to analyze species or genus-level data.  It does not scale, however, to the thousands of features often encountered when analyzing OTUs. 

Existing machine-learning and model-based clustering methods have the limitation that the number of clusters needs to be either pre-fixed or chosen in a post-hoc fashion.  When the number of clusters is not prespecified, one can treat the data as arising from an infinite mixture of a distribution and use Bayesian nonparametric methods, such as the Dirichlet process mixture (DPM) model.  However, it has been shown that the number of clusters will grow with the number of observations in a DPM model [@MillHarr14].  Alternatively, one can treat the data as from a finite mixture of a given distribution and use methods such as reversible jump Markov chain Monte Carlo [@RichGree97; @TadeSha05] to learn the number of clusters from the data.  @MillHarr18 proved that the mixture of finite mixtures model (MFM) can consistently estimate the number of clusters, while the number of clusters in Dirichlet process mixtures will increase with sample size.  However, when there are only a small number of observations, the two methods generate similar results [@MillHarr18]. 

In this paper, we use the mixture of finite mixtures model, which puts a prior on the number of clusters.  To efficiently handle the high dimensionality of the microbiome data, we exploit the conjugacy between the Dirichlet distribution and the multinomial distribution in the choice of mixture component distributions.  When there is no phylogenetic tree describing the relationships among OTUs, we use a Dirichlet multinomial distribution.  When the phylogenetic tree of the OTUs is provided, we use the Dirichlet tree multinomial distribution.  The Dirichlet tree multinomial distribution [@Denn91] can be thought of as an extension of the Dirichlet multinomial to account for phylogenetic tree topological structure.  It does not have a variance constraint, and can capture the dependency between OTUs.  This distribution has been applied in information retrieval and human machine interface [@HaffTeh09; @LiuChen16], and recently used in microbiome regression models [@WangZhao17; @TangNico17; @TangLi18]. 

We hypothesize that microbiome datasets can often contain “noise" OTUs that mask signal from informative OTUs and hinder successful clustering.  We therefore select informative OTUs or phylogenetic tree nodes by updating a latent indicator during the clustering processes.  Similar discriminatory feature selection techniques have been developed for unsupervised text clustering [@HuangYu13; @LiGuin17].  When the mixture kernel is a Dirichlet tree multinomial, the model selects a parsimonious set of discriminatory taxa, enabling biological insight into the cluster assignments.  We start with the @DeFiCava10 stool sample dataset, which consists of two natural clusters corresponding to 14 children with a vegetarian diet from Burkina Faso (Africa) and 15 Italian children with a Mediterranean diet (Italy).  The samples form the two clusters are relatively distinct due to great geographical separation and dietary differences.  This dataset serves as a verification of our method and the basis for simulation studies.  We then apply our method to a pancreatic adenocarcinoma (PDAC) dataset [@RiquZhan19], where the natural grouping of samples is less obvious.  The dataset consists of 68 pancreatic patients from two hospitals, where each patient contributes one tumor sample.  Among them, 36 patients survived longer than 5 years after the surgery (LTS: long-time survivors), while the rest died within five years after the surgery (STS: short-time survivors). 

The paper is structured as follows.  In Section \[sec:method\], we present the formulation of feature (OTU or phylogenetic tree node, depending on the data given) selection and the clustering method.  In Section \[sec:mcmc\], we describe the implementation of the proposed methods.  Section \[sec:application\] conducts simulations and implements our method on two application datasets with naturally existing clusters.  Section \[sec:discussion\] concludes the paper with a brief discussion. 

Formulation of the proposed method {#sec:method}
==================================

In this section, we first describe the likelihood of the data without feature selection and clustering.  Then, we introduce a latent indicator representing feature inclusion.  Finally, we discuss @MillHarr18’s mixture of finite mixtures (MFM) model in detail and demonstrate how this model achieves clustering without pre-specifying the number of clusters. 

Likelihood without feature selection and clustering
---------------------------------------------------

The simplest form of a microbiome dataset involves an $N \times d$ matrix, $\mathbf{Y}$, where $N$ is the number of observations, and $d$ is the total number of features.  The entry $y_{i,j}$ represents the number of genetic sequences observed for the $j$th OTU in the $i$th observation.  In a multinomial model, the $i$th observation is modeled as $\boldsymbol{y}_i \sim \mathrm{Multinomial}(q_1, q_2, \dots, q_d).$ Integrating out parameters under the conjugate prior $\mathrm{Dirichlet}(\alpha,\alpha,\dots,\alpha)$ yields a Dirichlet multinomial distribution:

$$P(\mathbf{y}_i|\alpha)=\frac{y_{i,.}!}{y_{i,1}!y_{i,2}!\dots y_{i,d}!}\frac{\Gamma(d\alpha)}{\Gamma(\alpha)^d}\frac{(y_{i,1}+\alpha)!(y_{i,2}+\alpha)!\dots (y_{i,d}+\alpha)!}{(y_{i,.}+d\alpha)!}.$$

The Dirichlet-multinomial distribution is frequently used to model microbiome data, as it better captures overdispersion than a simple multinomial likelihood [@HolmHarr12; @Chen2013].  Noting that related species within a tree may often have similar functions, we want to better model the structure present in the data by accounting for taxonomic or phylogenetic trees.  An extension of the Dirichlet distribution, the Dirichlet tree distribution, can help us achieve this goal while maintaining conjugacy to the multinomial distribution.  Figure \[Dirichlettree\] gives an example of the Dirichlet tree distribution, where the probability of choosing a leaf is the product of the branch probabilities leading to that leaf. 

![ An example of the Dirichlet tree distribution.[]{data-label="Dirichlettree"}](plots/capture2){width="0.5\linewidth"}

To show the conjugacy with a Dirichlet tree distribution, a multinomial sample can be represented as the outcome of a finite stochastic process.  Given the phylogenetic tree structure $T$ and the transilient probabilities $B$ between nodes and leaves, the probability of one genetic sequence $x$ can be written as: $$P(x|B,T)=\prod_{j \in J} \prod_{k \in K_j} b_{jk}^{\delta_{jk}(x)},$$ where $J$ is the set of parent nodes, $K_j$ is the set of nodes or leaves descending from node $j$, and $\delta_{jk}(x)$ is the indicator of whether the genetic sequence passes through the branch linking node $j$ and node $k$. Similarly, the probability of the $i$th sample is: $$P(X_i|B,T)=\prod_{j \in J}\frac{n_{j,.}(X_i)!}{n_{j,1}(X_i)!  n_{j,2}(X_i)!  \dots n_{j,|K_j|}(X_i)!}\prod_{k \in K_j} b_{jk}^{n_{jk}(X_i)},$$ where $n_{j,.}(X_i)$ is the sum of the counts of all the genetic sequences descending from node $j$ for the $i$th observation, and $n_{j,k}(X_i)$ is the number of genetic sequences descending from node $j$ to node $k$ for the $i$th observation, $|K_j|$ is the number of children of node $j$.

A Dirichlet tree distribution can be expressed as a product of Dirichlet densities, $P(B)=\prod_{j \in J} p(\boldsymbol{b}_j)$. If each $\boldsymbol{b}_j$ is given a Dirichlet prior, $P(\boldsymbol{b}_j|\alpha)=\mathrm{Dirichlet}(\alpha,\alpha,\dots,\alpha)$, one can obtain a Dirichlet tree multinomial distribution for observation $i$ by integrating out $B$:

$$\begin{aligned}
P(X_i|T)=&\prod_{j \in J}\left\{\frac{n_{j,.}(X_i)!}{n_{j,1}(X_i)!n_{j,2}(X_i)!\dots n_{j,|K_j|}(X_i)!}  \cdot \right.\\
\nonumber&\left.\prod_{k \in K_j} \frac{\Gamma(|K_j| \alpha)}{[\Gamma(\alpha)]^{|K_j|}}\frac{(\prod_{k \in K_j}\Gamma(\alpha+n_{jk}(X_i)))}{\Gamma(|K_j| \alpha+\sum_{k \in K_j} n_{jk}(X_i))}\right\}.\end{aligned}$$

Selection of OTUs and phylogenetic tree nodes
---------------------------------------------

In the current work, clustering of samples is achieved through identifying a parsimonious set of parameters which differ across groups, while the rest of the parameters are not cluster-specific.  To obtain a sparse model that selects informative OTUs or phylogenetic tree nodes, we introduce a binary vector $\boldsymbol{\gamma}$, whose entries indicate whether the corresponding features are useful in discriminating between groups of samples, where

$$\gamma_j=\left\{
\begin{array}{ll}
1 & \text{if }j\text{ is an informative feature}\\
0 & \text{otherwise.} \\ 
\end{array}
\right.$$

For simplicity, we use a Bernoulli prior for $\boldsymbol{\gamma}$: $p(\boldsymbol{\gamma}) \propto \prod_{j=1}^d w^{\gamma_j}(1-w)^{1-\gamma_j}$ [@GeorMcCu93; @MadiYork95], where $d$ is the total number of features, and the prior probability of a feature being informative is $w$.

### OTU selection in the mixture of Dirichlet-multinomials model

In the proposed parsimonious model, the parameter set of the $i$th observation can be divided into a common part shared across all observations, $p_1, p_2, ..., p_{d-d_{\boldsymbol{\gamma}}}$, and the part specific to the cluster the $i$th observation belongs to, $w_{c_i}$, $q_{{c_i},1}, q_{{c_i},2}, ..., q_{{c_i},d_{\boldsymbol{\gamma}}}$, where $c_i$ indicates cluster membership, and $d_{\boldsymbol{\gamma}}$ denotes the number of informative features.  Each observation is modeled as a multinomial variable with parameters $\mathbf{y}_i \sim \mathrm{Multinomial}(w_{c_i} p_1,w_{c_i} p_2,$ $\dots,$ $w_{c_i} p_{d-d_{\boldsymbol{\gamma}}} ,(1 - w_{c_i}) q_{{c_i},1}, (1 - w_{c_i}) q_{{c_i},2},\dots,$ $(1 - w_{c_i}) q_{{c_i},d_{\boldsymbol{\gamma}}})$, where $w_{c_i}$ controls the proportion of sequences belonging to the “noise" OTUs in cluster ${c_i}$. The length of the vectors $\boldsymbol{p}=(p_1, p_2, ..., p_{d-d_{\boldsymbol{\gamma}}})$ and $\boldsymbol{q}_{c_i}=(q_{{c_i},1}, q_{{c_i},2}, ..., q_{{c_i},d_{\boldsymbol{\gamma}}})$ will change with the number of OTUs selected as informative, but both $\boldsymbol{p}$ and $\boldsymbol{q}_{c_i}$ sum to $1$, thus guaranteeing the mean vector of the multinomial distribution will sum to $1$ as well.  Finally, the likelihood for observation $i$ can be written as:

$$\begin{aligned}
P(\mathbf{y}_i|c_i,\boldsymbol{p},\boldsymbol{q}_{c_i},w_{c_i},\boldsymbol{\gamma})=&\frac{y_{i}!}{y_{n,i,1}!\dots y_{n,i,d-d_{\boldsymbol{\gamma}}}!  y_{e,i,1}!\dots y_{e,i,d_{\boldsymbol{\gamma}}}!}\\
\nonumber& \times(w_{c_i} p_1)^{y_{n,i,1}}\dots (w_{c_i} p_{d-d_{\boldsymbol{\gamma}}})^{y_{n,i,d-d_{\boldsymbol{\gamma}}}}\\
\nonumber& \times [(1-w_{c_i}) q_{{c_i},1}]^{y_{e,i,1}}\dots [(1-w_{c_i}) q_{{c_i},d_{\boldsymbol{\gamma}}}]^{y_{e,i,d_{\boldsymbol{\gamma}}}}.\end{aligned}$$

For observation $i$, $y_{n,i,j}$ is the number of sequences of “noise" OTU $j$, $y_{e,i,l}$ is the number of sequences of the informative OTU $l$, and $y_{i,.}=y_{n,i,.}+y_{e,i,.}$ is the total number of sequences.  The likelihood for all the observations is $P(\mathbf{Y}|\boldsymbol{p},\boldsymbol{q},\mathbf{w},\boldsymbol{\gamma},\mathbf{c}) = \prod_{c\in C}\prod_{c_i=c} P(\mathbf{y}_i|c_i,\boldsymbol{p},\boldsymbol{q}_{c_i},w_{c_i},\boldsymbol{\gamma})$, where $C$ denotes the set of distinct cluster indicators. 

In order to obtain a more tractable posterior distribution, we use the simplest conjugate priors by setting

$$\begin{aligned}
   p_1, p_2, \dots, p_{d-d_{\boldsymbol{\gamma}}} &\sim \mathrm{Dirichlet}( \alpha , \alpha , \dots, \alpha )\\  
 \nonumber q_{c,1}, q_{c,2}, \dots, q_{c,d_{\boldsymbol{\gamma}}} &\sim \mathrm{Dirichlet}( \alpha , \alpha , \dots , \alpha)\\
 \nonumber w_c &\sim \mathrm{Beta}( \beta_1 , \beta_2). 
 \end{aligned}$$

For simplicity and objectivity, we assume all the parameters in the Dirichlet distribution are equal.  Since our primary interest is in learning the cluster assignments and identifying discriminating features, we integrate out the remaining parameters to speed up computation.  The resulting marginal likelihood is: $$\begin{aligned}
 P(\mathbf{Y}|\boldsymbol{\gamma},\mathbf{c})=&\prod_{i=1}^{N}\frac{y_{i}!}{y_{n,i,1}!\dots y_{n,i,d-d_{\boldsymbol{\gamma}}}!  y_{e,i,1}!  \dots y_{e,i,d_{\boldsymbol{\gamma}}}}  \frac{\Gamma((d-d_{\boldsymbol{\gamma}})\alpha)}{[\Gamma(\alpha)]^{d-d_{\boldsymbol{\gamma}}}} \\
 \nonumber&\times \frac{\Gamma(\sum_{i=1}^N y_{n,i,1}+\alpha)\dots\Gamma(\sum_{i=1}^N y_{n,i,d-d_{\boldsymbol{\gamma}}}+\alpha)}{\Gamma(y_{n,. ,.}+(d-d_{\boldsymbol{\gamma}})\alpha)}\\
 \nonumber&\times \prod_{c\in C} \frac{\Gamma(\beta_1+\beta_2)}{\Gamma(\beta_1)\Gamma(\beta_2)} \frac{\Gamma(\beta_1+\sum_{c_i=c} y_{n,i,.})\Gamma(\beta_2+\sum_{c_i=c} y_{e,i,.})}{\Gamma(\beta_1+\beta_2+\sum_{c_i=c} y_i)}\\
 \nonumber&\times \frac{\Gamma(d_{\boldsymbol{\gamma}}\alpha)}{[\Gamma(\alpha)]^{d_{\boldsymbol{\gamma}}}}\frac{\Gamma(\sum_{c_i=c} y_{e,i,1}+\alpha) \dots \Gamma(\sum_{c_i=c} y_{e,i,d_{\boldsymbol{\gamma}}}+\alpha)}{\Gamma(\sum_{c_i=c} y_{e,i,.}+d_{\boldsymbol{\gamma}}\alpha)}. 
 \end{aligned}$$

### Phylogenetic tree node selection in the mixture of Dirichlet tree multinomials model

We assume that nodes not contributing to clustering are “noise" nodes, whose allocation probabilities are the same across clusters.  Others are informative nodes, whose allocation probabilities are cluster-specific.  We denote the set of informative nodes $J^d$, and “noise" nodes $J^n$. If the tree given is the taxonomic tree, the method can show at which level the bacterial composition begins to differentiate between different clusters. 

Denote all the parameters associated with the observation $i$ as $B_i$, with the part belonging to the “noise" nodes as $\boldsymbol{b}_{j}$, $j \in J^n$, and the part shared only by the observations in the same cluster $\boldsymbol{b}_{j}(c_i)$, $j \in J^d$, $c_i \in \{1,2,\dots,C\}$. The likelihood of observation $X_i$ is: $$\begin{aligned}
 P(X_i|B,T,\boldsymbol{\gamma},c_i)=&\prod_j\frac{n_{j,.}(X_i)!}{n_{j,1}(X_i)!n_{j,2}(X_i)!\dots n_{j,|K_j|}(X_i)!}\\
 \nonumber&\prod_{j \in J^n} \prod_{k \in K_j} b_{jk}^{n_{jk}(X_i)}\prod_{j \in J^d} \prod_{k \in K_j} b_{jk} (c_i)^{ n_{jk}(X_i)}. 
 \end{aligned}$$ The likelihood for all the observations $\mathbf{X}$ is: $$\begin{aligned}
 P(\mathbf{X}|B,T,\boldsymbol{\gamma},C)=&\prod_i \left[\prod_j\frac{n_{j,.}(X_i)!}{n_{j,1}(X_i)!n_{j,2}(X_i)!\dots n_{j,|K_j|}(X_i)!}  \right]\\
 \nonumber& \prod_{j \in J^n} \prod_{k \in K_j} b_{jk}^{n_{jk}(\mathbf{X})}\prod_{c \in C} \prod_{j \in J^d}\prod_{k \in K_j} b_{jk} (c)^{ n_{jk}(\mathbf{X}^c)}.  
 \end{aligned}$$

Here $n_{j,k}(\mathbf{X^c})$ is the total number of genetic sequences descending from node $j$ to node $k$ for the observations in cluster $c$.

Again, for simplicity of computation, we use a strategy of integrating out parameters.  We consider the simplest conjugate prior $\mathbf{b}_{j}|\alpha \sim \mathrm{Dirichlet}(\alpha)$, and integrate out $B$. The marginal likelihood is $$\begin{aligned}
 P(\mathbf{X}|T,\boldsymbol{\gamma},\mathbf{c})=& \prod_i \left[\prod_j\frac{n_{j,.}(X_i)!}{n_{j,1}(X_i)!n_{j,2}(X_i)!\dots n_{j,|K_j|}(X_i)!}  \right]\\
 \nonumber&\prod_{j \in J^n} \frac{\Gamma(|K_j| \alpha)}{[\Gamma(\alpha)]^{|K_j|}}
 \frac{(\prod_{k \in K_j}\Gamma(\alpha+n_{jk}(\mathbf{X})))}{\Gamma(|K_j| \alpha+\sum_{k \in K_j} n_{jk}(\mathbf{X}))}\\
 \nonumber&\prod_{c \in C}\prod_{j \in J^d} \frac{\Gamma(|K_j| \alpha)}{[\Gamma(\alpha)]^{|K_j|}}
 \frac{(\prod_{k \in K_j}\Gamma(\alpha+n_{jk}(\mathbf{X}^c)))}{\Gamma(|K_j| \alpha+\sum_{k \in K_j} n_{jk}(\mathbf{X}^c))}. 
 \end{aligned}$$

Mixture of finite mixtures
--------------------------

The MFM model can achieve clustering without pre-defining the number of clusters.  One advantage of the MFM model over the Dirichlet process mixture model is that it can give a consistent estimate of the number of clusters.  The hierarchy of this model is: $$\begin{aligned}
  M&\sim p_m&\text{where }p_m\text{ is a p.m.f on }\{1,2,\dots\}\\
% &&\text{on }\{1,2,\dots\}\\
\nonumber (\pi_1,\dots,\pi_M)&\sim\mathrm{Dirichlet}_m(\eta,\dots,\eta)& \text{given } M=m\\
\nonumber c_1,\dots,c_N & \sim \boldsymbol{\pi}&\text{given }\boldsymbol{\pi}\\
\nonumber \theta_1,\dots,\theta_m &\sim G_0& \text{given }M=m\\
\nonumber \mathbf{X_i}&\sim F(\theta_{c_i})\text{ for } i=1,\dots,N &\text{given }\theta_{1:M},c_{1:N}.\end{aligned}$$

Here, $M$ is the underlying number of components from the population.  The vector $\boldsymbol\pi=(\pi_1,\dots,\pi_M)$ is the probability of a random sample belonging to a component.  $c_1,\dots,c_N$ are the component indicators for the samples.  $\theta_i$ is the set of corresponding parameters for observation $i$. $G_0$ is the base distribution.  $F$ is the mixing kernel introduced in the above subsections, i.e., Dirichlet multinomial or Dirichlet tree multinomial including feature selection. 

Similar to the Dirichlet process mixture model, the underlying discrete measure of the MFM model has a Pólya urn scheme representation, which enables sampling the parameters for each observation sequentially.  This close parallelism makes most sampling algorithms designed for the Dirichlet process mixture model directly applicable.  The parameter set for observation $i$, $\theta_i$, will either take the identical value of an existing parameter, or a newly generated value from the base distribution $G_0$ with the following probabilities: $$\theta_i|\theta_{-i} \sim \sum_{c \in C} (n_{c,-i}+\eta)\delta(\theta_c^*)+ \frac{V_N(|C|+1)}{V_N(|C|)}\eta G_0,$$ where $\theta_{-i}$ are the parameters for all the observations except for the $i$th observation; $\theta_c^*$s are the distinct values of $\theta_{-i}$, and $n_{c,-i}$s are the corresponding numbers of observations having the parameter $\theta_c^*$, except for the $i$th observation.  The function $V_N(R)$ is defined as:

$$V_N(R)=\sum_{m=R}^{\infty}\frac{\Gamma(m+1) \Gamma(\eta m)}{\Gamma(m-R+1) \Gamma(\eta m+N)}p_M(m).$$

This completes the specification of the mixture of finite Dirichlet multinomial mixtures (MFMDM) and mixture of finite Dirichlet tree multinomial mixtures (MFMDTM) models. 

Method implementation {#sec:mcmc}
=====================

Obtaining a sample from the posterior distribution of either model requires the use of Markov chain Monte Carlo (MCMC).  In each MCMC iteration, we first select OTUs or phylogenetic tree nodes by updating the latent indicator $\boldsymbol{\gamma}$ given the current cluster assignments, then fix $\boldsymbol{\gamma}$ and apply the split-and-merge algorithm to assign observations into clusters.  Here we provide a high-level description of the algorithm, with additional details provided in the supplemental material. 

Selection of OTUs or tree nodes
-------------------------------

The latent selection indicator $\boldsymbol{\gamma}$ is updated by repeating the following Metropolis step $t$ times (@KimTade06 suggested $t=20$).  A new candidate $\boldsymbol{\gamma}^{new}$ is generated by randomly choosing one of the two transition moves:

1.   **add/delete** by randomly picking one of the $d$ indices in $\boldsymbol{\gamma}^{old}$ and changing its value (from $0$ to $1$ or from $1$ to $0$);

2.   **swap** by randomly drawing a 0 and a 1 in $\boldsymbol{\gamma}^{old}$ and switching their values. 

The new candidate is accepted with probability $
\min\left\{1,\frac{f(\boldsymbol{\gamma}^{new}|\mathbf{X},\mathbf{c})}{f(\boldsymbol{\gamma}^{old}|\mathbf{X},\mathbf{c})}\right\}
$, where $\boldsymbol{c}$ is the cluster assignment vector.  As $f(\boldsymbol{\gamma}|\mathbf{X},\boldsymbol{c})\propto f(\mathbf{X}|\boldsymbol{\gamma},\boldsymbol{c}) \mathrm{P}(\boldsymbol{\gamma})$, the proposed acceptance probability can be calculated by $$f(\boldsymbol{\gamma}^{new}|\mathbf{X},\boldsymbol{c})/f(\boldsymbol{\gamma}^{old}|\mathbf{X},\boldsymbol{c}) = f(\mathbf{X}|\boldsymbol{\gamma}^{new},\boldsymbol{c}) \mathrm{P}(\boldsymbol{\gamma}^{new}) / f(\mathbf{X}|\boldsymbol{\gamma}^{old},\mathbf{c}) \mathrm{P}(\boldsymbol{\gamma}^{old}).$$

Cluster assignment
------------------

We update the latent sample allocation vector $\boldsymbol{c}$ using @JainNeal04’s split-and-merge algorithm by first selecting two distinct observations, $i$ and $l$ uniformly at random.  Let $\mathscr{C}$ denote the set of other observations that are in the same cluster with $i$ or $l$.

If $\mathscr{C}$ is empty, we use the simple random split-merge algorithm.  Otherwise, we use the restricted Gibbs sampling split-merge algorithm.  Both involve a Metropolis-Hasting sampling step, with acceptance probability $$a(\boldsymbol{c}^{merge},\mathbf{c})=\min\left\{1, \frac{q(\boldsymbol{c}|\boldsymbol{c}^{merge})\mathrm{P}(\boldsymbol{c}^{merge})L(\boldsymbol{c}^{merge}|\mathbf{X},\boldsymbol{\gamma})}{q(\boldsymbol{c}^{merge}|\boldsymbol{c})\boldsymbol{P}(\boldsymbol{c})L(\boldsymbol{c}|\mathbf{X},\boldsymbol{\gamma})}
\right\} \text{if $c_i \neq c_l$, and}$$

$$a(\mathbf{c}^{split}|\mathbf{c})=\min\left\{ 
1, \frac{q(\boldsymbol{c}|\boldsymbol{c}^{split})\mathrm{P}(\mathbf{c}^{split})L(\mathbf{c}^{split}|\mathbf{X},\boldsymbol{\gamma})}{q(\boldsymbol{c}^{split}|\boldsymbol{c})\mathrm{P}(\boldsymbol{c})L(\boldsymbol{c}|\mathbf{X},\boldsymbol{\gamma})}
\right\} \text{if $c_i = c_l$}.$$

For the simple random split-merge algorithm, $q(\mathbf{c}|\mathbf{c}^{merge})/q(\mathbf{c}^{merge}|\mathbf{c}) =1$, $  q(\mathbf{c}|\mathbf{c}^{split}) / q(\mathbf{c}^{split}|\mathbf{c}) = 1$. For the restricted Gibbs sampling, we first randomly create a launch state.  This launch state is modified by a series of “intermediate" restricted Gibbs sampling steps to achieve a reasonable split of the observations.  The last launch state is used for the calculation of the transient probabilities.  Details of the split-and-merge algorithm are described in the supplemental material. 

The prior ratio, $\mathrm{P}(\mathbf{c}^{merge})/\mathrm{P}(\mathbf{c})$ or $\mathrm{P}(\mathbf{c}^{split})/\mathrm{P}(\mathbf{c})$, relies on the partition distribution $\mathrm{P}(\mathbf{c})$. In a MFM model, the probability function of $\mathbf{c}$ is: $$\mathrm{P}(\mathbf{c})=V_N(|C|)\prod_{c \in C}\eta^{(n_c)}l$$ where the probabilities of spliting a cluster and combining two clusters are: $$\frac{\mathrm{P}(\mathbf{c}^{split})}{\mathrm{P}(\mathbf{c})}=\frac{V_N(|C|+1)}{V_N(|C|)}\frac{\Gamma(n_{c_1}+\eta)\Gamma(n_{c_2}+\eta)}{\Gamma(n_{c_1}+n_{c_2}+\eta)\Gamma(\eta)};$$ $$\frac{\mathrm{P}(\mathbf{c}^{merge})}{\mathrm{P}(\mathbf{c})}=\frac{V_N(|C|-1)}{V_N(|C|)}\frac{\Gamma(n_{c_1}+n_{c_2}+\eta)\Gamma(\eta)}{\Gamma(n_{c_1}+\eta)\Gamma(n_{c_2}+\eta)},$$ where $n_{c_1}$ and $n_{c_2}$ are the number of observations in the two clusters. 

Post-processing of MCMC samples
-------------------------------

The sampled values of the clustering indicators can only describe whether two observations belong to the same cluster, but are not comparable between iterations, as the same indicator value may represent different clusters due to the “label switching" issue.  In this paper, we adopt @FritIcks09’s method for summarizing posterior cluster labels from MCMC samples.  We denote the proposed clustering estimate as $c^*$, and estimate the probability that sample $i$ and $j$ belong to the same cluster from $M$ MCMC samples by $\zeta_{ij}=\frac{1}{M}\sum_{m=1}^M I(c_i^{(m)}=c_j^{(m)})$. The estimated cluster assignment can be obtained by maximizing the adjusted Rand index: $$AR(c^*,\zeta)=\frac{\sum_{i<j}\mathrm{I}_{\{c_i^*=c_j^*\}}\zeta_{ij}-\sum_{i<j}\mathrm{I}_{\{c_i^*=c_j^*\}}\sum_{i<j}\zeta_{ij}/C_2^N}
{\frac{1}{2}[\sum_{i<j}\mathrm{I}_{\{c_i^*=c_j^*\}}+\sum_{i<j}\zeta_{ij}]-\sum_{i<j}\mathrm{I}_{\{c_i^*=c_j^*\}}\sum_{i<j}\zeta_{ij}/C_2^N}.$$ This method can handle the label-switching issue, and can be simply implemented using the R package “mcclust" [@Frit12]. 

Simulation study and real data applications {#sec:application}
===========================================

In this section, we plan to first introduce @DeFiCava10’s dataset, which serves as the basis for our simulation study and a verification of our method, then apply the MFMDM method to @RiquZhan19’s clinical dataset. 

@DeFiCava10 performed 16S rRNA gene sequencing, which produced 406,431 reads split over 29 individuals.  The sequencing depth varies from 7,730 to 22,987 with a median of 13,523.  We applied the UNOISE2 function in usearch [@Edga16] and identified 2,803 OTUs in this dataset.  Figure S3 in supplemental material shows the PCoA plots of this dataset using different metrics: Euclidean distance, Bray Curtis distance, unweighted UniFrac distance, and weighted UniFrac distance.  In this case, unweighted Unifrac distance offers a better visual separation between the two geographic clusters. 

Simulation study
----------------

We compare our models to several machine-learning based methods, including PAM and hierarchical clustering (i.e., hcut) with complete linkage using the Euclidean, Bray Curtis, unweighted UniFrac, and weighted UniFrac metrics.  We generated data with structure similar to the De Filippo dataset, consisting of two groups (A and B) and 2803 OTUs.  Each group has 15 observations and each observation has 15000 sequences.  We simulated 5 scenarios, with decreasing levels of complexity, and generate the OTU counts for the $z$th scenario, $z=1,2,\dots,5$, in the following way. 

1.   Find two non-overlapping subsets of OTUs, $\Psi$ and $\Lambda$. In our simulation set-up, one subset $\Psi$ accounts for $13\%$ of the total genetic sequences and $356$ OTUs, while the other subset $\Lambda$ accounts for $15\%$ of the genetic sequences and $595$ OTUs. 

2.   Set the expected abundance of the two groups using the marginal distribution of the De Filippo dataset.  Move $z/5$ of the probabilities from the descending OTUs of $\Psi$ to $\Lambda$ for group A, and $z/5$ of the probabilities from the descending OTUs of $\Lambda$ to $\Psi$ for group B.

3.   Generate the number of genetic sequences from a Dirichlet multinomial with the sum of parameters to be $200$. (In the De Filippo dataset, if we assume the data are from two Dirichlet multinomial distributions for the two geographic groups, the sums of the parameters are estimated to be $191.81$ for the African group and $234.60$ for the Italian group.)  The phylogenetic tree information used in the mixture of Dirichlet tree distributions is the phylogenetic tree of the De Filippo dataset. 

4.   Repeat the above steps 200 times. 

The parameters of the Beta distribution in the mixture of Dirichlet multinomial distributions are set to be $\beta_1=\beta_2=1$, which corresponds to a uniform prior on the number of genetic sequences which are informative for clustering.  Similarly, the parameter of the Dirichlet distribution for both the MFMDM and MFMDTM are set to be $\alpha=1$, which is a uniform prior in the multinomial case.  The prior probabilities of OTUs being informative is $50\%$, i.e., $w=0.5$ for MFMDM model.  We give $M-1$ a $\mathrm{Poisson}(1)$ distribution, which expresses a preference for a small number of clusters. 

In both models, a large number of observed sequences inflates the factorial terms in the likelihood, which tends to support finer clusters.  To temper this effect and better achieve meaningful clustering, we take an approach similar to @GrieMcDa18 who “normalize" the data by first dividing the observed counts by a scaling parameter.  For the simulated data, we found $50$ is a reasonable scaling parameter.  In general, the maximum sequencing depth divided by $300$ is recommended. 

![ Comparison with Machine Learning Methods in Terms of Rand Indices, and the AUC of the High Abundance OTUs[]{data-label="randAUC"}](plots/capture){width="\linewidth"}

[0.32]{}

![ Comparison with Machine Learning Methods in Terms of Rand Indices, and the AUC of the High Abundance OTUs[]{data-label="randAUC"}](plots/othermethodsbar_10002){width="\textwidth"}

[0.32]{}

![ Comparison with Machine Learning Methods in Terms of Rand Indices, and the AUC of the High Abundance OTUs[]{data-label="randAUC"}](plots/othermethodsbar_10004){width="\textwidth"}

[0.32]{}

![ Comparison with Machine Learning Methods in Terms of Rand Indices, and the AUC of the High Abundance OTUs[]{data-label="randAUC"}](plots/AUC2_1000){width="\textwidth"}

For each dataset, we run 20,000 iterations, the first 10,000 of which are discarded as burn-in, and then apply a thinning of $10$ and keep 1,000 samples for inference.  For machine learning methods, the number of clusters is determined by the silhouette width, which is appropriate for non-Euclidean distances.  We measure clustering performance using the adjusted Rand index.  The expected value of the adjusted Rand index is 0 when clustering is done at random, while 1 reflects perfect recovery of the true underlying clusters in the data.  Figure \[randAUC\] (a) and (b) show the performance of our proposed methods, compared to some machine learning methods for scenarios 2 and 4, where the heights of the bars in the plot represent the median over 200 datasets and the black intervals represent the empirical estimate of the $95\%$ confidence interval.  Barplots for scenario 1, 3 and 5 can be found in supplementary Figure S4.  When the separation between two clusters is relatively small, the proposed MFMDTM method, which performs variable selection accounting for the phylogenetic tree, shows a big advantage over the machine learning based methods including PAM using UniFrac distances, which incorporate phylogenetic information (Figure \[randAUC\] (a)).  When the separation between the clusters becomes larger, the performance of MFMDM still shows significant improvement over that of competing methods that do not account for phlyogeny (Figure \[randAUC\] (b)).  In general, the methods which incorporate the tree information outperform those that do not, with MFMDTM achieving the highest adjusted Rand indices across all methods considered. 

An advantage of the proposed methods over the alternatives is that they enable the selection of informative features.  The inference about informative vs. noisy OTUs is based on the marginal posterior distribution of the latent indicator $\boldsymbol{\gamma}$, which is estimated from the selection frequencies in the MCMC output [@KimTade06].  It is worth mentioning that $\Psi$ and $\Lambda$ contain OTUs with low abundance, whose effects are negligible compared with the simulation noise.  To set a meaningful goal for selection, we consider the $37$ high abundance OTUs which differ across groups, from among the $197$ high abundance OTUs in the dataset, as the true discriminatory features, where “high abundance" is defined as marginal abundance greater than $0.001$. As shown in Figure \[randAUC\] (c), the AUCs of the ROC curves suggest that MFMDM’s selections are successful even when the Rand index is low.  The decrease in the AUC for Scenario 5 is due to the fact that the larger separation enables detection of informative OTUs that are below $0.001$, which reduces the specificity.  These results show that the proposed method is able to accurately recover the informative features. 

Application to De Filippo data as a validation
----------------------------------------------

We now summarize the application of the MFMDM and MFMDTM methods to the @DeFiCava10 dataset.  All parameter choices and MCMC options for the Bayesian methods are as in the simulation studies provided above.  Again, for the machine-learning methods, the number of clusters is determined by the silhouette width.  Assuming the geographic location as the ground truth for cluster assignment, the resulting Rand indices, which summarize clustering success, demonstrate that our proposed parametric models outperform commonly used machine-learning methods, such as partitioning around medoids (PAM) or hierarchical clustering with complete linkage (hcut) using Euclidean, Bray Curtis and UniFrac distances (Figure \[realdata\] (a)).  The clusters given by our methods are finer than the geographical groupings, but the observations from different countries have a very low probability of being assigned in the same cluster (Figure \[realdata\] (b)-(c)). 

[0.32]{}

![ MFMDTM Heatmap](plots/randrealdata){width="\textwidth"}

[0.32]{}

![ MFMDTM Heatmap](plots/MFMDM_1000){width="\textwidth"}

[0.32]{}

![ MFMDTM Heatmap](plots/MFMDTM_1000){width="\textwidth"}

![ MFMDTM Heatmap](plots/keys){width="\textwidth"}

The tree plot in Figure \[realdata\] (d) shows the node “Bacteria" is selected with a high probability, which indicates a difference between the African samples and Italian samples at the phylum level.  This confirms @DeFiCava10’s finding that African children’s microbiota showed a significant enrichment in Bacteroidetes and Actinobacteria and depletion in Firmicutes.  The original paper also found Italian children’s microbiota lacked genus Prevotella and Xylanibacter, whose ancestor node Bacteroidales is selected in our model.  Also, they identified that family Enterobacteriaceae were significantly underrepresented in African children vs. Italian children, and our model selects the corresponding phylum Proteobacteria and class Gammaproteobacteria as informative nodes for clustering.  On a computer with an Intel Core i5-6500 3.20GHz processor and 16GB memory, it takes Rcpp 25 minutes for 1000 MCMC iterations with the MFMDM model, and MATLAB R2016b 26 minutes for 1000 MCMC iterations with the MFMDTM model.  We also implemented the Dirichlet process mixture of Dirichlet multinomials (DPDM) and Dirichlet tree multinomials (DPDTM).  The results from the simulation and real data application are similar to the results from the MFM models (Supplemental Figure S1 and S2). 

[0.7]{}

![ Performance Measurements Based on Proposed Methods for De Filippo Dataset\
[]{data-label="realdata"}](plots/selectedTree_1000){width="\textwidth"}

Application to Riquelme data
----------------------------

Pancreatic ductal adenocarcinoma (PDAC) is one of the most deadly cancer types.  Most patients are diagnosed with an advanced stage, and the 5 year overall survival probability is below $10\%$ [@SiegMill19].  Though the prognosis is dismal, a minority of patients can survival more than 5 years after surgery.  The goal of this case study is to identify microbiome features that characterize a specially identified set of long-term survivors in contrast to short term survivors. 

In the data set provided by @RiquZhan19, we identified 2,410 OTUs corresponding to 1,095 taxonomic units above the species level.  We applied the MFMDM and MFMDTM methods to this data set, using the same parameter settings as in the simulation and previous case study.  For the 43 patients treated at MD Anderson (MDA) Cancer Center included in the study, we have the exact survival times for 42 of them, while for patients treated at John Hopkins, only a binary indicator of whether a patient survived more than 5 years was provided.  For patients with exact survival times, we plot a heatmap showing the posterior probability of any two observations belonging to the same cluster, where samples are sorted by the exact survival times (Figure \[RiquHeatmap\]).  The heatmap shows that long-time survivors’ microbiota are more homogeneous than that of short-time survivors.  This finding suggests that it may be fruitful for researchers to investigate the specific bacteria present in long-term survivors, as reflecting a distinctive protective microbiome state.  A similar conclusion can be drawn from the heatmap of all the samples using both the MFMDM and MFMDTM methods (Supplemental Figure S5).  Though some pre-clinical models [@PushHund18; @AykuPush19] have suggested that certain microbial species are positively associated with tumor progression, our finding is aligned with @RiquZhan19’s conclusion that a protective microbiome induces anti-tumor immunity in long-time survivors and those protective species are the key for future interventions. 

![ Heatmap Showing the Posterior Probability of Being Assigned to the Same Cluster for Patients with Exact Survival Time[]{data-label="RiquHeatmap"}](plots/realdataMDAOnly "fig:"){width="\textwidth"} ![ Heatmap Showing the Posterior Probability of Being Assigned to the Same Cluster for Patients with Exact Survival Time[]{data-label="RiquHeatmap"}](plots/keys "fig:"){width="\textwidth"}

Our model also identifies 40 phylogenetic tree nodes that drive the clustering of the samples.  Figure \[featureSelection\] shows the selected tree nodes with red color, the majority of which were also identified in @RiquZhan19.  The original paper showed the predominance of Clostridia in STS patients and Alphaproteobacteria in LTS patients at the class level, while our method shows that the two corresponding phyla, Firmicutes and Proteobacteria, are differential across clusters with posterior probability greater than $95\%$. @RiquZhan19 identified the species Bacillus clausii as predictive of survivorship, and our method selects the genus it belongs to as a relevant feature.  Some of the taxa identified by MFMDTM and @RiquZhan19 were discovered by other previous research on PDAC.  For example, [@FarrZhan12] identified that the abundance of the genus Corynebacterium is lower in PDAC patients than healthy individuals, while our method specifically points out that the sequence allocation probability at the family Corynebacteriaceae is differential across clusters. [ @GellBarz17] discovered that Proteobacteria producing cytidine deaminase are most associated with pancreatic cancer and our model identified several features under the phylum Proteobacteria, including the family Porphyromonadaceae and Enterobacteriaceae.  In addition, MFMDTM identifies features that have not been thoroughly discussed in pancreatic cancer literature before, such as the order Rhizobiales, which was found in higher abundance among patients with Helicobacter pylori-negative intestinal metaplasia than those with Helicobacter pylori-negative chronic superficial gastritis or cancer [@ParkLee19]. 

Compared with the LEfSe method used in the original paper, which tends to select nested features, our method can identify the exact taxonomic level at which clusters are different.  For completeness, we plot the heatmap of the conditional probability of allocating a sequence to child nodes for selected nodes (Supplement Figure S6) and the PCoA plot of the samples colored by the cluster assignment derived from the posterior of MFMDTM (Supplement Figure S7). 

Discussion {#sec:discussion}
==========

We have proposed two novel approaches for clustering of microbiome samples.  Unlike existing approaches, our methods perform variable selection, which enables biological understanding of features that differentiate clusters present in the data, without the need for pre-specification of the number of clusters.  The simulation results and application to real data demonstrate that our methods can outperform commonly used unsupervised clustering algorithms in terms of Rand indices, suggesting that the sparse models are not only more interpretable, but also more robust to noise.  Our application to gut microbiome profiling of pancreatic cancer patients enhances the originally published analysis of this dataset: while [@RiquZhan19] applied LefSe to identify microbiome features assuming known group membership, through our approach we show that the long-term survivors comprise a more natural cluster, while the short-term survivors are more heterogeneous, and identify additional features that are differential across the inferred clusters.  These findings could guide the development of future microbiome interventions to improve cancer outcomes, which are an active and exciting area of current medical research [@Ret2019; @Mcq2019]. 

![ MFMDM Results on @RiquZhan19 Data](plots/selectedTree2){width="85.00000%"}

Software {#software .unnumbered}
========

The Dirichlet multinomial mixture model code is included in the R package BayesianMicrobiome, while the Dirichlet tree multinomial mixture is implemented using Matlab.  Both are available at <https://github.com/YushuShi/BayesianMicrobiome.git>. 

Acknowledgements {#acknowledgements .unnumbered}
================

KAD is partially supported by NIH/NCI CCSG grant P30CA016672, MD Anderson Moon Shot Programs, Prostate Cancer SPORE grant P50CA140388, CCTS grant 5UL1TR000371, and CPRIT RP160693.  CBP is partially supported by NIH/NCI CCSG grant P30CA016672 and MD Anderson Moon Shot Programs.  RRJ is partially supported by NIH grant R01 HL124112 and CPRIT grant RR160089. 

[1]{}

Aykut, B., Pushalkar, S., Chen, R., et al.(2019).  The fungal mycobiome promotes pancreatic oncogenesis via activation of [M]{}[B]{}[L]{}. ,  574(7777):264–267. 

Bray, J. R. and Curtis, J. T. (1957).  An ordination of the upland forest communities of southern [W]{}isconsin. ,  27(4):325–349. 

Chen, J. and Li, H. (2013).  Variable selection for sparse [D]{}irichlet-multinomial regression with an application to microbiome data analysis. ,  7(1). 

De Filippo, C., Cavalieri, D., Di Paola, M.,et al. ( 2010).  Impact of diet in shaping gut microbiota revealed by a comparative study in children from [E]{}urope and rural [A]{}frica. ,  107(33):14691–14696. 

Dennis III, S. Y. (1991).  On the hyper-[D]{}irichlet type 1 and hyper-[L]{}iouville distributions. ,  20(12):4069–4081. 

D’Argenio, V. (2018).  Human microbiome acquisition and bioinformatic challenges in metagenomic studies. ,  19(2):383. 

Edgar, R. C. (2016). :  improved error-correction for [I]{}llumina 16[S]{} and its amplicon sequencing. . 

Farrell, J. J., Zhang, L., Zhou, H.,et al. ( 2012).  Variations of oral microbiota are associated with pancreatic diseases including pancreatic cancer. ,  61(4):582–588. 

Fritsch, A. (2012). .  R package version 1.0. 

Fritsch, A. and Ickstadt, K. (2009).  Improved criteria for clustering based on the posterior similarity matrix. ,  4(2):367–391. 

Geller, L. T., Barzily-Rokni, M., Danino, T., et al. ( 2017).  Potential role of intratumor bacteria in mediating tumor resistance to the chemotherapeutic drug gemcitabine. ,  357(6356):1156–1160. 

George, E. I. and McCulloch, R. E. (1993).  Variable selection via gibbs sampling. ,  88(423):881–889. 

Gilbert, J. A., Blaser, M. J., Caporaso, J. G., et al. ( 2018).  Current understanding of the human microbiome. ,  24(4):392–400. 

Goodrich, J. K., Di Rienzi, S. C., Poole, A. C., et al. ( 2014).  Conducting a microbiome study. ,  158(2):250–262. 

Gopalakrishnan, V., Spencer, C., Nezi, L.,et al. ( 2018).  Gut microbiome modulates response to anti–pd-1 immunotherapy in melanoma patients. ,  359(6371):97–103. 

Grier, A., McDavid, A., Wang, B., et al. ( 2018).  Neonatal gut and respiratory microbiota: coordinated development through time and space. ,  6(1):193. 

Haffari, G. and Teh, Y. W. (2009).  Hierarchical [D]{}irichlet trees for information retrieval.  In [*Human Language Technologies: The 2009 Annual Conference of the North American Chapter of the ACL, pages Boulder, Colorado, June 2009*]{} 173–181. 

Holmes, I., Harris, K., and Quince, C. (2012).  Dirichlet multinomial mixtures: Generative models for microbial metagenomics. ,  7(2):1–15. 

Huang, R., Yu, G., Wang, Z.,et al. ( 2013).  Dirichlet process mixture model for document clustering with feature partition. ,  25(8):1748–1759. 

Jain, S. and Neal, R. M. (2004).  A split-merge [Markov]{} chain [Monte Carlo]{} procedure for the [Dirichlet]{} process mixture model. ,  13(1):158–182. 

Kaufman, L. and Rousseeuw, P. J. (2008). ,  chapter 2, pages 68–125.  Wiley-Blackwell. 

Kim, S., Tadesse, M. G., and Vannucci, M. (2006).  Variable selection in clustering via [Dirichlet]{} process mixture models. ,  93(4):877–893. 

Knight, R., Callewaert, C., Marotz, C., et al. ( 2017).  The microbiome and human biology. ,  18:65–86. 

Li, Q., Guindani, M., Reich, B. J., et al. ( 2017).  A [B]{}ayesian mixture model for clustering and selection of feature occurrence rates under mean constraints. ,  10(6):393–409. 

Liu, Y., Chen, J., Su, Z., et al. ( 2016).  Robust head pose estimation using [D]{}irichlet-tree distribution enhanced random forests. ,  173:42–53. 

Love, M. I., Huber, W., and Anders, S. (2014).  Moderated estimation of fold change and dispersion for [RNA]{}-seq data with [DESeq]{}2. ,  15(12):550. 

Lozupone, C. and Knight, R. (2005).  Unifrac: a new phylogenetic method for comparing microbial communities. ,  71(12):8228–8235. 

Lozupone, C. A., Hamady, M., Kelley, S. T., et al. ( 2007).  Quantitative and qualitative $\beta$ diversity measures lead to different insights into factors that structure microbial communities. ,  73(5):1576–1585. 

MacQueen, J. (1967).  Some methods for classification and analysis of multivariate observations.  In [*Proceedings of the fifth Berkeley symposium on mathematical statistics and probability*]{}, volume 1, pages 281–297.  Oakland, CA, USA. 

Madigan, D., York, J., and Allard, D. (1995).  Bayesian graphical models for discrete data. ,  63(2):215–232. 

Mart[í]{}nez, I., Stegen, J. C., Maldonado-G[ó]{}mez, et al. ( 2015).  The gut microbiota of rural [P]{}apua [N]{}ew [G]{}uineans: composition, diversity patterns, and ecological processes. ,  11(4):527–538. 

McQuade, J. L., Daniel, C. R., Helmink, B. A., et al. ( 2019).  Modulating the microbiome to improve therapeutic response in cancer. ,  20(2):e77–e91. 

Miller, J. W. and Harrison, M. T. (2014).  Inconsistency of [P]{}itman-[Y]{}or process mixtures for the number of components. ,  15:3333–3370. 

Miller, J. W. and Harrison, M. T. (2018).  Mixture models with a prior on the number of components. ,  113(521):340–356. 

Park, C. H., Lee, A., Lee, Y., et al. ( 2019).  Evaluation of gastric microbiome and metagenomic function in patients with intestinal metaplasia using 16[S]{} r[RNA]{} gene sequencing. ,  24(1):e12547. 

Price, M. N., Dehal, P. S., and Arkin, A. P. (2010).  Fast[T]{}ree 2 [-]{} approximately maximum-likelihood trees for large alignments. ,  5(3):e9490. 

Pushalkar, S., Hundeyin, M., Daley, D.,et al. ( 2018).  The pancreatic cancer microbiome promotes oncogenesis by induction of innate and adaptive immune suppression. ,  8(4):403–416. 

Reticker-Flynn, N. E. and Engleman, E. G. (2019).  A gut punch fights cancer and infection. ,  565:573–574. 

Richardson, S. and Green, P. J. (1997).  On [B]{}ayesian analysis of mixtures with an unknown number of components (with discussion). ,  59(4):731–792. 

Riquelme, E., Zhang, Y., Zhang, L., et al. ( 2019).  Tumor microbiome diversity and composition influence pancreatic cancer outcomes. ,  178(4):795 – 806.e12. 

Rousseeuw, P. J. (1987).  Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. ,  20:53–65. 

Siegel, R. L., Miller, K. D., and Jemal, A. (2019).  Cancer statistics, 2019. ,  69(1):7–34. 

Tadesse, M. G., Sha, N., and Vannucci, M. (2005).  Bayesian variable selection in clustering high-dimensional data. ,  100(470):602–617. 

Tang, Y., Ma, L., and Nicolae, D. L. (2018).  A phylogenetic scan test on a [D]{}irichlet-tree multinomial model for microbiome data. ,  12(1):1–26. 

Tang, Y. and Nicolae, D. L. (2017).  Mixed effect [D]{}irichlet-tree multinomial for longitudinal microbiome data and weight prediction. 

Wang, T. and Zhao, H. (2017).  A [D]{}irichlet-tree multinomial regression model for associating dietary nutrients with gut microorganisms. ,  73(3):792–801. 

[^1]: Department of Statistics, University of Missouri, Columbia, Columbia, MO , 65201, USA. 

[^2]: Department of Biostatistics, University of Texas MD Anderson Cancer Center, Houston, TX , 77030, USA. 

[^3]: Department of Genomic Medicine, University of Texas MD Anderson Cancer Center, Houston, TX , 77030, USA. 

[^4]: To whom correspondence should be addressed.