4 Case studies
To better understand the potential impact of taxonomic bias on DA analysis in practice, we conducted several case studies spanning a range of biological scenarios and sequencing technologies.
4.1 Foliar fungi experiment
Measurement of control communities along with the primary experimental samples can enable researchers to directly measure and remove the effect of bias prior to or concurrent with downstream DA analysis (McLaren, Willis, and Callahan (2019)). Gnotobiotic community experiments are well suited to this form of calibration via community controls since, unlike most natural ecosystems, it is possible to assemble ‘mock communities’ containing all species in known relative abundances.
Leopold and Busby (2020) are the first to use mock-based calibration in a gnotobiotic microbiome experiment. In a study of host-commensal-pathogen interactions, Leopold and Busby (2020) inoculated plants with 8 commensal fungal species and subsequently exposed plants to a fungal pathogen. To investigate the joint effect of colonization order and host genotype, plants varied in source location and in which commensal colonized first. The authors used ITS amplicon sequencing to measure communities before and after pathogen infection, along with mock communities with quantified genome concentrations of DNA from the 9 species. These DNA mock communities allowed the authors to directly measure the bias due to the sequencing workflow following DNA extraction and correct for it when analyzing the primary experimental samples following the method of McLaren, Willis, and Callahan (2019). Although extraction is not accounted for, considerable bias might still be expected due to preferential PCR binding and amplification and the substantial ribosomal copy-number-variation (CNV) in fungi (Lofgren et al. (2019)). The authors found a 13X difference between the most and least efficiently measured commensal, while the pathogen was measured 40X more efficiently than the least efficiently measured commensal. We reanalyzed this experiment to understand the impact that bias might have had on the study’s analyses were it not accounted for.
Leopold and Busby (2020) used the ITS sequencing profiles of the pre-infection samples to understand the impact of host genotype and an experimental treatment—which species was allowed to colonize first—on the commensal community composition prior to infection. Two DA analyses were conducted to address these questions. First, the authors used linear regression of log species proportion (via negative-binomial regression to account for random read-sampling error) to ask whether host genotype, treatment, and their interaction significantly impacted overall community composition. Second, the authors quantified the strength of priority effects—the advantage of being allowed to colonize first—in different host genotypes by the LFC in the proportion of a species when allowed to colonize first versus later with the majority of commensals.
Since these analysis are both based on proportions, bias could in principle impact their results. We re-ran both analyses with and without bias correction, finding nearly identical results regardless of whether bias is accounted for6. To understand why, we compared the multiplicative variation taxa proportions with that of the mean efficiency. The proportions of each taxon varied substantially across samples, much more so than the mean efficiency, which had a geometric range of 1.62X and a geometric standard deviation of just 1.05X (SI Figure D.1). Why does the mean efficiency vary so little despite substantial multiplicative variation in the efficiency among taxa and in the taxa proportions across samples? This answer relates to the fact that the mean efficiency is a weighted arithmetic average of species proportions and so is insensitive to large multiplicative variation in low-proportion species. The pre-infection samples in the Leopold and Busby (2020) study are always dominated by the three species with the highest efficiencies, whose efficiencies vary only by 3X (and just 1.5X between the two most dominant). The lowest efficiency species, despite having large variation in their proportion on the log scale, are always rare (proportions \(\ll 1\), even after bias correction) and hence don’t significantly affect the mean efficiency. Because the log mean efficiency varies much less than the log proportions for each taxon, there is little room for it to distort the results of DA analyses based on variation in log proportions.
We additionally considered the potential for bias to impact analysis of how commensal taxa responded to infection. In particular, it is interesting to consider whether any commensals are able to increase in absolute concentration in response to infection. The experiment performed by Leopold and Busby (2020) does not allow absolute-abundance measurement; nevertheless, the results of the previous sections allow us to consider the impact that bias would have on an analysis of absolute genome concentration measured using the total-abundance normalization method (Equation (2.10)). It is first useful to consider how bias impacts the estimated change in log proportion of each commensal following infection. The pathogen is absent in pre-infection samples but tends to dominate the community post-infection, raising the mean efficiency of post-infection samples (SI Figure D.2). Across different host genotypes, the average increase in mean efficiency ranges from 2.5X to 5.2X. We used simple linear regression to measure the average change in log proportion of each commensal for each host genotype, with and without bias correction. Regardless of whether bias is accounted for, the log proportion of each commensal decreased, which is unsurprising given the pathogen’s growth and the sum-to-one constraint of species proportions. Without bias correction, however, the LFCs were lower by an amount corresponding to the inverse change in log mean efficiency. The magnitude errors created were substantial for a large fraction of commensal-host pairs; in several cases, the LFC estimates without bias correction were clearly negative, whereas the estimates with bias correction were close to 0.
Now, consider the impact this shift in mean efficiency might have on a regression analysis with log absolute genome concentration as the response, if measured using the total-abundance normalization method (Equation (2.10)). We consider two scenarios. Scenario 1: Total genome concentration in each sample is perfectly known and used to measure the concentrations of individual taxa via Equation (2.10). In this case, bias in the MGS measurements would create absolute errors in the LFCs for genome concentration identical to those for proportions. The scientific error, however, may be much worse. Since the sum-to-one constraint of proportions no longer applies, the LFCs will necessarily be larger (less negative), such that the absolute error caused by the change in mean efficiency could plausibly create substantial biological errors where some taxa that persist or increase in abundance instead appear to decrease. Scenario 2: Total ITS concentration is measured by qPCR with the same primers and PCR protocol used for ITS sequencing, and used without any correction to measure individual taxa concentration with the total-abundance normalization method. In this case, we should expect the increase in mean efficiency to have little to no impact on LFCs, since it would be offset by the increase in mean efficiency of the total-concentration measurement. Moreover, since the total ITS concentration and ITS sequencing are both performed on the extracted DNA, we can also expect the effect of taxonomic bias caused by DNA sequencing to be offset. Thus, the systematic shift in mean efficiency following infection may or may not significantly impact the LFC estimates depending on how the total-abundance measurement is made prior to normalization.
4.2 Vaginal microbiomes of pregnant women
The vaginal microbiome during pregnancy has been a source of intensive study due to its apparent connection with the health of both mother and her developing child. Many MGS studies have found associations of specific microbial species and community characteristics with rates of urinary tract and sexually-transmitted infections, bacterial vaginosis (BV), and preterm birth. Yet these associations vary across studies of different populations and using different MGS methods. DA analyses of the vaginal microbiome are commonly based on proportions, creating an opportunity for taxonomic bias to impact results. A number of studies have experimentally demonstrated substantial taxonomic bias among MGS protocols and individual steps (such as extraction and PCR amplification) in vaginal samples or in vitro samples of vaginally-associated species (Yuan et al. (2012),Brooks et al. (2015),Gill et al. (2016),Graspeuntner et al. (2018)). Bias has been proposed as a potential explanation for discrepancies across studies (Callahan et al. (2017), Others?), but there has so far been little quantitative analysis of this possibility. Here we use empirical bias measurements from control samples to investigate the role of bias in proportion-based DA analyses of vaginal microbiomes from a recent large-cohort study of pregnant women.
As part of the Multi-Omic Microbiome Study: Pregnancy Initiative (MOMS-PI) study, Fettweis et al. (2019) collected longitudinal samples from over 1500 pregnant women, including nearly 600 that were measured by amplicon sequencing of the 16S V1-V3 hypervariable region to yield species-level bacterial taxonomic profiles. Taxonomic bias of this MGS protocol was previously investigated by Brooks et al. (2015) and McLaren, Willis, and Callahan (2019), using measurements by Brooks et al. (2015) of cellular mock communities of seven common, clinically-relevant vaginal bacterial species. Of these, Lactobacillus iners had the highest efficiency, which was nearly 30X larger than that of the species with the lowest efficiency, Gardnerella vaginalis. A second Lactobacillus species, L. crispatus, had an efficiency that was approximately 2X less than L. iners and 15X greater than G. vaginalis. These species, along with the unculturable Lachnospiraceae BVAB1, are the most common top (most-abundant) species (SI analysis) and can reach high proportions in individual samples, indicating that shifts between them might drive large changes in the mean efficiency which might in turn impact DA results.
We sought to assess this possibility using a joint analysis of the control measurements from Brooks et al. (2015) and the microbiome measurements from the MOMS-PI study. To obtain the taxonomic bias among all species identified in the MOMS-PI measurements, we used the taxonomic relationships with the seven control species to impute the efficiencies for the remaining species. We used these imputed efficiencies to calibrate (correct the effect of bias in) the MOMS-PI measurements, examine variation in the mean efficiency varied across samples, and compare the results of a DA analysis with and without bias correction.
The mean efficiency varies substantially across vaginal samples (Figure 4.1. This variation appears to be primarily driven by variation in which species is most abundant, samples in which a Lactobacillus species is most abundant (after calibration) typically have an efficiency that is 3-20X greater than samples in which G. vaginalis is most abundant. Shifts between Lactobacillus-dominance and Gardnerella-dominance are common in between-women comparisons and occasionally occur between consecutive visits in individual women (SI Figure D.4). These shifts typically result in substantial fold changes in mean efficiency (SI Figure D.4) and can cause spurious fold changes in the trajectories of lower-abundance species (SI Figure D.5).
Figure 4.1: The mean efficiency in vaginal samples from the MOMS-PI study varies with the most abundant species.
We next considered whether the observed variation in mean efficiency could cause systematic error in a DA analysis. In particular, we hypothesized that DA analysis of species proportions versus a covariate that is associated with Lactobacillus and/or Gardnerella would be particularly prone to spurious results. Patient metadata was not available due to privacy restrictions; we therefore sought a clinically relevant covariate to use in a regression analysis that could be determined directly from the microbiome profiles and that we a priori expected to be associated with the proportions of Lactobacillus and Gardnerella. Alpha diversity metrics such as species richness, the Shannon index, and the (Inverse) Simpson index have been repeatedly found to be strongly positively associated with bacterial vaginosis (BV; Srinivasan et al. (2012), Cartwright et al. (2018)). Cartwright et al. (2018) found that an observed Simpson Diversity Index of 0.82 (corresponding to an order-2 effective number of species of 5.6) classified high diversity samples as BV positive with a sensitivity of 100% and specificity of 85.1%. In addition, it is commonly observed that samples from women with and without BV that are dominated by Lactobacillus spp. tend to have higher diversity, whereas samples dominated by Gardnerella tend to have lower diversity. We therefore chose to perform a DA analysis of species proportion versus alpha diversity, hypothesizing that Lactobacillus and Gardnerella would drive a negative association of mean efficiency with diversity and thereby distort DA estimates for all species. We split samples into low, medium, and high diversity groups based on Shannon diversity in observed (uncalibrated) microbiome profiles. We then estimated the LFC in proportion from low- to high-diversity samples for a 30 common species by simple linear regression, using calibrated (bias-corrected) and observed (uncorrected) microbiome profiles. As expected, the mean efficiency was higher in low-diversity samples due to a larger fraction of samples dominated by Lactobacillus and a lower fraction with samples dominated by Gardnerella. The decline in mean efficiency in the high diversity group caused a concomitant increase in the LFCs of all species when bias was not accounted for. This error due to bias resulted in serious magnitude errors in nearly all species and sign errors in over one third, with errors of both types seen for several clinically relevant species.
4.2.1 Notes
- It might be cleaner story-wise to redo the diversity partitioning using the threshold identified by Cartwright et al. (2018)
- For the DA analysis of proportions vs. diversity, I tried other types of DA analysis methods (Gamma-Poisson regression and regression on ranks, similar to Wilcoxon test) to understand how the impact of bias differs across methods. In each case bias has a big effect but the species whose estimates are most affected varies. I would need to revise and verify this analysis before including its results. One reason for doing so would be to show that our qualitative finding about how bias distorts results applies beyond the basic linear regression model.
- Decreases in mean efficiency during transitions from Lactobacillus to Gardnerella dominance can be expected to be even more extreme for commonly-used vaginal microbiome primers that fail to amplify Gardnerella.
4.3 Microbial growth in marine sediments
Our route for mean efficiency to become associated with the covariate is if the covariate quantifies a biological process that preferentially selects for a microbial trait that also tends to increase or decrease measurement efficiency. We illustrate this potential mechanism with a study of microbial growth in marine sediments.
The surface layers of marine sediments harbor diverse and abundant microbiomes. Total cell density and species richness decrease with depth as resources are consumed in older, deeper sediment layers; however, some taxa are able to grow and increase in density as they are slowly buried. Lloyd et al. (2020) performed a systematic assessment of growth rates of bacterial and archaeal taxa over a depth of 10 cm (corresponding to ~40 years of burial time) in sediment of the White Oak River estuary. To estimate growth rate, the authors first measured absolute cell density of microbial taxa using the total-abundance normalization method (Equation (2.10)), with taxa proportions measured with 16S amplicon sequencing and total community density measured by directly counting cells using epifluorescence microscopy (CARD-FISH). The authors refer to these absolute densities as FRAxC measurements, for ‘fraction of 16S reads times total cell counts’. The FRAxC measurements were used to infer growth rate from the slope of a simple linear regression of log cell density against burial time over the first 3 cm below the bioirrigation layer (corresponding to \(\sim 8\) years of burial). To validate the inference of positive growers, the authors compared the growth rates from FRAxC-based inference from two sediment cores, qPCR measurements in these cores for a few reference taxa, and FRAxC-based inference in two replicate laboratory incubation experiments.
Taxonomic bias could lead to systematic error in FRAxC-derived growth rates if sample mean efficiency tends to systematically vary with burial time (or equivalently, depth). One possibility is that microbes with tougher cell walls tend to persist longer (alive or dead) in the sediment, while at the same time being more difficult to extract DNA from than microbes with weaker cell walls. In this case, we would expect the relative abundance of tougher species to increase with depth and hence the mean extraction efficiency to decrease, which in turn would lead to inflated growth-rate estimates for all taxa; a possible end result would be that taxa that decay sufficiently slowly would be mistakenly inferred to have positive growth.
We can test the hypothesis that systematic variation in log mean efficiency with depth distorts inferred growth rates using the qPCR measurements of two clades also collected by Lloyd et al. (2020). In Section B, we argue that taxonomic bias in targeted qPCR measurements is expected to create a roughly constant fold error and yield fold changes in absolute density that are relatively unaffected by changes in mean efficiency. Thus comparing qPCR to FRAxC growth rates on the same reference taxa allows us to estimate the systematic error in FRAxC growth rates. Because the systematic error in regression slopes due to bias is the same for each taxon (Section 3), these comparisons allow us to draw conclusions about the accuracy of FRAxC growth rates for all taxa. The first soil core included qPCR measurements of a single archaeal clade, Bathyarchaeota, for which growth rates by qPCR and FRAxC were nearly identical (doubling rates of 0.099/yr by FRAxC and 0.097/yr by qPCR). The second soil core included qPCR measurements of Bathyarchaeota and a second clade, Thermoprofundales/MBG-D. In this core, FRAxC and qPCR growth rates differed more substantially, with growth rates from FRAxC being larger by 0.012/yr for Bathyarchaeota (0.112/yr by FRAxC and 0.1/yr by qPCR) and by 0.086/yr for Thermoprofundales/MBG-D (0.294/yr by FRAxC and 0.208/yr by qPCR). A low number of experimental samples and noise in both the FRAxC and qPCR measurements place significant uncertainty in these measurements; however, the fact that FRAxC-derived growth rates are larger than qPCR-derived rates in all three cases is consistent with our hypothesis that mean efficiency decreases with depth in a manner that systematically biases FRAxC-derived rates to higher values. The differences in growth rate are small in absolute terms; however, the maximum observed difference of 0.086/yr suggests an error large enough to impact results for some taxa classified as positive growers, whose FRAxC growth rates ranged between 0.04/yr and 0.5/yr. Overall, the comparison between FRAxC and qPCR measurements gives support to the study conclusions, but suggest that species at the lower end of this range of positive FRAxC-derived rates may in fact be merely persisting or even slowly declining in abundance.
4.4 Summary and conclusions
The impact of bias can depend on protocol, biological system, and type of DA analysis being done. Though these case studies span a highly limited range of possibilities, when combined with the theoretical results of section X help suggest some general conclusions about how and when bias will impact DA analyses based on fold changes in proportions.
First, the error caused by consistent taxonomic bias can mostly cancel in cross-sample comparisons and so not impact DA analyses of fold changes in species proportions. Our theoretical results from Section 3 indicate that when the mean efficiency is roughly constant across samples, the error in proportions cancels in fold change calculations and is absorbed by the intercept term in regression models. We observed this scenario in the analysis of pre-infection fungal microbiome samples of Leopold and Busby (2020) and when analyzing the trajectories of the vaginal microbiomes of women when the dominant species remained constant7, and a stable mean efficiency within marine soil cores is plausible and consistent with (though unable to be fully determined by) the different growth rate estimates in the Lloyd et al. (2020) experiment. Section 3 showed that absolute DA analysis using total-abundance normalization are susceptible to errors as with DA analysis of proportions; yet our analysis of two (hypothetical) approaches to analyzing absolute changes in response to infection in the fungal microbiome experiment indicates that the increasingly common approach of pairing marker-gene sequencing with qPCR of the total marker density can largely mitigate this effect, since here what matters is the variation in the ratio of mean efficiencies of the MGS measurement to the total-density measurement and this ratio may remain roughly constant if bias is largely shared by the two measurements.
Yet in other cases, the mean efficiency can vary substantially and create substantial error in fold changes between pairs of samples. We saw examples when comparing fungal microbiomes pre- and post-infection in the foliar-fungi experiment and comparing vaginal microbiomes with different dominant species in the MOMS-PI experiment. The impact of these errors on results of DA regression analysis across many samples depends on whether the mean efficiency varies systematically with the covariate of interest. In these two examples, systematic variation of the mean efficiency arose as high-efficiency species tended to dominate samples in one of the sample conditions (post-infection foliar samples or low-diversity vaginal samples). Thus one way for significant error in DA regression results to arise is when a small number of (or just one) species that have particularly high or low mean efficiencies, form a large proportion of the community in a substantial fraction of samples, and are associated with the covariate of interest.
For many experimental systems, however, there may not be any one species that frequently forms a large fraction of the community. In these systems, systematic variation of the mean efficiency can still arise through the collective change of many species that are associated with the covariate of interest. We described a hypothesized scenario in the marine-sediment case study in which increased burial time (the covariate) selects for lysis-resistant (and so lower efficiency) species. As another example, treatment with a certain antibiotic might selectively kill easier-to-lyze Gram negative species (or harder-to-lyse Gram positives) and thereby decreasing (or increasing) the mean efficiency in samples collected from a host post-treatment. Such situations can arise whenever there is a microbial trait that is associated with both measurement efficiency and the biological processes of interest. An example besides cell-wall structure is ribosomal copy number, which increases a species’ efficiency in ribosomal amplicon measurements and is positively associated with metabolic rate. A more generic mechanism by which mean-efficiency associations might arise is from the evolutionary relationships among species. Species with more recent common ancestry are expected to be more similar across a wide range of heritable traits, which include traits that affect measurement efficiency (such as cell wall structure, genome size, ribosomal copy number, PCR binding sequence) as well as traits that affect the biological processes under study. For example, differences in the phylum-level composition of samples from two conditions might be driven by many species that show phylum-level conservation in a relevant biological trait. If these species also show phylum-level conservation in (potentially different) traits that affect measurement efficiency, an association of mean efficiency with the condition can arise.
So long as the mean efficiency does not vary systematically with the regression covariates, the error in fold changes serve primarily to reduce the precision and power of regression analyses without systematically distorting estimates and so does not necessarily lead to invalid inferences. But the reduction in precision could be substantial, particularly given the small sample sizes common in many microbiome studies, and thereby substantially limit the ability to draw meaningful conclusions from a study8.
These observations provide reasons to both worry and hope. It seems likely that in many experiments the mean efficiency is consistent or is at least not associated with the covariate, so that DA inferences remain valid. Yet it is not obvious a priori in which studies this condition holds, and there are plausible mechanisms that can create problematic associations of the mean efficiency even in ecosystems with high species diversity. Thus while we should not discount the large set of existing DA results, we should seek ways to better assess the robustness of results from previous studies and to measure and correct the error caused by bias in future ones.
References
The two sets of analyses can be seen at the following links: negative-binomial regression analysis; priority-effects analysis↩︎
This point could use more support in the MOMSPI case study↩︎
I suspect this non-systematic variation in the mean efficiency is the more typical situation, but we currently lack an example of it in the case studies. We can find a semi-real example by taking real microbiome data from a case-control-type design and simulating a strong degree of bias. I did this with the HMP2 IBD study, but thinnk it might be better to find a different example where there are some clear DA results.↩︎