A negative binomial latent factor model for paired microbiome sequencing data

A negative binomial latent factor model for paired microbiome sequencing data | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results A negative binomial latent factor model for paired microbiome sequencing data Hyotae Kim , Nazema Siddiqui , Lisa Karstens , Li Ma doi: https://doi.org/10.1101/2024.12.01.626246 Hyotae Kim 1 Department of Biostatistics & Bioinformatics, Duke University Roles: Postdoctoral Associate Find this author on Google Scholar Find this author on PubMed Search for this author on this site Nazema Siddiqui 2 Department of Obstetrics and Gynecology, Duke University Roles: Associate Professor Find this author on Google Scholar Find this author on PubMed Search for this author on this site Lisa Karstens 3 Department of Medical Informatics and Clinical Epidemiology and Department of Obstetrics and Gynecology, Oregon Health & Science University Roles: Associate Professor Find this author on Google Scholar Find this author on PubMed Search for this author on this site Li Ma 4 Department of Statistical Science and Department of Biostatistics and Bioinformatics, Duke University Roles: Professor Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: li.ma{at}duke.edu Abstract Full Text Info/History Metrics Preview PDF Abstract Motivation Microbiome compositional data are often collected from several body sites and exhibit dependency among them. Analyzing microbial compositions from different sites jointly allows for effective borrowing of information by exploiting the underlying cross-site correlation, which can lead to more effective statistical analysis, especially when the sample size at one or both sites is limited. To this end, we introduce a joint model for microbiome compositions at two (or more) sites within the same subjects. Our model incorporates (i) latent factors shared across two body sites to explain the common subject effects and to serve as the source of correlation between the two sites; and (ii) mixtures of latent factors to allow heterogeneity among the samples in their level of cross-site association. The model is illustrated with synthetic data and we apply it in a case study involving samples of the urinary and vaginal microbiome collected from women. Results Simulation studies show how common subject effects influence regression analysis results; a stronger association between two sites in the data causes a greater degree of bias in the analysis. The model with latent factors mitigates the bias present in the model without latent factors, whereas the two models perform comparably for the data set without paired associations. In a case study involving samples collected from a study on the female urogenital microbiome with aging (e.g., the UMICRO study), our model leads to the detection of covariate associations of the vaginal and urinary microbiome composition that are otherwise not statistically significant under a similar regression model applied to the two sites separately. Our model also enables prediction of the microbial abundance at one site based on observations from another site. We also consider a model extension that allows the clustering of subjects (samples) and cluster-specific levels of paired association. Under the extended modeling framework, the clusters can be classified according to their association strengths. 1 Introduction The advent of next-generation sequencing (NGS) technology enables the identification of a wide range of microbes from environmental samples without the need for cultivation, thus facilitating the exploration of microbial communities. The two most widely used methods for sequencing microbial communities are universal marker gene amplicon sequencing, such as the 16S rRNA gene, and whole metagenome shotgun sequencing (WMS) of all microbial genomes. These methods produce sequencing reads, which are subsequently mapped to taxa at various taxonomic levels using bioinformatic preprocessing pipelines, such as DADA2 [ Callahan et al., 2016 ]and MetaPhlAn [ Blanco-Míguez et al., 2023 ]. As a result of this taxonomic profiling, a large, highly sparse count table of taxa per sample is produced, with typically a finer taxonomic resolution in WMS than in 16S amplicon sequencing. To gain a comprehensive understanding of the human microbiome, researchers often collect data from multiple body sites for each individual and compare the composition and functions of the microbial communities in different parts of the body (e.g., [ HMPC, 2012a , HMPC, 2012b ]). The UMICRO data set, which we will use as a case study, includes vaginal-urine paired samples obtained by vaginal swabbing of the distal vagina and urine collection by transurethral catheterization, with the aim of identifying the microbial composition of the two communities and analyzing their variations with respect to the menopausal status of the participants. As the sample pairs were collected from two different body sites of the same subject, common effects associated with the subject are likely to be present in both vaginal and urine samples. This data set motivates the development of a model to capture the potential associations between the vaginal and urine niches. For the analysis of microbial compositional data in the form of a count table with rows for samples and columns for taxa, modeling methods originally developed for RNA sequencing (RNA-seq) or single-cell RNA sequencing (scRNA-seq) data can be adopted. Similarly to microbial compositional data, RNA sequencing data provide sparse count tables that report the number of sequence fragments assigned to each gene per sample or cell, for which the following models have been devised: negative binomial models [ Robinson et al., 2010 , Love et al., 2014 ],zero-inflated negative binomial models [ Risso et al., 2018 ], Poisson zero-inflated log-normal models [ Wang et al., 2018 ],and truncated Gaussian hurdle model [ Finak et al., 2015 ]. In addition, [ Martin et al., 2020 , Morton et al., 2019 , Paulson et al., 20 Sohn et al., 2015 ]introduced beta-binomial, multinomial, and zero-inflated Gaussian models, which were specifically developed for microbiome data. Although zero-inflated models were created to accommodate the pronounced sparsity of microbiome data, their validity in count-based models for sequencing reads remains controversial. As discussed in [ Silverman et al., 2020 , Sarkar and Stephens, 2021 ], in many common settings involving sparse sequencing count data, the abundance of zeros can often be adequately accommodated by simply incorporating overdispersion into count-based sampling models, such as in negative binomial models, without needing an additional zero-inflation component. We share this viewpoint, and because we are considering overdispersed count-based models in this paper, we do not by default incorporate an additional zero-inflation component, but in cases where such a component is indeed justified, incorporating it into the model is straightforward. We propose a model-based approach to jointly account for microbiome compositions at two (or more) sites, adopting a negative binomial regression model for each site while incorporating a shared latent factor to parsimoniously capture potential correlation in paired samples from different body sites in some taxa. Specifically, two negative binomial distributions, one for each site, share a set of latent factors, which are interpreted as unobserved common effects that contribute to their correlation. Under suitable choices of priors on the model parameters, posterior samples can be drawn efficiently using a Gibbs sampler that employs a data-augmentation technique called Pólya-Gamma augmentation. We extend our model on the latent factors to accommodate the common assumption that the samples are not homogeneous but rather form subgroups or clusters, each with their own cross-site correlation patterns. We present versions of the model that work both when the subgrouping is observed and when it is not. In the latter case the subgrouping is inferred from the data. We carry out a case study on the UMICRO data in which the study subjects come from three subgroups by study design. The rest of the paper is organized as follows. Section 2 introduces our negative binomial latent variable model with Section 2.1 for the base modeling framework and Section 2.2 for some model extensions. In Section 2.3 , we present a hierarchical representation of our model with prior specifications for the model parameters, followed by a posterior prediction approach. The proposed model is illustrated through synthetic data in Section 3.1 and the UMICRO study in Section 3.2. Finally, Section 4 concludes. 2 Methods 2.1 Joint Negative Binomial Model For a given taxon (e.g., genus), we use y si to denote the observed count in sample i at body site s , with i = 1, ⋯, n and s = 1, 2. Note that the data are actually also indexed by taxa, but for simplicity, throughout our description of the model and the computational recipes, we suppress the index for taxa, as our model is taxon-specific and is applied to each of the taxa separately. Let NB( μ, α ) be a negative binomial distribution with mean μ and variance (1 + μα ) μ . We consider the following joint negative binomial model (JNBM) that relates the counts from two sampling sites through a taxon-specific latent factor, γ i , in the form of a multiplicative factor on the mean count for the taxon, where X i denotes the vector of covariates for sample i , and β s is the regression coefficient vector for body site s. α s is the overdispersion parameter of the negative binomial distribution. N si indicates the total number of read counts, that is, the sum of y si over all taxa of interest. The latent factor γ i (or its exponential exp( γ i ))is an unobserved variable that represents a site-invariant sample-specific effect. The multiplicative effect of exp( γ i ) on the mean, μ si , has E(exp( γ i )) = 1 and Var(exp( γ i )) = exp( ϕ 2 ) − 1. This distributional assumption with the mean of 1 enables the random effects of exp( γ i ) to have no preference for a positive impact ( > 1) or negative impact ( < 1). The parameter ϕ 2 measures the strength of the association between the two body sites; when ϕ 2 becomes zero, the model is reduced to two independent negative binomial distributions for each body site, that is, NB for s = 1, 2. We will refer to it as the separate negative binomial model (SNBM). Section 3.2 will compare our joint two-site model with SNBM in a case study on the UMICRO data. To perform JNBM-based Bayesian inference, we need priors for ( α s , β s , ϕ 2 ); a description of the complete Bayesian hierarchical model can be found in Section 2.3.1. A fully conjugate sampling recipe for posterior inference based on the Pólya-Gamma data augmentation technique is detailed in Supplementary Material B. For brevity in the model description, we assume balanced paired data, i.e., both sites have the same number of samples. However, the model can also be used for unbalanced data where some samples are available only for one site. 2.2 Model Extension Although JNBM assumes a constant level of cross-site association throughout all paired samples, in real-world examples, this assumption is often unrealistic. However, allowing each sample to have its own level of association may lead to an overly flexible model. We thus compromise and assume that there are subgroups (either observed or unobserved) among the samples for which the extent of cross-site association is comparable. Consequently, we extend the base joint model by allowing different sample groups to exhibit different levels of paired association, which can easily be achieved by replacing the normal distribution for the latent factors with a mixture distribution. The following is a twocomponent normal mixture-based model, which randomly assigns paired samples to two groups: It should be noted that we enforce an inequality constraint on the mixing parameter, and , to ensure model identifiability. The general form of the mixture prior is given by with and for L > 1. In the paper, we will call this joint negative binomial model with the two-component mixture distribution JNBM Mix. Supplementary Material A.4 presents JNBM Mix’s superior predictive performance over JNBM, along with random clustering results for Lactobacillus samples as an example. The clustering of paired samples can also be accomplished by using a given variable g ( i ) = { 1, ⋯, L} rather than random assignment with {ν l } , where differences in association levels across L groups may serve as an indicator of the proximity among them. The mixture distribution for γ i is adjusted as follows: where δ l ( x ) is a delta function with δ l ( x ) = 1 for x = l and 0 otherwise. For the UMICRO case study, we use Study Group variable for clustering, such that g ( i ) = { Postmenopausal with no estrogen, Postmenopausal on estrogen, Premenopausal } , which model is named JNBM SG. Figure S11 shows the posterior distributions of for each study group under the model, revealing obvious similarities between the two postmenopausal groups and differences in location or/and dispersion between postmenopausal and premenopausal for Lactobacillus, Aerococcus, Prevotella and Corynebacterium . JNBM SG also outperforms SNBM and JNBM in prediction, which will be discussed in Section 3.2 and Supplementary Material A.2. 2.3 Posterior inference 2.3.1 Hierarchical model representation Below is a full description of the joint negative binomial model (JNBM) presented in (1) of Section 2.1. Let y = {y si : s = 1, 2 and i = 1, ⋯, n} , γ = {γ i : i = 1, ⋯, n} , and β s = {β sp : p = 1, ⋯, P} . The model is defined as where . For simplicity, we place exponential priors on their rate parameters were chosen empirically, which are conservative choices and lead to substantial prior-to-posterior learning. For example, in Section 3.1, we set for a simulation study with true regression coefficients of ( β 11 , β 12 , β 13 ) = (0.05, 0.001, 0.03) and ( β 21 , β 22 , β 23 ) = (0.03, 0.002, 0.01). With the prior means of , the 95% prior uncertainty bands for β s 1 , β s 2 , and β s 3 are given by ( − 2.5, 1.5), ( − 0.7, 0.6), and ( − 2.5, 1.5), respectively. The conservative choices result in much tighter 95% posterior uncertainty bands with substantial shifts of the posterior means (from the prior means) toward the true values. Likewise, we chose conservative hyperparameters of . For the case study in Section 3.2, we set the hyperparameters to , as well as for coefficients of the in tercept, two continuous covariates, and five binary covariates. We observed that this prior specification gives rise to such prior-to-posterior learning, although the true underlying parameter values are unknown in the case study. As in γ i , the regression coefficients β sp are assigned normal distribution priors with mean and variance , allowing the multiplicative effects of exp( β sp ) on μ si to have their means equal to 1 with their variances of exp . The following section discusses the augmented likelihood under the negative binomial modeling framework, which enables the regression coefficients and the latent factors to have full conditionals in closed form. Similarly, we can represent the extended models – JNBM Mix and JNBM SG – with variations in distributional assumptions on γ i as follows, [ JNBM Mix ] : with auxiliary variables {ξ i } for a hierarchical representation of the mixture distribution on γ i , where Dir( p 1 , ⋯, p L ) denotes a Dirichlet distribution with p l = 1 /L . [ JNBM SG ] : where g ( i ) is an observed variable, e.g., the study group clinical variable of the UMICRO data in Section 3.2. The hyperparameter of the exponential prior for is chosen in the same manner as a ϕ 2 of ϕ 2 in JNBM. 2.3.2 Cross-site prediction In making predictive inferences about unknown observables y ∗ , we consider the posterior predictive distribution below, from which we draw predictive samples and compute their average. The posterior predictive distribution is where y denotes the observed data and θ is a vector of model parameters: ( α, β ) for SNBM and ( α, β, γ ) for JNBM. Posterior predictive samples for y ∗ are drawn from the negative binomial sampling distribution p ( y ∗ | θ ) with parameters that are substituted with the posterior samples of θ taken from p ( θ | y ). Then, the mean of the posterior predictive samples becomes the predictive value for y ∗ . One useful feature of our joint model is the ability to predict counts at one site using observations from the other site. Suppose, for example, that y 1 i ′ for i ′ ∈ I ⊂ A ≡ { 1, ⋯, n} are unobserved and prediction targets. In SNBM, posterior sampling of model parameters ( α 1 , β 1 ) is based only on observations {y 1 i : i ∈ A − I} , while, in JNBM, not only {y 1 i : i ∈ A − I} but also y 2 i ′, i ′ ∈ I , are used for ( α 1 , β 1 , γ i ′ ) posterior sampling, specifically for latent factors γ i ′ . By plugging the posterior samples of the model parameters into the sampling distribution NB , with γ i ′ = 0 for SNBM, we can draw posterior predictive samples for y 1 i ′ . Let be a posterior predictive sample from NB with acquired at b -th MCMC iteration for b = 1, ⋯, B . Then, the predictive value for y 1 i ′ is defined as the posterior mean, that is, . Section 3.2 discusses the predictive performance of the models, which is based on the leave-one-out cross-validation (LOOCV) approach; we select one observation from the 2 n observations across n pairs of samples as a test set (i.e., an unknown observable for prediction), fit the models to the remaining data, and predict the test set value. We repeated the LOOCV 2 n times to obtain the predictive values for all sample pairs of a given taxon. 3 Results 3.1 Simulation We simulate paired samples with varying levels of association, comparing separate and joint negative binomial models to empirically grasp the role of latent factors in the joint model. Synthetic data were generated by drawing n = 300 pairs of samples from negative binomial distributions with latent factors defined as, The dispersion parameters and regression coefficients are arbitrarily chosen to be α 1 = 2, α 2 = 1, β 1 = ( β 11 , β 12 , β 13 ) = (0.05, 0.001, 0.03), and β 2 = ( β 21 , β 22 , β 23 ) = (0.03, 0.002, 0.01),= where X i = ( x i 1 , x i 2 , x i 3 ) ′ is a vector consisting of an intercept x i 1 , a continuous covariate x i 2 ∼ Unif(20, 80) and a binary covariate x i 3 ∼ Bern(0.3). Unif( a, b ) and Bern( p ) are (continuous) uniform and Bernoulli distributions with means ( a + b ) / 2 and p , respectively. Again, the dispersion parameter ϕ 2 of the latent factors in the underlying distribution determines how strong the association between pairs of samples is, and it is set to 0 for independent pairs and to 2 or 10 for dependent pairs, with 10 indicating a stronger association. N 1 i and N 2 i are offsets, which in real-world applications are used to normalize for sequencing depth (or the sum of abundances across selected taxa of interest). We set sampling of N si from the negative binomial distribution results in count data with a right-skewness as with the sequencing depths in the UMICRO data set (but on a smaller scale). By repeatedly sampling with the parameters, covariates, and offsets defined above, we produced 300 sets of 300 sample pairs for the following analysis. We evaluate the separate and joint models with regard to the accuracy of their estimated regression coefficients. Table 1 shows the mean square error (MSE) between the estimated and true coefficients for the two covariates under each model, demonstrating that JNBM is superior to SNBM for the data with non-zero associations ( ϕ 2 ≠ 0). This indicates that ignoring the latent sample effect ( γ i ) as in SNBM can severely reduce the efficiency in estimating the underlying relationships between the response variables and predictors. As expected, the difference in estimation accuracy between the two models becomes greater as the paired association in the underlying distribution becomes stronger, and therefore the common effects become larger. View this table: View inline View popup Download powerpoint Table 1. MSE, denoted as with site s = 1, 2 and covariate p = 2, 3 for regression coefficients {β sp } , with its standard error in parentheses. 3.2 Case Study: the UMICRO Data We apply JNBM to the UMICRO data set and compare the resulting inferences to those from SNBM in the estimation of regression coefficients and in the prediction of taxa abundance at one of the sites. Again, data were collected from two different body sites by catheterization for urine samples and vaginal swabbing for vaginal samples. Full-length contigs received from Loop Genomics (Element Biosciences, San Diego, CA), which we refer to as LoopSeq, were processed with DADA2 (v 1.24.0) to generate amplicon sequence variants (ASV), using parameters as recommended for synthetic full-length 16S LoopSeq data in [ Callahan et al., 2021 ]. Taxonomic classifications were assigned using BLCA (v 2.2) and the 16S NCBI database (downloaded on 11/16/2021). The processed data contain 68 sample pairs with 151 common taxa (genera) found in both vaginal and urinary samples, while the following analysis focuses on nine taxa of interest: Gardnerella, Streptococcus, Lactobacillus, Aerococcus, Anaerococcus, Bifidobacterium, Corynebacterium, Fannyhessea, Prevotella . We incorporate several clinical and demographic metadata for the subjects as covariates in the models: age, body mass index (BMI), diabetes (yes/no), daily yogurt or probiotic consumption (yes/no), race ethnicity (0: White, 1: Black, 2: Others), the presence or absence of overactive bladder (OAB) (yes/no). In addition, the study group variable on menopausal status (Postmenopausal with no estrogen; Postmenopausal on estrogen; Premenopausal) is used to cluster the samples for JNBM SG, introduced in Section 2.2. Figure 1 shows examples of differences between SNBM and JNBM when it comes to the estimation of regression coefficients. The results indicate that JNBM has relatively pronounced positive effects of BMI and age on Gardnerella and Streptococcus counts in both vaginal and urinary samples, while these effects are negligible under SNBM. In the previous section, we saw that the joint model provides more accurate estimation results when paired associations are present. Here, the posterior means (and standard deviations) of ϕ 2 for Gardnerella and Streptococcus are given by 17 (4) and 14 (4), respectively. It favors the results of JNBM for the regression coefficients. Furthermore, the positive effects of BMI and age on the taxa under JNBM are consistent with the biological findings reported in [ Brookheart et al., 2019 , Xu et al., 2020 ]. More results on regression coefficients can be found in Supplementary Material A.1. Download figure Open in new tab Figure 1. Posterior distributions of regression coefficients for BMI in the vaginal and urine data sets for Gardnerella (first two panels); and regression coefficients for Age in the data sets for Streptococcus (last two panels). The dashed lines indicate the posterior means. JNBM SG clusters the samples according to the study group variable, allowing each group to have a different level of paired association. Figure 2 illustrates that, for the three genera, the two postmenopausal groups have similar association strengths but differ from the premenopausal group, which has a weaker association. Figure S11 displays the posterior distributions of the parameters for the nine taxa; interestingly, Gardnerella and Fannyhessea exhibit remarkable differences in strength among the three groups, with postmenopausal with no estrogen showing the strongest association, followed by postmenopausal with estrogen. Download figure Open in new tab Figure 2. Posterior distributions of , the dispersion hyperparameters of latent factors, from JNBM SG for the three different genera. To assess the predictive performance of the models for the microbial composition, we calculated the residuals between the observed and predicted relative abundances of each microbial genus, denoted as , where y sik is an observed count for the k -th taxon in the i -th sample of the s -th body site and the corresponding predicted count under model M . Figure 3 presents the difference in the absolute residuals between SNBM and JNBM SG, that is, . The greater the difference, the better the predictive performance of JNBM SG. The boxplots represent the distributions of DARs per taxon for each body site, with the positive medians (in most taxa) indicating that JNBM SG enhances the predictive efficiency of SNBM. The improvement in prediction under the joint model is attributed to the borrowing of information between sites. As such, the base joint model, JNBM, also outperforms SNBM, which has been further improved by accommodating different levels of association for clusters in JNBM SG ( Figure S10 ) and in JNBM Mix ( Figure S12 ). Download figure Open in new tab Figure 3. Boxplots (from the 1st to 3rd quartiles) of , the differences in the absolute residuals between SNBM and JNBM SG, for each taxon. The genera on the y-axis are sorted by their median DARs in each data set. 4 Conclusions We have proposed a joint negative binomial model (JNBM), a latent variable model based on negative binomial distributions, for paired microbiome data. The UMICRO data, gathered from two different body sites for each of the subjects, is a motivating real-world example, and the proposed model incorporates a set of latent factors to capture associations between the two body sites. JNBM has been extended by modifying the distribution for the latent factors, which permits varying levels of paired association across groups of samples. Our joint negative binomial model provides more accurate regression estimates than the separate negative binomial model without latent factors when there exists a non-zero paired association in the data. In the UMICRO case study, the joint model reveals conspicuous positive effects of BMI and age on Gardnerella and Streptococcus abundances, respectively, which are in agreement with some previous research findings ([ Brookheart et al., 2019 , Xu et al.,2020 ]). Moreover, the joint model outperforms the separate model in prediction, which can be attributed to the use of observations from the opposite body site for prediction via latent factors. With extended models, prediction performance can be further enhanced due to their ability to assign different association strengths to each group. In addition, estimates of the association strengths can also be used to characterize sample groups. For example, we found that the two postmenopausal cohorts (with and without estrogen) had similar associations but were different from the premenopausal cohort, which had a relatively weaker association, for Lactobacillus, Aerococcus , and Prevotella in the UMICRO case study. Although the case study focuses on nine taxa of interest, we have applied these models to other genera. Unsurprisingly, the models struggled to predict taxa that are prevalent in one site but rare in the other. For example, Escherichia, Klebsiella , and Pseudomonas are pathogens that are common in the urine but scarce in the vagina, showing a high sparsity (high proportion of samples with zero count) in the vaginal data set. It is particularly challenging to predict non-zero counts (urine) using zero counts (vaginal) compared to the reverse scenario. The shared latent factor assumption under our model also becomes questionable when the taxon is rarely present in one of the sites but common in the other. Supplementary Material for “Negative binomial latent variable model for paired data in microbiome analysis” A. Additional figures for real data analysis A.1 Posterior distributions of β s Download figure Open in new tab Figure S1. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Gardnerella . Download figure Open in new tab Figure S2. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Streptococcus . Download figure Open in new tab Figure S3. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Lactobacillus . Download figure Open in new tab Figure S4. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Aerococcus . Download figure Open in new tab Figure S5. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Anaerococcus . Download figure Open in new tab Figure S6. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Bifidobacterium . Download figure Open in new tab Figure S7. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Corynebacterium . Download figure Open in new tab Figure S8. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Fannyhessea . Download figure Open in new tab Figure S9. Posterior distributions of regression coefficients in the vaginal (first row) and urine (second row) data sets for Prevotella . A.2 Boxplots of DARs Download figure Open in new tab Figure S10. Boxplots of (top row) and (bottom row). A.3 Posterior distributions of under JNBM SG Download figure Open in new tab Figure S11. Posterior distributions of colored by Study Group g ( i ). A.4 Prediction results of JNBM Mix Download figure Open in new tab Figure S12. Boxplots of . Download figure Open in new tab Figure S13. Lactobacillus: On the left is a scatter plot of observed relative abundances colored by posterior estimates of the auxiliary variable ξ i for the mixture prior; red for and green for , where is the posterior median estimate. The next two panels show marginal density estimates of the relative abundances, with observed values at the bottom. The last panel presents the posterior distributions of two dispersion hyperparameters and . B. Computational details on posterior inference B.1 Pólya-Gamma augmentation for negative binomial models This section describes the Pólya Gamma augmentation scheme for the proposed negative binomial model. For sample i at body site s , the negative binomial probability mass function is given by where under JNBM. Following [ Polson et al., 2013 ], the last term in the square brackets can be expressed as, where Z si ≡ log( μ si ) + log( α s ). PG( ω | a, b ) denotes the probability density function of the Pólya-Gamma distribution PG( a, b ). Using the (Pólya-Gamma) auxiliary variable ω si , the equation (6) can be represented hierarchically without the integration. Then, the joint probability function for y si and ω si is derived as, with for i = 1, ⋯, n . The augmented likelihood with the Pólya-Gamma auxiliary variables ensures (normal) prior conjugacy for ( β s , γ i ), facilitating posterior inference, which will be discussed in the following section. B.2 Posterior simulation The Gibbs sampler is used to draw posterior samples of model parameters for inference. The Pólya-Gamma augmentation enables us to obtain the full conditionals in closed form for most of JNBM parameters except α = {α 1 , α 2 }, ϕ 2 , and , for which we employ the Metropolis–Hastings algorithm. The following are the full conditionals for the parameters of our joint models. Let y s· = {y si : i = 1, ⋯, n} , γ = {γ i : i = 1, ⋯, n} , and ω s· = {ω si : i = 1, ⋯, n} for s = 1, 2. With the normal prior , the full conditional for regression coefficients, β s = ( β s 1 , β s 2 , ⋯, β sP ) ′ , is derived as, where and U s = ( X ′ Ω s X + T − 1 ) − 1 . X is a design matrix consisting of the intercept and the covariates for all samples, that is, X = ( X 1 , ⋯, X n ) ′ with X i = (1, x i 2 , ⋯, x ip ) ′ . Ω s indicates a diagonal matrix of ω si , such that Ω s = diag( ω s 1 , ⋯, ω sn ). T is also a diagonal matrix of diag . Similarly, the full conditional for the latent factors γ , with the normal distribution assumption of with , is given by where V = (Ω 1 + Ω 2 + Φ −1 ) − 1 and . Φ is a diagonal matrix of diag . The full conditional for the auxiliary variables of ω si , with the Polya-Gamma distribution PG prior, takes the form of According to Theorem 1 of [ Polson et al., 2013 ],this is proportional to a density function of a Polya-Gamma distribution PG ; using the prior conjugacy, the posterior samples of ω si can be taken from the updated Pólya-Gamma distribution. Other parameters, ( α , ϕ 2 , τ 2 ), have no closed-form full conditionals with the joint full conditional density Hence, the parameters are updated with Metropolis-Hastings steps in MCMC, using log-normal proposal distributions. For the extended models, the dispersion parameters , of latent factors are updated using the Metropolis-Hastings algorithm, too. Unlike JNBM SG, JNBM Mix has additional parameters {ν l } and {ξ i } , given by (4). As Dir is a conjugate prior for ( ν 1 , ⋯, ν L ), the mixture weight parameters are updated using the Dirichlet distribution with parameters , where | A | indicates the cardinality of set A . Finally, posterior samples of auxiliary variables {ξ i } can be drawn from the updated discrete probability function: , n and l = 1, ⋯, L . As our modeling strategy is taxon-by-taxon, posterior estimates of model parameters for the microbial community (of multiple taxa) can be obtained by repeating the above Gibbs samplers multiple times, but can be done in parallel. Acknowledgement The research was supported in part by NIGMS grant R01-GM135440, NSF grant EEC-2133504, and NIA grant R03 AG060082. Footnotes * Hyotae Kim ( hyotae.kim{at}duke.edu ), Nazema Siddiqui ( nazema.siddiqui{at}duke.edu ), Lisa Karstens( karstens{at}ohsu.edu ) References ↵ [ Blanco-Míguez et al., 2023 ] Blanco-Míguez , A. , Beghini , F. , Cumbo , F. , McIver , L. J. , Thompson , K. N. , Zolfo , M. , Manghi , P. , Dubois , L. , Huang , K. D. , Thomas , A. M. , et al. ( 2023 ). Extending and improving metagenomic taxonomic profiling with uncharacterized species using metaphlan 4 . Nature Biotechnology , pages 1 – 12 . ↵ [ Brookheart et al., 2019 ] Brookheart , R. T. , Lewis , W. G. , Peipert , J. F. , Lewis , A. L. , and Allsworth , J. E. ( 2019 ). Association between obesity and bacterial vaginosis as assessed by nugent score . American journal of obstetrics and gynecology , 220 ( 5 ): 476 – e1 . OpenUrl CrossRef ↵ [ Callahan et al., 2021 ] Callahan , B. J. , Grinevich , D. , Thakur , S. , Balamotis , M. A. , and Yehezkel , T. B. ( 2021 ). Ultra-accurate microbial amplicon sequencing with synthetic long reads . Microbiome , 9 ( 1 ): 130 . OpenUrl CrossRef PubMed ↵ [ Callahan et al., 2016 ] Callahan , B. J. , McMurdie , P. J. , Rosen , M. J. , Han , A. W. , Johnson , A. J. A. , and Holmes , S. P. ( 2016 ). Dada2: High-resolution sample inference from illumina amplicon data . Nature methods , 13 ( 7 ): 581 – 583 . OpenUrl CrossRef PubMed ↵ [ Finak et al., 2015 ] Finak , G. , McDavid , A. , Yajima , M. , Deng , J. , Gersuk , V. , Shalek , A. K. , Slichter , C. K. , Miller , H. W. , McElrath , M. J. , Prlic , M. , et al. ( 2015 ). Mast: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell rna sequencing data . Genome biology , 16 ( 1 ): 1 – 13 . OpenUrl CrossRef PubMed ↵ [ HMPC, 2012a ] HMPC ( 2012a ). A framework for human microbiome research . nature , 486 ( 7402 ): 215 – 221 . OpenUrl CrossRef PubMed Web of Science ↵ [ HMPC, 2012b ] HMPC ( 2012b ). Structure, function and diversity of the healthy human microbiome . nature , 486 ( 7402 ): 207 – 214 . OpenUrl CrossRef PubMed Web of Science ↵ [ Love et al., 2014 ] Love , M. I. , Huber , W. , and Anders , S. ( 2014 ). Moderated estimation of fold change and dispersion for rna-seq data with deseq2 . Genome biology , 15 ( 12 ): 1 – 21 . OpenUrl CrossRef PubMed ↵ [ Martin et al., 2020 ] Martin , B. D. , Witten , D. , and Willis , A. D. ( 2020 ). Modeling microbial abundances and dysbiosis with beta-binomial regression . The annals of applied statistics , 14 ( 1 ): 94 . OpenUrl PubMed ↵ [ Morton et al., 2019 ] Morton , J. T. , Marotz , C. , Washburne , A. , Silverman , J. , Zaramela , L. S. , Edlund , A. , Zengler , K. , and Knight , R. ( 2019 ). Establishing microbial composition measurement standards with reference frames . Nature communications , 10 ( 1 ): 2719 . OpenUrl CrossRef PubMed ↵ [ Paulson et al., 2013 ] Paulson , J. N. , Stine , O. C. , Bravo , H. C. , and Pop , M. ( 2013 ). Differential abundance analysis for microbial marker-gene surveys . Nature methods , 10 ( 12 ): 1200 – 1202 . OpenUrl CrossRef PubMed ↵ [ Polson et al., 2013 ] Polson , N. G. , Scott , J. G. , and Windle , J. ( 2013 ). Bayesian inference for logistic models using pólya–gamma latent variables . Journal of the American statistical Association , 108 ( 504 ): 1339 – 1349 . OpenUrl CrossRef ↵ [ Risso et al., 2018 ] Risso , D. , Perraudeau , F. , Gribkova , S. , Dudoit , S. , and Vert , J.-P. ( 2018 ). A general and flexible method for signal extraction from single-cell rna-seq data . Nature communications , 9 ( 1 ): 284 . OpenUrl CrossRef PubMed ↵ [ Robinson et al., 2010 ] Robinson , M. D. , McCarthy , D. J. , and Smyth , G. K. ( 2010 ). edger: a bioconductor package for differential expression analysis of digital gene expression data . Bioinformatics , 26 ( 1 ): 139 – 140 . OpenUrl CrossRef PubMed Web of Science ↵ [ Sarkar and Stephens, 2021 ] Sarkar , A. and Stephens , M. ( 2021 ). Separating measurement and expression models clarifies confusion in single-cell rna sequencing analysis . Nature genetics , 53 ( 6 ): 770 – 777 . OpenUrl CrossRef PubMed ↵ [ Silverman et al., 2020 ] Silverman , J. D. , Roche , K. , Mukherjee , S. , and David , L. A. ( 2020 ). Naught all zeros in sequence count data are the same . Computational and structural biotechnology journal , 18 : 2789 – 2798 . OpenUrl CrossRef ↵ [ Sohn et al., 2015 ] Sohn , M. B. , Du , R. , and An , L. ( 2015 ). A robust approach for identifying differentially abundant features in metagenomic samples . Bioinformatics , 31 ( 14 ): 2269 – 2275 . OpenUrl CrossRef PubMed ↵ [ Wang et al., 2018 ] Wang , J. , Huang , M. , Torre , E. , Dueck , H. , Shaffer , S. , Murray , J. , Raj , A. , Li , M. , and Zhang , N. R. ( 2018 ). Gene expression distribution deconvolution in single-cell rna sequencing . Proceedings of the National Academy of Sciences , 115 ( 28 ): E6437 – E6446 . OpenUrl Abstract / FREE Full Text ↵ [ Xu et al., 2020 ] Xu , J. , Bian , G. , Zheng , M. , Lu , G. , Chan , W.-Y. , Li , W. , Yang , K. , Chen , Z.-J. , and Du , Y. ( 2020 ). Fertility factors affect the vaginal microbiome in women of reproductive age . American Journal of Reproductive Immunology , 83 ( 4 ): e13220 . OpenUrl CrossRef Back to top Previous Next Posted December 05, 2024. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following A negative binomial latent factor model for paired microbiome sequencing data Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share A negative binomial latent factor model for paired microbiome sequencing data Hyotae Kim , Nazema Siddiqui , Lisa Karstens , Li Ma bioRxiv 2024.12.01.626246; doi: https://doi.org/10.1101/2024.12.01.626246 Share This Article: Copy Citation Tools A negative binomial latent factor model for paired microbiome sequencing data Hyotae Kim , Nazema Siddiqui , Lisa Karstens , Li Ma bioRxiv 2024.12.01.626246; doi: https://doi.org/10.1101/2024.12.01.626246 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7957) Biochemistry (18606) Bioengineering (14747) Bioinformatics (44059) Biophysics (22432) Cancer Biology (19561) Cell Biology (26713) Clinical Trials (138) Developmental Biology (13878) Ecology (20847) Epidemiology (2067) Evolutionary Biology (25270) Genetics (16086) Genomics (23369) Immunology (18563) Microbiology (42166) Molecular Biology (17927) Neuroscience (92710) Paleontology (693) Pathology (2964) Pharmacology and Toxicology (5054) Physiology (8043) Plant Biology (15886) Scientific Communication and Education (2090) Synthetic Biology (4532) Systems Biology (10169) Zoology (2371) window.__CF$cv$params={r:'a356b533fb396c9c',t:'MTc4ODQ1ODk5MA==',u:'01a06876901779638b7f50e02d69643d',ut:'lOZmSyf7PxftvtPU5cmye5HWbGZkgu_p3IUmu.RN9Ow-1788458995-1.2.1.1-bdVvwPrjQAwHBUWltmME7.hLwEA65ko2fG2hEUwMdyKIWIdPwI4ChA9K2gV657ghR4.yAa21NBtWgNH_mKHq90.NyEaYS9DDY1epBAchyyQ',i:60};(function(){if(!document.body)return;var s=document.createElement('script');s.src='/cdn-cgi/challenge-platform/scripts/precursor/main.js';document.head.appendChild(s);})();