Bayesian bacterial GWAS model reveals threshold-dependent genetic changes in antimicrobial resistance phenotypes

Genome-wide association studies (GWASs) utilize the association between genetic variants and phenotypic changes to identify causal variants underlying a phenotype of interest, and in some cases employ these correlations to perform phenotype prediction for novel genomes given their repertoire of genetic variants. Often utilizing a logistic or linear regression framework, bacterial GWAS methods have frequently been applied to infer variants associated with antimicrobial resistance (AMR) phenotypes measured using serial dilutions and recorded as minimum inhibitory concentrations (MICs). Beyond their increased granularity compared to resistant-sensitive classifications, MICs provide a more robust measure of AMR because they are not dependent on threshold values that can change over time and between data sources. One prominent challenge posed by both existing inference and prediction GWAS frameworks is the representation of MICs as censored ordered intervals rather than as binary resistant-sensitive or continuous approximations, a more suitable statistical model that has the potential to improve the discovery of lower-penetrance variants and MIC prediction. We show that Bayesian ordered logistic regression GWAS models can recover known antimicrobial resistance variants, identify variants with threshold-dependent effects without employing computationally-costly pairwise comparisons, and predict categorical MIC phenotypes. Specifically, the models recovered known resistance-conferring variants for {beta}-lactams (in the penicillin-binding protein genes pbp2x, pbp1a, and pbp2b) and trimethoprim (in the dihydrofolate reductase gene folA) in S. pneumoniae, and for rifampicin (in the rpoB gene encoding the beta subunit of bacterial RNA polymerase) for M. tuberculosis, in addition to numerous lower-penetrance resistance variants. We found that ordered categorical models that allow for non-proportionality in the effect of each genetic variant across different MIC categories replicate the pattern of stepwise accumulation of resistance mutations in pbp genes in their relative effects at different MIC thresholds, and can improve balanced accuracy for MIC prediction compared to a proportional-effects model, suggesting that nonlinear effects may play a role in missing heritability. We also describe the use of lineage clustering as an alternative method of reducing false positives induced by genetic relatedness within a bacterial population. The methods described here incorporate a more accurate representation of minimum inhibitory concentration phenotypes into a highly flexible bacterial GWAS methodology capable of capturing nonlinear variant effects and suitable for both causal variant inference and prediction of MIC phenotypes.