SNP-bacteria association screening

Our discovery dataset included variants from four genes (NOD2, CARD9, LRRK2, and ATG16L1, S1 Table) and bacterial levels from the gut obtained through 16S sequencing of 182 well-phenotyped IBD cases (Table 1). After stringent quality control (S1 Text, S1 Fig), 168 bacterial taxa remained for association testing. The relative proportion of these 168 taxa within each level is summarized in Fig 1A. A total of 6, 12, 17, 36, 63, and 34 taxa were included in each of the 6 hierarchical levels (phylum, class, order, family, genus and species, respectively). The bacterial microbiota was dominated by bacteria from the Firmicutes, Bacteroidetes and Proteobacteria phyla. When compared to a similar cohort of 38 healthy controls (who had complete microbiome data but no SNP data) using standard logistic regression (S2 Fig), we found that 73% of those bacteria displayed a negative association with IBD (Fig 1B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Overview of the bacterial taxa distribution in IBD cases and NOD2 association.

Panel a) shows the relative proportion of the 168 bacterial taxa analyzed across the six hierarchical levels. Only taxa representing more than 0.5% of the total bacterial loading are labelled. Grey areas correspond to unknown, unmeasured, or underrepresented taxa. Panel b) show the log of the odds ratio (OR) of the labelled bacteria on IBD case-control status and corresponding 95% confidence interval, as derived using standard logistic regression while adjusting for confounding factors.

https://doi.org/10.1371/journal.pgen.1008018.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Sample characteristics.

https://doi.org/10.1371/journal.pgen.1008018.t001

We performed systematic screening for associations between the four IBD-susceptibility genetic variants and bacterial quantification. Due to substantial correlation across hierarchical levels (S3 Fig), we conducted analyses within each level and applied Bonferroni correction to the p-value within each level to select the best candidate associations. An association test was performed using standard linear regression after adjusting for established confounding factors as well as covariates selected by the recently developed CMS approach [19] to increase the statistical power (S1 Text). We identified seven candidate associations (S2 Table), six of which involved NOD2 (with c_Bacteroidia, f_Bacteroidaceae, g_Bacteroides, g_Roseburia, R. faecis, and F. prausnitzii) and one that involved CARD9 (with p_Firmicutes). However, these associations corresponded to four independent signals (Table 1) with the same association being projected along branches of the phylogenic tree (see S4–S7 Figs for hierarchical plots of the associations). Indeed, g_Roseburia and R. faecis, which were found to be associated with NOD2, had a correlation of 0.84 in our dataset. Similarly, f_Bacteroidaceae and g_Bacteroides, which were also associated with NOD2, were almost identical and had a correlation of 0.93 with c_Bacteroidia. The two remaining associations were observed between NOD2 and F. prausnitzii and between CARD9 and p_Firmicutes. All four signals showed that the IBD risk alleles (of either NOD2 or CARD9) were negatively associated with the bacteria in question.

To validate these associations, we performed in silico replication using three independent datasets of adult human subjects with IBD including a total of 451 individuals [15, 21]. Genetic variant and bacteria data were only partly available in these datasets. CARD9 was available in only one of the three replication cohorts, and we used the genus g_Faecalibacterium as a proxy for F. prausnitzii in all three replication analyses. We preprocessed and analyzed the data similarly to that of the discovery sample and combined the two stages using a standard inverse-variance meta-analysis (Table 1). The replication analysis confirmed the associations between the NOD2 risk allele and decreased abundance of g_Roseburia and F. prausnitzii (Fig 2A); and both signals passed a stringent Bonferroni correction accounting for all tests performed at the meta-analysis stage (P < 7x10^-5). Firmicutes-CARD9 and Bacteroides-NOD2 replication showed a signal concordant with those observed during discovery (i.e., a negative association) but were not significant, indicating either a false positive in the discovery or a lack of statistical power in the replication. To explore further potential heterogeneity of the SNP effects, we conducted association analyses for the following disease subtypes: CD overall, Ileal CD, non-ileal CD, and ulcerative colitis (Table 1). The effect estimates for CD and ileal CD were consistent with the signal observed for IBD during both discovery and replication. Conversely, we observed heterogeneity of effects for both non-ileal CD and UC. Finally, we confirm further the NOD2 and F. prausnitzii association, which is a bacterium of particular interest for our group[22], using real-time quantitative PCR [23] in the discovery samples (Fig 2B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Summary of association signals between NOD2 variants and F. prausnitzii.

Panel a) presents the boxplots of rank-based inverse normalized F. prausnitzii from fecal microbiota in the discovery samples and the three replication datasets for non-carrier and carriers of NOD2 risk alleles. Panel b) presents F. prausnitzii in the healthy controls and cases from the discovery sample quantified using 16S qPCR and normalized to the global bacterial population (mean ± s.e.m.).

https://doi.org/10.1371/journal.pgen.1008018.g002

Model fit while accounting for selection bias

The associations between the IBD risk variants and two bacteria that are also associated with the IBD risk raise questions about potential reverse causation or conversely potential mediation of the genetic risk through bacteria. In an attempt to address the question of causality, we performed a series of analyses to assess the fitness of models matching existing knowledge, which can be summarized as follows. First, the effects of both NOD2 and CARD9 on IBD are well established [3]. Second, SNPs can only be explanatory variables, which implies a (direct or indirect) unidirectional relationship. Third, the four taxa from Table 1 (Roseburia genus, Bacteroidia class, F. prausnitzii, and Firmicutes phylum) have been consistently found to be negatively correlated with the IBD case status [8, 24]; however, the directionality of the effect remains a topic of investigation [25]. In regards of these three constraints, four (potentially overlapping) parsimonious underlying causal models can fit the data (Fig 3A–3D): (a) the effect of the risk allele is mediated by bacterial level; (b) the genetic variant influences both the disease and the bacteria, inducing a correlation between the two latter variables; (c) the risk allele and bacteria are independent risk factors of IBD; and (d) IBD mediates the effect of the risk allele on bacterial level (a so-called reverse causation). Note that these models remain relevant when adding intermediate factors. For example, NOD2 variants have been shown to correlate with lower concentrations of defensins [26], which in turn can be responsible for changes in the microbiota composition. Alternatively, IBD may be associated with increased epithelial oxygenation [27], which itself can be responsible for changes in the microbiota composition. Such scenarios will be in full agreement with model a) and model d), respectively, and therefore will not alter the conclusion of the analysis.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Hypothetical causal models.

Expected correlation between an IBD risk variant g and a bacteria B negatively correlated with IBD, across four potential underlying models. For illustrative purposes, g and B are assumed to be continuous and normally distributed and all effects are larger than those observed in real data. Top panels (a, b, c, and d) present the hypothetical causal diagrams and bottom panels (e, f, g, and h) present the corresponding scatterplots of B as a function of g in the population (dark grey points, and trend in black), and in cases only (light grey points and trend in red). In model a) the effect of g on IBD is mediated by B; in cases the effect of g on B is underestimated because of the oversampling of participants carrying risk alleles (e). In model b) the genetic variant influences both IBD and B, inducing a correlation between IBD and B which is observed in both the whole sample and cases only (f). In model c), g and B act independently on IBD and are therefore not associated in the population, however g and B are positively correlated in case-only samples because of biased selection (g). Finally, in model d) the effect of g on B is mediated by IBD; the indirect association between g and B observed in the general population is canceled when looking at cases only as the sample is stratified by the mediator (h).

https://doi.org/10.1371/journal.pgen.1008018.g003

Importantly, secondary phenotype analysis in highly ascertained samples (such as in the present study which includes IBD cases only) can introduce selection bias [28, 29], although previous work has shown that this bias will be minimal for rare diseases (i.e., those with a prevalence ≤ 1%, such as IBD [30]) and for cases with low to moderate associations between variables [29] (S3 Table). To assess the fitness of the data with the four causal models while accounting for selection bias, we derived an estimator of Δ, the bias induced by the IBD case-only sampling in the regression coefficient between the risk allele tested and the bacteria. As shown in S1 Text, Δ can be approximated in some specific situations (common diseases and modest effects of the predictors), but estimation in more general cases is extremely challenging. Nevertheless, this theoretical framework allowed us to determine the expected sign of the risk allele-bacteria association across the scenarios considered without relying on assumption on the data distributions. Following the counterfactual principle [20], we leveraged this property to infer which model had the highest fit for the data based on the signs of the observed associations. As illustrated in the toy example from Fig 3E–3H, our analytical derivation (see Materials and Methods) showed that the ascertainment could induce a positive bias of the SNP-bacteria association in model (c) and lead to the removal of the SNP-bacteria association in model (d). The results from this analysis are twofold. First, it shows strong evidence that among the four causal models considered, models (a) and (b), which included a true gene-bacteria effect, were the only models in which we expected a persistent negative SNP-bacteria effect, which was in agreement with the results from Table 1. Second, it rules out selection bias as a likely explanation for the observed association.

We then asked whether the proposed approximation of Δ could be also informative for the remaining gene-bacteria association results. Indeed, besides the main signals from Table 1, our data showed both enrichment for a negative bacterium-disease association (Fig 1B) and significant enrichment for negative effects of IBD risk alleles on all bacteria levels, with 65%, 71%, 70%, and 60% of the associations negative for ATG16L1, CARD9, LRRK2, and NOD2, respectively (Fig 4A–4D, p-values for enrichment equal 7.1x10^-5, 3.3x10^-8, 1.3x10^-8, and 5.3x10^-3, respectively). After aligning the bacterium-disease and bacterium-SNP associations, we observed a strong and significant (P = 5.0x10^-8) concordance of effects (i.e., most bacteria showed an association with SNPs and IBD with the same sign) (Fig 4E and 4F). However, our derivation and simulation study (S1 Text and S8–S12 Figs) demonstrated that, regardless of the error distribution of the bacterial level, the correlation between the bacterium-disease and the bacterium-SNP effect estimates in the case-only samples, is expected to be negative under model c), and null under model d). Therefore, the observed positive correlation (Fig 4F) provides further evidence against model c) (the risk allele and bacteria act independently on IBD) and d) (IBD mediates the effect of the risk allele on the bacteria) as generative models, and again indicates that models a) and b) are the best fits for our data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Correlation between SNP-taxa and IBD-taxa effect estimates.

The 168 bacterial taxa were tested for association with the variants from each of the four genes considered: (a) ATG16L1 (rs12994997), (b) CARD9 (rs10781499), (c) LRRK2 (rs11564258), and (d) NOD2 (rs2066844, rs2066845, and rs2066847). The histograms on the left panel show the distribution of IBD risk alleles-bacteria association (i.e. of , the regression coefficients) and the enrichment for negative effects (in blue, p-values equal 7.1x10^-5, 3.3x10^-8, 1.3x10^-8, and 5.3x10^-3, respectively). Panel (e) shows a similar histogram while merging the per-risk allele change in bacteria level of the four genes (i.e. summing for each bacteria the of the four genes). Panel (f) shows the distribution of bacteria-IBD association derived in a IBD cases-controls dataset () for each bin from panel (e). Together, panels (e) and (f) show the strong concordance of the gene-bacteria and bacteria-IBD effects, in agreement with a mediation effect of the risk allele on IBD through the microbiome. In particular, bacteria displaying lower level in carrier of IBD risk alleles are more likely to be negatively associated with the risk of IBD.

https://doi.org/10.1371/journal.pgen.1008018.g004

Fully deciphering models a) and b) would be extremely challenging because of the complex interrelationships between the variables in question and the untestable assumptions that would have to be assumed. However, we found some indications supporting model a) against model b) in our data. In particular, as illustrated in Fig 3, the association effect estimates between a genetic variant and bacteria are expected to be slightly deflated in the cases alone compared to that of the estimates derived in a population including controls in model a) but should remain unchanged in model b). To assess the fitness of the two models, we use the special case of NOD2, which is associated almost exclusively with CD [3]; thus, the NOD2-bacteria associations are expected to be unchanged between all IBD samples (including UC cases) and CD alone under model b) but should be different in model a). As shown in S13 Fig for all bacteria, we found that the effect estimates derived for the CD cases alone were significantly (P = 5x10^-3) smaller than those derived for the complete dataset, which matched the expected deflation observed in the presence of mediation and therefore provided an argument in favor of model a).

Sensitivity analyses

As most causal inference approaches, the proposed counterfactual strategy we used in this study has some limitations. In order to assess the validity of our simulation framework, we performed a series of sensitivity analyses that closely mimicked current knowledge on IBD associations. First, empirical bacterial levels follow a non-negative, overdispersed distribution and often harbor a large number of zero values; therefore, use of the Poisson distribution, negative binomial distributions, or hurdle models has been suggested for analyses of microbiome data [31]. Moreover, we applied a rank-based inverse-normal transformation in our real data analysis, although sparsity remained (after QC, the bacteria had a minimum of 20% nonzero values). To assess the impacts of these parameters, we performed a simulation study in which the bacterial levels were drawn from a negative binomial distribution with varying percentages of zero-values [0%; 95%]. For each replicate, the SNP-bacteria effect was estimated similar to the procedure used in the real data analysis. As shown in S14 Fig, no qualitative difference was found between the untransformed and transformed analyses over the four scenarios considered (i.e., Fig 3A–3D). Overall, both the true effect (models (a) and (b)) and bias (model (c)) decreased toward the null with the increased data sparsity.

Another potential limitation of our inference was that we assumed that all effects were homogeneous, and Crohn’s disease and ulcerative colitis cases were treated as a single outcome variable (IBD). However, as aforementioned, NOD2 variants and F. prausnitzii are mostly associated with CD, and ileal CD in particular, but show little evidence of an association with UC [3]. Such effect heterogeneity can be treated as a misclassification problem (i.e., UC cases are misclassified as CD cases). To assess the impact of genetic heterogeneity, we performed a second simulation in which the SNP-disease and bacteria-disease associations were present only for a subset of the cases (thus mimicking UC-CD heterogeneity). In this simulation, all parameters were set to ensure that the magnitudes of the known effects were similar to those reported in the literature for the NOD2-F. prausnitzii example and in our analysis. As described above, we considered the four causal models from Fig 3 and estimated for each replicate and each disease subtype strata the association between the genotype and the bacteria. As shown in S15 Fig, in this scenario, some of the bias from models (c) and (d) illustrated in Fig 3G and 3H decreased toward the null. Moreover, our simulations showed that under models (c) and (d), we expected no association signal from either the disease subtype strata or the whole IBD-cases sample. Again, only models a) and b) fit the observed results from both Table 1 and Fig 4.

Finally, the validity of some of our analyses relies on the assumption that no intermediate variable buffers the associations with IBD. For example, genetic variants and bacteria may be associated with IBD through the disease severity instead of the disease itself. In S16 Fig, we repeat the simulation from Fig 3 after replacing IBD with the disease severity drawn from a binomial distribution with n = 5. When assessing the relationships among IBD (here defined as severity ≥ 1), bacteria and genetic variants, we observed the same qualitative results for models a), b), and c) as those from our primary example. However, in model d), substantial correlation between genetics and bacteria persisted in the case sample. Nevertheless, if the severity is measured, a naïve but efficient solution consists of adjusting for that variable. To the best of our knowledge, no standard severity score is available for IBD. However, severity is commonly quantified by ineffective treatment, remission, or the number and location of inflamed sites. However, only the flare/remission status was available in our study. To assess potential confounding, we also reran the SNP-bacteria and bacteria-IBD associations. As shown in S17 Fig, we did not observe any qualitative change in the association pattern.

Patients and samples collection

All subjects were recruited in the Gastroenterology Department of the Saint Antoine Hospital (Paris, France) and provided informed consent, and approval was obtained from the local ethics committee (Comite de Protection des Personnes Ile-de-France IV, IRB 00003835, Suivitheque study). Their microbiota composition analysis has been published[44]. A diagnosis of IBD was defined by clinical, radiological, endoscopic and histological criteria. None of the study participants had taken antibiotics or used colon-cleansing products for at least 2 months prior to enrolment. Fecal samples were collected from 182 patients with IBD. Whole stools were collected in sterile boxes and immediately homogenized, and 0.2 g aliquots were frozen at −80°C for further analysis. Details of the genomic DNA extraction, 16S rRNA gene sequencing, and additional real-time quantitative PCR are presented in the S1 Text.

16S rRNA gene sequencing

Genomic DNA was extracted from 200 mg of feces as described previously[45]. Following microbial lysis involving both mechanical and chemical step, nucleic acids were precipitated by isopropanol for 10 minutes at room temperature, followed by incubation for 15 minutes on ice, and centrifugation for 30 minutes at 15,000g and 4°C. Pellets were suspended in 112 μL of phosphate buffer and 12 μL of potassium acetate. After the RNase treatment and DNA precipitation, nucleic acids were be recovered by centrifugation at 15,000g and 4°C for 30 minutes. The DNA pellet were suspended in 100 μL of TE buffer. After extraction, the total DNA concentration was measured using PicoGreen (Invitrogen), and global 16S gene DNA copy numbers were measured using a qPCR method adapted from Maeda et al[46] allowing for inhibition effect estimation and DNA concentration adjustment.

The sequence region of the 16S rRNA gene spanning the variable region V3-V5 was amplified using the broad-range forward primer For16S_519 (CAGCMGCCGCGGTAATAC) and reverse primer Rev16S_926 (CCGTCAATTCMTTTGAGTTT). Amplification reaction (initial activation step at 94°C for 1 min followed by 30 cycles of 94°C for 15 s, 43°C for 15 s and 68°C for 45 s plus final incubation at 68°C for 1 min) was performed in a total volume of 100 μL containing 1X PCR buffer, 2 mM MgSO₄, 1 U of DNA High Fidelity Taq Polymerase (Invitrogen, Carlsbad, CA), 625 nM of each barcoded primer (IDT), 250 μM of each dNTP (Invitrogen) and the concentration-adjusted DNA sample. A bidirectional library was prepared using the One Touch2 Template Kit and sequenced on PGM Ion Torrent using the Ion PGM Sequencing 400 Kit (Life Technologies, Carlsbad, CA).

The sequences were demultiplexed and quality filtered using the Quantitative Insights Into Microbial Ecology (QIIME, version 1.8.0) software package. The sequences were trimmed for barcodes and PCR primers and were binned for a minimal sequence length of 200 pb. The sequences were then assigned to Operational Taxonomic Units (OTUs) using the UCLUST algorithm[47] with a 97% threshold pairwise identity and taxonomically classified using the Greengenes reference database[48]. Rarefaction was performed (10,000 and 2,000 sequences per sample respectively for discovery and replication cohorts). Rarefied reads counts were then transpose into per-individual proportion and used to compare OTUs abundances across samples. We identified a total of 897 bacterial taxa across 6 hierarchical levels (Phylum, Class, Order, Family, Genus and Species). However 577 taxa were detected in less than 20% of the participants (S1 Fig), and were therefore removed to limit sparsity issue in further analyses.

Bacterial level pre-processing

Large-scale genomics data (e.g. microarray, proteomic, metabolomic) are commonly pre-process, applying various transformation and normalization procedures to remove confounding effects[49] but also to increase power[50]. Pre-processing can strongly impact both power and robustness[51] and should therefore be conducted carefully. To our knowledge, there is no established consensus on how to perform optimal pre-processing of microbiome data. The hierarchical relationship between taxa (Fig 1), along with strong pairwise correlation within each level (S3 Fig), makes pre-processing particularly challenging. Moreover, as we recently showed[52], standard approaches such as adjusting outcomes for the top principal components (PCs) of the taxa level (as commonly done in gene expression analysis[53]) can potentially introduce bias in our analysis–a phenomenon that would be amplified in our data because PCs can capture the shared variability of specific bacterial families. Therefore, we focused the pre-processing on transformation and filtering that will ensure the validity of the statistical tests. In particular, linear regression is sensitive to the presence of outliers, assumes a normal distribution of the outcome residual, and requires sufficient sample size for reliable effect estimation when multiple predictors are considered. The raw distributions of bacterial levels show widespread and strong outliers and skewed distribution.

In practice, we first filtered out all elements quantified in <20% of individuals or poorly annotated. We did not to use a more stringent threshold to allow for a screening as broad as possible, while we validated the most significant results through replication analysis in independent datasets. However, as showed in our sensitivity analyses (S13 Fig) we expect very low power to identify association with taxa quantified in less than 70% of the individuals in our dataset. Indeed, all top results from our assocaition analysis (Table 2) involve taxa were present in more than 85% of the individuals (e.g. Fprau was quantified in more than 95% of the individuals). Second, we applied rank-based inverse normal transformation to address both non-normal distribution of the bacteria levels and avoid false signal due to outliers. We acknowledge that such non-linear transformation can affect the inference [54]. However, rank-based inverse normal transformation is very commonly used in marginal genetic effect analyses of ‘omics data (e.g. microarray, RNA-seq, metabolites, which have properties similar to our microbiome data) and the current consensus is that, despite drawbacks, it remains a simple and efficient solution in many settings as compared to more complex approaches [55, 56]. Moreover, we also noted that rank-based inverse normal transformation is also now used in some microbiome data analysis[11, 57].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Bacteria-genetic variant associations.

https://doi.org/10.1371/journal.pgen.1008018.t002

Statistical analysis

We performed a systematic screening for association between the selected genetic variants and bacterial level. Association test were performed using standard univariate linear regression including seven potential confounding factors: gender, age, smoking, and treatment for 5-aminosalicylic acid, corticosteroids, anti-TNF agent, thiopurine or methotrexate. Importantly, we did not adjust the model for disease type (UC versus CD) as we demonstrated in a recent study that adjusting for variables that have a genetic basis might invalid association test [58]. To increase statistical power we also employed CMS (Covariates for Multiphenotype Studies), an approach recently developed by our group [19] to select additional covariates from the available bacterial data. Hence the final test we applied consisted in evaluating , the estimated effect of the genetic variant G on the quantification of a bacteria Y, in the following model: where β_Z and β_C are vectors of effect of Z, the counfounding factors, and C, a set of bacteria selected by CMS. To account for the hierarchical structure of the data, CMS search for covariates C of an outcome Y was conducted only within taxa of the same hierarchical levels as Y.

Replication study

We performed a replication of our main results using two independent studies. Both studies are themselves combination of multiple cohorts. The first study included 474 adult IBD cases from three cohorts. We only considered two of them, MGH (Mass General Hospital, N = 170) and MSH (Mont Sinai Hospital, N = 65), both from the United States. We discarded the last one because of missing data on the main confounding factors. In these two cohorts mucosa associated microbiota was analyzed from intestinal biopsies. Because of major differences in the composition and treatment between MGH and MSH (S4 Table), we analyzed the two cohorts separately and refer to them as replication 1 and replication 2. Moreover, all participants that undergone proctocolectomy with ileoanal anastomosis were excluded a priori, as the microbiota composition might be dramatically impacted by this operation and therefore not reflect the overall population. The second study included 216 IBD cases from a Netherland cohort with almost complete covariates data [21]. We refer to this dataset as replication 3. All replication analyses (1, 2, and 3) were adjusted for gender, age, smoking, mesalamine, antibiotics, and immunosuppressant treatments. To ensure our discovery signal was not biased because of the used of covariates identified by the CMS approach, we did not use CMS at the replication stage. Finally, we performed an inverse-variance meta-analysis of all three replications dataset using the R metafor package. Study-specific results are presented in S5 Table.

Evaluation of ascertainment bias

The established association between IBD and both the genetic variants and the bacterial levels implies that the joint distribution of these variables in IBD cases only might not representative of the overall population. This is a known issue that has been widely discussed in the context of secondary trait analysis [28, 29, 59–62]. Indeed, with the surge of case-control genome-wide association study, a number of investigators are facing the problem of analyzing secondary phenotypes, measured on behalf of the disease status, while accounting for the sample ascertainment. These works show that when mishandled, secondary trait analysis can be biased. In particular, if a predictor X (here, a genetic variant) and a secondary phenotype Y (here, a bacterial level) are associated with the outcome (here, IBD), the two variables might display a false association. While studies showed that for rare diseases (i.e. prevalence ≤ 1%, as for IBD), and small effect, the Y~X association test in disease cases only, as done in this study, is expected to be small [28]. Estimating the changes in expected effect estimates depending on sampling design is a non-trivial problem. As discussed in the literature for logistic models [29, 63],it has a direct relationship with the parametrization of the underlying generating models.

Consider a biallelic genetic variant g, generated from a binomial distribution with minor allele frequency p, a bacterial level b following a distribution Ω and a disease status d with probability π. For mathematical convenience we further consider D, B and G for standardized d, b and g, respectively. We are interested in estimating Δ = β_{d = 1}−β, the difference between the regression coefficient between G and B in the general population (β). Without assuming any particular distribution of the variables, Δ has the general form: where is the variance of G conditional on d = 1. As , the expectation conditional on d can we be rewritten as: So without loss of generality, the bias Δ, equals: Note that the above formulation of Δ is valid for any generating model and data distribution.

Model fit using estimates of selection bias

We are interested in assessing the impact of the bias Δ on the SNP-bacteria association for the four hypothetical causal models (Fig 3A–3D) in order to determine their fitness to our data. Here we assumed IBD status follows a binomial distribution with n = 1 and is linked to other variable through a logistic model, and we constrained the model to fit known associations: the positive association between the genetic variants and IBD [3, 4], and the negative association between the four taxa (Roseburia genus, Bacteroidia class, F. prausnitzii, and Firmicutes phylum) and IBD [8, 24, 64]. In model a) the effect of g on D is mediated by B. In model b) the effect of g on B is mediated by D. In model c), g and B act independently on IBD and are not associated in the population. Finally in model d), B and D are both influenced by U, the unmeasured risk factors.

The data can be expressed through the following generative models: Building on our estimator of Δ and leveraging previous work relating logistic and linear model [63, 65] we further demonstrate that when the disease is common, under the assumption of normality of B and G and small effects of the independent variables, Δ for the four hypothetical causal models can be approximated as (S1 Text and S8 Fig):

However, the normality assumption will not hold in general for both G (being drawn from a binomial) and B (e.g. drawn from a negative binomial), and the two other hypothesis (common disease and small effect size) are also violated in our data. Nevertheless, as previously showed [65], and as confirmed by our simulations (S1 Text and S9–S12 Figs), parameters from the linear approximation have the proper sign, so that the above Δ approximations is still informative about the direction of the bias, whatever the residual distribution of B and G. It follows that:

Thus, in model a) the effect estimate will be biased toward the null, i.e the magnitude of the genetic effect on bacteria will be underestimated. In model b), we expect the effect to be the same in case only sample as in the whole population. In model c), the ascertainment will induce a bias which direction has opposite sign to the product of the genetic effect on IBD and the bacterium effect on IBD. Finally In model d), we expect the indirect effect of the genetic variant on the bacterium to be null in a case only sample.