CALDERA: finding all significant de Bruijn subgraphs for bacterial GWAS

Abstract

Motivation: Genome-wide association studies (GWAS), aiming to find genetic variants associated with a trait, have widely been used on bacteria to identify genetic determinants of drug resistance or hypervirulence. Recent bacterial GWAS methods usually rely on k-mers, whose presence in a genome can denote variants ranging from single-nucleotide polymorphisms to mobile genetic elements. This approach does not require a reference genome, making it easier to account for accessory genes. However, a same gene can exist in slightly different versions across different strains, leading to diluted effects.

Results: Here, we overcome this issue by testing covariates built from closed connected subgraphs (CCSs) of the de Bruijn graph defined over genomic k-mers. These covariates capture polymorphic genes as a single entity, improving k-mer-based GWAS both in terms of power and interpretability. However, a method naively testing all possible subgraphs would be powerless due to multiple testing corrections, and the mere exploration of these subgraphs would quickly become computationally intractable. The concept of testable hypothesis has successfully been used to address both problems in similar contexts. We leverage this concept to test all CCSs by proposing a novel enumeration scheme for these objects which fully exploits the pruning opportunity offered by testability, resulting in drastic improvements in computational efficiency. Our method integrates with existing visual tools to facilitate interpretation.

Availability and implementation: We provide an implementation of our method, as well as code to reproduce all results at https://github.com/HectorRDB/Caldera_ISMB.

Supplementary information: Supplementary data are available at Bioinformatics online.

Fig. 1.

Example of de Bruijn graphs. (a) A general example with two genes, each with some variability, resulting in a mostly linear sequence only at the coarse level. More details in Section 4. (b) A simpler setting with two samples and four nodes, leading to three CCSs: $\left\{v_1\right\},\hspace{0.17em}\left\{v_2\right\}$ and $\{v_0,v_1,v_2,v_3\}$

Fig. 2.

A short illustration of CALDERA’s exploration. A simple graph with four nodes and four samples. Nodes v₁, v₂ and v₃ are linked to node v₄. To construct the CCS $\{v_1,v_2,v_3,v_4\}$ from $\{v_3,v_4\}$, we can either first add v₁—and thus construct $\{v_1,v_3,v_4\}$—and then add v₂, or directly add v₂ and then we get $\{v_1,v_2,v_3,v_4\}$ by closure. To avoid enumerating $\{v_1,v_2,v_3,v_4\}$ twice, we therefore need local itemtables

Fig. 3.

Results of CALDERA. (a) Run times for CALDERA and COIN+LAMP on graphs with various values of p. In this setting, n = 100. (b) Proportion of all unitigs associated with the resistant phenotype that are found to be significant by CALDERA, the LAMP2 procedure on all unitigs and DBGWAS, as the value of α changes

Fig. 4.

Akkermansia dataset: Tubulin/FtsZ_GTPase gene. Screenshot from the output of CALDERA. We select the first component, which is the one which contains the most significant CCS. (a) Unitigs belonging to the most significant CCS are colored in green (darker gray in the black and white version of the document). Other unitigs linked to the CCS in the DBG are colored in gray. Size denotes overall frequency, while a black contour denotes that the sequence of the unitig has a match in the database, here RefSeq (Tatusova et al., 2016). (b) All significant hits on that database are listed in a panel, usually on top of the subgraph. (c) User can click on nodes to display information, or right-click to select all nodes from the same component. This contains info on the node, such as the frequency, the pattern of the associated CCS, or any match to the database (A color version of this figure appears in the online version of this article)