# Predicting susceptibility to tuberculosis based on gene expression profiling in dendritic cells

John D Blischak, Ludovic Tailleux, Marsha Myrthil, Cécile Charlois, Emmanuel Bergot, Aurélien Dinh, Gloria Morizot, Olivia Chény, Cassandre Von Platen, Jean-Louis Herrmann, et al.

► To cite this version:

John D Blischak, Ludovic Tailleux, Marsha Myrthil, Cécile Charlois, Emmanuel Bergot, et al. Predicting susceptibility to tuberculosis based on gene expression profiling in dendritic cells. Scientific Reports, 2017, 7 (1), pp.5702. 10.1038/s41598-017-05878-w. pasteur-02046440

HAL Id: pasteur-02046440

https://pasteur.hal.science/pasteur-02046440

Submitted on 22 Feb 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_2.png" /></figure>

OPEN

# Predicting susceptibility to tuberculosis based on gene expression profiling in dendritic cells

Received: 13 February 2017

Accepted: 5 June 2017

Published online: 18 July 2017

John D. Blischak<sup>1,2</sup>, Ludovic Tailleux<sup>3</sup>, Marsha Myrthil<sup>1</sup>, Cécile Charlois<sup>4</sup>, Emmanuel Bergot<sup>5</sup>, Aurélien Dinh<sup>6</sup>, Gloria Morizot<sup>7</sup>, Olivia Chény<sup>8</sup>, Cassandre Von Platen<sup>8</sup>, Jean-Louis Herrmann<sup>9,10</sup>, Roland Brosch<sup>1,3</sup>, Luis B. Barreiro<sup>11,12</sup> & Yoav Gilad<sup>1,13</sup>

Tuberculosis (TB) is a deadly infectious disease, which kills millions of people every year. The causative pathogen, *Mycobacterium tuberculosis* (MTB), is estimated to have infected up to a third of the world's population; however, only approximately 10% of infected healthy individuals progress to active TB. Despite evidence for heritability, it is not currently possible to predict who may develop TB. To explore approaches to classify susceptibility to TB, we infected with MTB dendritic cells (DCs) from putatively resistant individuals diagnosed with latent TB, and from susceptible individuals that had recovered from active TB. We measured gene expression levels in infected and non-infected cells and found hundreds of differentially expressed genes between susceptible and resistant individuals in the non-infected cells. We further found that genetic polymorphisms nearby the differentially expressed genes between susceptible and resistant individuals are more likely to be associated with TB susceptibility in published GWAS data. Lastly, we trained a classifier based on the gene expression levels in the non-infected cells, and demonstrated reasonable performance on our data and an independent data set. Overall, our promising results from this small study suggest that training a classifier on a larger cohort may enable us to accurately predict TB susceptibility.

Tuberculosis (TB) is a major public health issue. Worldwide, over a million people die of TB annually, and millions more currently live with the disease<sup>1-3</sup>. Successful treatment requires months of antibiotic therapy<sup>4</sup>, and drug-resistant strains of *Mycobacterium tuberculosis* (MTB) continuously emerge<sup>5</sup>. Approximately a third of the world's population is estimated to be infected with MTB, but most are asymptomatic. While these naturally resistant individuals are able to avoid active disease, MTB might persist in a dormant state, known as latent TB<sup>6</sup>. In contrast, approximately 10% of individuals will develop active TB after infection with MTB<sup>7,8</sup>. Unfortunately, we are currently unable to predict if an individual is susceptible. While twin and family studies have indicated a heritable component of TB susceptibility<sup>9-12</sup>, genome wide association studies (GWAS) have only identified a few loci with low effect size<sup>13-19</sup>. Due to the highly polygenic architecture, it may be informative to examine differences between susceptible and resistant individuals at a higher level of organization, e.g. gene regulatory networks.

<footnote n="1"><sup>1</sup>Department of Human Genetics, University of Chicago, Chicago, Illinois, USA. <sup>2</sup>Committee on Genetics, Genomics, and Systems Biology, University of Chicago, Chicago, Illinois, USA. <sup>3</sup>Integrated Mycobacterial Pathogenomics, Institut Pasteur, Paris, France. <sup>4</sup>Centre de Lutte Antituberculeuse de Paris, DASES Mairie de Paris, 75013, Paris, France. <sup>5</sup>Service de pneumologie et oncologie thoracique, CHU Côte de Nacre, 14033, Caen, France. <sup>6</sup>Maladies Infectieuses, AP-HP, Hôpital Universitaire Raymond-Poincaré, Garches, 92380, France. <sup>7</sup>Clinical Investigation & Access Biological Resources (ICAReB), Institut Pasteur, Paris, France. <sup>8</sup>Clinical Core, Centre for Translational Science, Institut Pasteur, Paris, France. <sup>9</sup>INSERM, U1173, UFR Simone Veil, Université de Versailles Saint Quentin, Saint Quentin en Yvelines, France. <sup>10</sup>APHP, Groupe Hospitalo-Universitaire Paris Île-de-France Ouest, Garches et Boulogne-Billancourt, France. <sup>11</sup>Department of Genetics, CHU Sainte-Justine Research Center, Montreal, Québec, Canada. <sup>12</sup>Department of Pediatrics, University of Montreal, Montreal, Québec, Canada. <sup>13</sup>Department of Medicine, University of Chicago, Chicago, Illinois, USA. John D. Blischak and Ludovic Tailleux contributed equally to this work. Correspondence and requests for materials should be addressed to L.T. (email: tailleux@pasteur.fr) or L.B.B. (email: luis.barreiro@umontreal.ca) or Y.G. (email: gilad@uchicago.edu)</footnote>
Using this approach, previous studies have characterized gene expression profiles in innate immune cells isolated from individuals known to be susceptible or resistant to infectious diseases, including those with latent or active TB²⁰ and acute rheumatic fever²¹.

We hypothesized that gene expression profiles in innate immune cells may be used to classify individuals with respect to their susceptibility to develop active TB. To test this hypothesis, we differentiated dendritic cells (DCs) from monocytes isolated from individuals that had recovered from a past episode of active TB, which we refer to as susceptible, and from individuals with confirmed latent TB, which we refer to as putatively resistant (this group is enriched in resistant individuals but we cannot exclude that some still have the potential to develop active TB²²). We infected the DCs with MTB because these innate immune cells help shape the adaptive immune response, which is critical for fighting MTB²³, ²⁴. We discovered that the gene expression differences between innate immune cells from resistant and susceptible individuals were present primarily in the non-infected state, that these differentially expressed genes were enriched for nearby SNPs with low p-values in TB susceptibility GWAS, and furthermore, that these gene expression levels could be used to classify individuals based on their susceptibility status.

## Results

**Susceptible individuals have an altered transcriptome in the non-infected state.** We obtained whole blood samples from 25 healthy male Caucasian individuals (Supplementary Data S1). Six of the donors had recovered from active TB, and are thus putatively susceptible. The remaining 19 tested positive for latent TB without ever experiencing symptoms of active TB, and are thus putatively resistant. We isolated dendritic cells (DCs) and treated them with *Mycobacterium tuberculosis* (MTB) or a mock control for 18 hours. To measure genome-wide gene expression levels in infected and non-infected samples, we isolated and sequenced RNA using a processing pipeline designed to minimize the introduction of unwanted technical variation (Supplementary Fig. S1). We obtained a mean (±SEM) of 48 ± 6 million raw reads per sample. We performed quality control analyses to remove non-expressed genes (Supplementary Fig. S2; Supplementary Data S2), identify and remove outliers (Supplementary Figs S3, S4 and S5), and check for confounding batch effects (Supplementary Figs S6 and S7). Ultimately, data from 6 of the 50 samples (25 individuals × 2 treatments) failed the quality checks and were removed from all downstream analyses (Supplementary Fig. S5). The remaining 44 high-quality samples included 17 putatively resistant individuals and 5 putatively susceptible individuals.

We performed a standard differential expression analysis using a linear modeling framework (Supplementary Data S3), defined in equation (1). As expected, there was a strong response to infection with MTB in both resistant and susceptible individuals (Supplementary Fig. S8). Considering the putatively resistant individuals, we identified 3,486 differentially expressed (DE) genes between the non-infected and infected states at a q-value of 10% and an arbitrary absolute log-fold change greater than 1. Similarly, 3,789 genes were classified as DE between the non-infected and infected states in the putatively susceptible individuals. In both classes of samples, the DE genes included the important immune response factors *IL12B*, *REL*, and *TNF*. While the treatment effect was obvious in all individuals, of most interest were the patterns of gene expression differences between the putatively susceptible and resistant individuals in either the non-infected or infected states (Fig. 1). We identified 645 DE genes between putatively resistant and susceptible individuals in the non-infected state at a q-value of 10%, including *ATPV1B2*, *FEZ2*, *PSMA2*, *TNFRSF25*, and *TRIM38*. The log₂ fold change in gene expression for these 645 genes was relatively small (less than 2), which was expected given that susceptibility to TB is a polygenic trait influenced by many genes. In contrast, no genes were DE between putatively resistant and susceptible individuals in the infected state (at a q-value of 10%).

**Differentially expressed genes are enriched with TB susceptibility loci.** We next sought evidence that genes classified as DE in our *in vitro* experimental system play a role in determining susceptibility to TB. To do this, we intersected our data with results from TB susceptibility GWAS conducted in Russia¹⁸, The Gambia¹³, Ghana¹³, and Uganda and Tanzania¹⁹. We also included data from a height GWAS conducted in individuals of European ancestry²⁵ as a negative control. To perform a combined analysis of our gene expression data and the GWAS results, we had to define pairs of genes (for which we have expression data) and SNPs (for which we obtained GWAS *P* values). Thus, each gene in our expression data was coupled with the GWAS SNP with the lowest *p*-value among all tested SNPs located within 50 kb of the gene's transcription start site (Supplementary Data S4; this gene-SNP definition was performed separately for each GWAS data set).

Once we defined gene-SNP pairs, we asked whether differences in gene expression levels between putatively resistant and susceptible individuals could help us identify genetic variation that is associated with susceptibility to TB. In other words, we asked whether increasing evidence for DE genes is associated with low GWAS *p*-values. To do so, we calculated the fraction of SNPs with a GWAS *p*-value lower than 0.05 among SNPs that were paired with ranked subsets of genes whose expression profiles show increasing effect size of expression differences between putatively resistant and susceptible individuals. In order to assess the significance of the observations, we performed 100 permutations of the enrichment analysis to derive an empirical *p*-value.

Using this approach, we observed a clear enrichment (empirical *P* < 0.01) of low *p*-values for TB GWAS SNPs that are paired with genes that are differentially expressed between susceptible and resistant individuals in the non-infected state (Fig. 2a). In fact, we observed significant enrichments of lower GWAS *p* values (empirical *P* < 0.01) in all 4 TB susceptibility GWAS (Russia, The Gambia, Ghana, Uganda and Tanzania) (Supplementary Fig. S10) for all 4 differential expression contrasts, namely resistant vs. susceptible individuals in the non-infected state (Fig. 2c), resistant vs. susceptible individuals in the infected state (Fig. 2d), effect of treatment in resistant individuals (Fig. 2e), and effect of treatment in susceptible individuals (Fig. 2f). Reassuringly, we did not observe an enrichment of low *p* values (empirical *P* > 0.01) when we used the same approach to consider data from the height GWAS (Fig. 2b–f).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Results of differential expression analysis. The top panels show the distributions of unadjusted p-values for testing the null of no differential expression between susceptible and resistant individuals in the (a) non-infected or (b) infected state. The bottom panels show the corresponding volcano plots for the (c) non-infected and (d) infected states. The x-axis is the log fold change in gene expression level between susceptible and resistant individuals and the y-axis is the −log₁₀ p-value. Red indicates genes that are classified as differentially expressed with a q-value less than 10%.</figcaption></figure>

**Susceptibility status can be predicted based on gene expression data.** Next we attempted to build a gene expression-based classifier to predict TB susceptibility status (Supplementary Data S5). We focused on the gene expression levels measured in the non-infected state both because this is where we observed the largest gene regulatory differences between putatively susceptible and resistant individuals (Fig. 1a,c), and also because, from the perspective of an ultimate translational application, it is more practical to obtain gene expression data from non-infected DCs. We trained a support vector machine using the 99 genes that were differentially expressed between resistant and susceptible individuals in the non-infected state at a q-value less than 5% (see Methods for a full description of how we selected this model). Encouragingly, we observed a clear separation between putatively susceptible and resistant individuals when comparing the predicted probability of being susceptible to TB for each sample obtained from leave-one-out-cross-validation (Fig. 3a). Using a cutoff of 0.25 for the predicted probability of being susceptible to TB, we obtained a sensitivity of 100% (5 out of 5 susceptible individuals classified as susceptible), a specificity of 88% (15 out of 17 resistant individuals classified as resistant), and a positive predictive value (PPV) of 71% (5 of 7 individuals classified as susceptible were susceptible).

Unfortunately our current data set is too small to properly split into separate training and testing sets (it is challenging to collect samples from previous TB patients, who are healthy and have no medical reason to go back for a GP visit). To our knowledge, there are also no other suitable data sets available with which to test out classifier (that said, see Supplementary Fig. S17 for the results of applying the classifier to a non-ideal data set, which measured gene expression in macrophages from a small number of individuals²⁰). Thus, in order to further assess the plausibility of our model, we applied the classifier to data from an independent study, which collected
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

**Figure 2.** Comparison of differential expression and TB susceptibility GWAS results. (a,b) The y-axis is the fold enrichment of SNPs with p-value less than 0.05 from the (a) GWAS of TB susceptibility in Russia¹⁸ or (b) height in individuals of European descent²⁵. The x-axis is bins of genes with increasingly stringent effect size cutoffs of the absolute expression log fold change between putatively susceptible and resistant individuals in the non-infected state. The effect size cutoffs were chosen such that each bin from left to right contained approximately 25 fewer genes. The red line shows the results from the actual data. The grey lines are the results from 100 permutations. The dashed blue line at y = 1 represents the null expectation. (c–f) Boxplots of the area under the curve of the fold enrichment (red line in (a,b)) minus the background level (blue y = 1 line in (a,b)) for each of the 5 GWAS¹³, ¹⁸, ¹⁹ considered for the 4 differential expression contrasts: (c) resistant vs. susceptible individuals in the non-infected state, (d) resistant vs. susceptible individuals in the infected state, (e) effect of treatment in resistant individuals, (f) effect of treatment in susceptible individuals. The boxplot is the result of the 100 permutations, and the red point is the result from the actual data. As a reference, the leftmost boxplot in (c) corresponds to the enrichment plot in (a), and the rightmost boxplot in (c) corresponds to the enrichment plot in (b).

genome-wide gene expression levels in DCs from 65 healthy individuals²⁴, none with a previous history of TB. Using the cutoff of 0.25 for the probability of being susceptible to TB (determined to be optimal in the training set), 11% (7 of 65) of the individuals were classified as susceptible to TB (Fig. 3b). Adjusting for the PPV obtained from the training set (71%), our model predicted that 7.7% of the healthy individuals were susceptible. While we cannot confirm this result (the true susceptibility status of these 65 individuals is unknown), this observation is encouraging because our estimate is similar to the commonly used inference that roughly 10% of the general population is susceptible to progression of infection from latent to active TB.
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Classifying TB susceptible individuals using a support vector machine model. (a) The estimates of predicted probability of TB susceptibility from the leave-one-out-cross-validation for individuals in the current study. The blue circles represent individuals known to be susceptible to TB, and orange those resistant to TB. The horizontal dashed red line at a probability of 0.25 separates susceptible and resistant individuals. (b) The estimates of predicted probability of TB susceptibility from applying the classifier trained on the data from the current study to a test set of independently collected healthy individuals²⁴.</figcaption></figure>

Discussion

We obtained dendritic cells (DCs) from individuals that were known to be putatively susceptible or resistant to developing active tuberculosis (TB) and measured genome-wide gene expression levels in non-infected DCs and DCs infected with *Mycobacterium tuberculosis* (MTB) for 18 hours. As expected, there were large changes in gene expression due to MTB infection in the DCs from both putatively resistant and susceptible individuals (Supplementary Fig. S8). We identified 645 genes which were differentially expressed (DE) between susceptible and resistant individuals in the non-infected state; whereas, we did not observe any DE genes between susceptible and resistant individuals in the infected state (Fig. 1). This suggests that the differences in the transcriptomes between DCs of resistant and susceptible individuals are present pre-infection. Yet, 18 hours after infection gene expression profiles in both susceptible and resistant individuals have converged to a similar gene regulatory network, presumably to fight the infection. We confirmed that the absence of DE genes in the infected state is not caused by a decrease in statistical power due to an overall increase in gene expression variance upon infection (Supplementary Fig. S9). We chose to measure gene expression 18 hours post-infection because this time point was previously associated with a large change in genome-wide gene expression levels²⁶. Given our observations, however, future studies investigating the difference in the innate immune response between individuals resistant and susceptible to TB may want to focus on earlier time points post-infection.

It is important to note that our study was not designed to uncover the mechanisms underlying susceptibility or resistance to TB, but to try and find a gene regulatory signature that might allow us to classify individuals as either susceptible or resistant. That said, among the 645 DE genes between resistant and susceptible individuals in the non-infected state, there were many interesting genes involved in important innate immune activities critical for fighting MTB and other pathogens such as autophagy²⁷, ²⁸, phagolysosomal acidification, and antigen processing. In particular, *FEZ2*, a suppressor of autophagosome formation²⁹, was down-regulated when DCs were infected with MTB; however, in the non-infected DCs, this gene has elevated expression level in susceptible compared with resistant individuals. In turn, *ATP6V1B2*, a gene coding for a subunit of the proton transporter responsible for acidifying phagolysosomes³⁰⁻³², has increased expression in susceptible individuals compared to resistant in the non-infected state. Lastly, genes coding for nine subunits of the proteasome, which is critical for processing of MTB antigens to be presented via major histocompatibility complex (MHC) class I molecules³³⁻³⁶, have increased expression in susceptible individuals compared to resistant in the non-infected state. These genes are candidates for future functional studies investigating the mechanisms of TB susceptibility.

We observed that DE genes in our *in vitro* experimental system were enriched for lower GWAS p-values (Fig. 2). This suggests that such *in vitro* approaches are informative for interrogating the genetic basis of disease susceptibility. That being said, we recognize a major caveat with this analysis is that assigning SNPs to their nearest gene on the linear chromosome is problematic because regulatory variants can have longer range effects. Nevertheless, considering this limitation, it was encouraging that we were able to detect evidence of the genetic basis of TB susceptibility in this system.

Not only did this analysis identify a global enrichment of TB susceptibility loci, but by intersecting the expression and GWAS data, we were able to identify a few interesting candidate genes, which were only marginally significant in the differential expression analysis (Supplementary Data S3) and in the original GWAS (Supplementary Data S4). Here we highlight two genes (*CCL1* and *UNC13A*), which have been previously
implicated in the response to MTB infection. CCL1 is a chemokine that stimulates migration of monocytes³⁷. In our study, it was upregulated in susceptible individuals compared to resistant in both the non-infected and infected states (but did not reach statistical significance in either) and was statistically significantly upregulated with MTB treatment. Furthermore, the nearby SNP assigned to CCL1 had a p-value less than 0.01 in the TB susceptibility GWAS from The Gambia and Ghana. A previous differential expression study of TB susceptibility (discussed in more detail below) found that CCL1 was upregulated to a greater extent 4 hours post-infection with MTB in macrophages isolated from individuals with active TB (i.e. susceptible) compared to individuals with latent TB (i.e. resistant)²⁰. Additionally they performed a candidate gene association study and found that SNPs nearby CCL1 were associated with TB susceptibility. In our previous study, we discovered that CCL1 was one of only 288 genes that were differentially expressed in macrophages 48 hours post-infection with MTB and related mycobacterial species but not unrelated virulent bacteria³⁸. UNC13A is involved in vesicle formation³⁹. In our study, it was downregulated in susceptible individuals compared to resistant in both the non-infected and infected states (but did not reach statistical significance in either) and was statistically significantly upregulated with MTB treatment. Furthermore, the nearby SNP assigned to UNC13A had a p-value less than 0.01 in the TB susceptibility GWAS from Russia, The Gambia, and Ghana. In our past study mapping expression quantitative trait loci (eQTLs) in DCs 18 hours post-infection with MTB, UNC13A was one of only 98 genes which were associated with an eQTL post-infection but not pre-infection, which we called MTB-specific eQTLs²⁴. Thus our new results increased the evidence that CCL1 and UNC13A may potentially impact TB susceptibility.

Previous attempts to use gene expression based classifiers in the context of TB have focused on predicting the status of an infection rather than the susceptibility status of an individual⁸, ⁴⁰, ⁴¹. In other words, the goal of most previous studies was to detect individuals in the early stages of active TB when antibiotic intervention would be most effective or to monitor the effectiveness of a treatment regimen⁴². In contrast, our goal was not to distinguish between active and latent TB, but instead to be able to determine susceptibility status before individuals are infected with MTB. Even with our small sample size, we were able to successfully train a classier with high sensitivity and reasonable specificity. Because such a classification of susceptibility status could affect the decision of whether or not to take antibiotics to treat latent TB⁶, false negatives (susceptible individuals mistakenly classified as resistant) would be much more harmful than false positives (resistant individuals mistakenly classified as susceptible). For that reason, we emphasized sensitivity over specificity.

While the small sample size was the main concern in interpreting the classifier results, we recognize multiple other potential caveats of our study design. First, because the cells used in the experiments were collected from Caucasians, it is unknown if the classifier would perform well when applied to other human populations. Second, because the classifier focused solely on the transcriptional response of the host immune cell to infection with MTB H37Rv, it may not be suitable for the detection of differences in susceptibility with respect to other, more or less virulent MTB strains⁴³. Third, because the susceptible individuals had already experienced active TB, it is possible that their innate immune cells retained memory of the infection⁴⁴, ⁴⁵. If this innate immune memory affected how the monocytes differentiated into DCs in response to cytokine stimulation, or in general affected gene regulation, this would drastically change the interpretation of our results.

With respect to the third caveat, while we cannot entirely exclude the possibility of trained immunity, we do not believe there is sufficient evidence to support this possibility in our case. Innate immune memory is known to be a short-term phenomenon⁴⁵, its affects may be erased when immune cells are moved to a new microenvironment⁴⁶, and the whole blood transcriptional signature of active TB disappears 6–12 months after the initiation of treatment⁴⁰, ⁴⁷. Moreover, the GWAS enrichment results (Fig. 2) suggest that trained immunity cannot explain our observations. Indeed, the differences in gene expression between susceptible and resistant individuals in the non-infected state were enriched for low p-values from TB susceptibility GWAS. This observation cannot be explained by possible downstream consequences of experiencing active TB because enrichment of DE genes in GWAS data supports a causal relationship. In other words, if most of the DE genes we observed could be explained by trained immunity, one would not expect these genes to be associated with susceptibility to the disease. That said, we recognize that the ideal study would be a long-term prospective design in which blood samples were obtained from a large cohort of individuals prior to MTB infection.

To our knowledge, our study was only the second to collect data from *in vitro* MTB-infected innate immune cells isolated from individuals known to be putatively susceptible to MTB²⁰. However, there were substantial differences between our study and that of Thuong *et al*.²⁰. First, they derived and infected macrophages, the primary target host cell in which MTB resides; whereas, we derived and infected DCs, which play a larger role in stimulating the adaptive immune response to MTB. Second, we collected samples from a larger number of putatively resistant individuals (19 versus 4), increasing our power to distinguish between the gene expression profiles of susceptible and resistant individuals. Third, they measured gene expression with microarrays; whereas, we used RNA-sequencing. Considering the substantial technical differences between the methods used and the biological differences between DCs and macrophages²⁶, ⁴⁸, unsurprisingly, we were unable to identify the susceptible individuals from Thuong *et al*.²⁰ using our classifier (Supplementary Fig. S17).

Indeed, at this time, we are not aware of any other data set from healthy individuals known to be susceptible to TB, with which we can further test our classifier. When we applied our classifier to an independent set of non-infected DCs isolated from healthy individuals of unknown susceptibility status, our model predicted that 7.7–11% of the individuals were susceptible to TB, which reassuringly is similar to the average in the general population (10%). Despite this, our results must be interpreted cautiously; at best as a proof-of-principle, due to our very small sample size of only 5 susceptible individuals. That said, our promising results in this small study suggest that collecting blood samples from a larger cohort of susceptible individuals would enable building a gene expression based classifier able to confidently assess risk of TB susceptibility. By reducing the number of resistant individuals receiving treatment for latent TB, we can eliminate the adverse health effects of a 6 month regimen of antibiotics for these individuals and also reduce the selective pressures on MTB to develop drug resistance.
## Methods

**Ethics statement.** We recruited 25 subjects to donate a blood sample for use in our study. All methods were carried out in accordance with relevant guidelines and regulations. All participants gave written informed consent in accordance with the Declaration of Helsinki principles. Peripheral human blood was collected from patients at ICAREB platform of Institut Pasteur Paris and at the Centre for Infectious Disease Prevention, University hospital Caen. The Protocol has been approved by French Ethical Committee (CPP North Ouest III, no. A12 - D33 -VOL.13), and by the Institutional Review Boards of the University of Chicago (10-504-B) and the Institut Pasteur (IRB00006966).

**Sample collection.** We collected whole blood samples from healthy, HIV-negative Caucasian male individuals living in France. The putatively resistant individuals tested positive for latent TB in an interferon-γ release assay, but had never developed active TB. The putatively susceptible individuals had developed active TB in the past, but were currently healthy after having completed at least 6 months of antibiotic chemotherapy. All putatively susceptible individuals had finished their treatment at least 6 months prior to sample collection and were considered fully recovered based on a chest X-ray, sputum culture, and lack of clinical symptoms. We do not know which strain of MTB had infected any of the individuals in the study.

**Isolation and infection of dendritic cells.** We performed these experiments as previously described²⁴. Briefly, we isolated mononuclear cells from the whole blood samples using Ficoll-Paque centrifugation, extracted monocytes via CD14 positive selection, and differentiated the monocytes into dendritic cells (DCs) by culturing them for 5 days in RPMI 1640 (Invitrogen) supplemented with 10% heat-inactivated FCS (Dutscher), L-glutamine (Invitrogen), GM-CSF (20 ng/mL; Immunotools), and IL-4 (20 ng/mL; Immunotools). Next we infected the DCs with *Mycobacterium tuberculosis* (MTB) H37Rv at a multiplicity of infection of 1-to-1 for 18 hours.

**RNA extraction and sequencing.** We extracted RNA using the Qiagen miRNeasy Kit and prepared sequencing libraries using the Illumina TruSeq Kit. We sent the master mixes to the University of Chicago Functional Genomics Facility to be sequenced on an Illumina HiSeq 4000. We designed the batches for RNA extraction, library preparation, and sequencing to balance the experimental factors of interest and thus avoid potential technical confounders (Supplementary Fig. S1).

**Read mapping.** We mapped reads to human genome hg38 (GRCh38) using Subread⁴⁹ and discarded non-uniquely mapping reads. We downloaded the exon coordinates of 19,800 Ensembl³⁰ protein-coding genes (Ensembl 83, Dec 2015, GRCh38.p5) using the R/Bioconductor⁵¹ package biomaRt⁵²,⁵³ and assigned mapped reads to these genes using featureCounts⁵⁴.

**Quality control.** First we filtered genes based on their expression level by removing all genes with a transformed median log₂ counts per million (cpm) of less than zero. This step resulted in a set of 11,336 genes for downstream analysis (Supplementary Fig. S2, Supplementary Data S2). Next we used principal components analysis (PCA) and hierarchical clustering to identify and remove 6 outlier samples (Supplementary Figs S3, S4 and S5). We did this systematically, by removing any sample whose data projections did not fall within two standard deviations of the mean for any of the first six PCs (for the first PC, which separated the samples by treatment, we calculated a separate mean for the non-infected and infected samples). The remaining 44 high-quality samples included 17 putatively resistant individuals and 5 putatively susceptible individuals in both the non-infected and infected states (Supplementary Fig. S5). Note that after outlier removal, the study design was no longer completely paired because a subset of individuals was only represented by available data from either the non-infected or the infected state.

After filtering lowly expressed genes and removing outliers, we performed the PCA again to check for any potential confounding technical batch effects (Supplementary Fig. S6). Reassuringly, the major sources of variation in the data were from the biological factors of interest. PC1 was strongly correlated with the effect of treatment, and PCs 2–6 were correlated with inter-individual variation. The only concerning technical factor was the infection experiments, which were done in 12 separate batches (Supplementary Fig. S1). Infection batch correlated with PCs 3 and 5; however, we verified that this variation was not confounded with our primary outcome of interest, TB susceptibility (Supplementary Fig. S7).

**Differential expression analysis.** We used limma+voom⁵⁵⁻⁵⁷ to implement the following linear model to test for differential expression:

$$Y \sim \beta_0 + X_{treat} \beta_{treat} + X_{status} \beta_{status} + X_{treat, status} \beta_{treat, status} + I + \varepsilon \quad (1)$$

where $β_0$ is the mean expression level in non-infected cells of resistant individuals, $β_{treat}$ is the fixed effect of treatment in resistant individuals, $β_{status}$ is the fixed effect of susceptibility status in non-infected cells, $β_{treat, status}$ is the fixed interaction effect of treatment in susceptible individuals (i.e. modeling the interaction between treatment and susceptibility status), and $I$ is the random effect of individual. The random individual effect was implemented using the limma function duplicateCorrelation⁵⁸. To jointly model the data with voom and duplicateCorrelation, we followed the recommended best practice of running both voom and duplicateCorrelation twice in succession⁵⁹.

We used the model to test different hypotheses (Supplementary Data S3). We identified genes which were differentially expressed (DE) between infected and non-infected DCs of resistant individuals by testing $β_{treat} = 0$, genes which were DE between infected and non-infected DCs of susceptible individuals by testing $β_{treat} + β_{treat, status} = 0$, genes which were DE between susceptible and resistant individuals in the non-infected
state by testing $\beta_{status} = 0$, and genes which were DE between susceptible and resistant individuals in the infected state by testing $\beta_{status} + \beta_{treat status} = 0$. We corrected for multiple testing using q-values estimated via adaptive shrinkage⁶⁰ and considered differentially expressed genes as those with a q-value less than 10%.

Note that we also tested the interaction term, $\beta_{treat, status} = 0$, to identify genes in which the difference in expression level between the infected and non-infected states was significantly different between susceptible and resistant individuals. However, as expected since no DE genes were identified between susceptible and resistant individuals in the infected state (see Results), the results of testing the interaction term were partially redundant with the results of testing differences between susceptible and resistant individuals in the non-infected state, and thus we ignored these results throughout this study.

**Combined analysis of gene expression data and GWAS results.** The GWAS p-values were from previously published studies of TB susceptibility conducted in Russia¹⁸, The Gambia¹³, Ghana¹³, and Uganda and Tanzania¹⁹ (and a height GWAS in individuals of European descent²⁵). To perform a combined analysis of the gene expression and the summary statistics from each GWAS, we assigned each gene to the SNP with the minimum GWAS p-value out of all the SNPs located within 50 kb up or downstream of its transcription start site. Specifically, we obtained the genomic coordinates of the SNPs with the R/Bioconductor⁵¹ package SNPlocs. Hsapiens.dbSNP144.GRCh38 and matched SNPs to nearby genes using GenomicRanges⁶¹. 10,265 to 11,060 of the 11,336 genes were assigned an association p-value depending on the GWAS (Supplementary Data S4). For each of the 4 differential expression contrasts we tested (resistant vs. susceptible individuals in the non-infected state, resistant vs. susceptible individuals in the infected state, effect of treatment in resistant individuals, effect of treatment in susceptible individuals), we also performed an enrichment analysis. To do so, we calculated the fraction of genes assigned a GWAS SNP with p-value less than 0.05 for bins of genes filtered by increasingly stringent cutoffs for the observed differential expression effect size (the absolute value of the log fold change). The effect size cutoffs were chosen such that on average each subsequent bin differed by 25 genes. To measure enrichment, we calculated the area under the curve using the R package flux⁶² (we also subtracted the background area under the line y=1 because the number of genes assigned a SNP varied across the GWAS). In order to assess significance, we calculated the area under the curve for 100 permutations of the data. All differential expression tests were statistically significantly enriched for SNPs with low GWAS p-values in every TB susceptibility GWAS (empirical P<0.01) and not enriched for the height GWAS (empirical P>0.01) (Fig. 2; Supplementary Fig. S10).

**Classifier.** The training set included data from the 22 high-quality non-infected samples from this study with known susceptibility status. The test set included the 65 non-infected samples from one of our previous studies in which the susceptibility status is unknown²⁴, and thus assumed to be similar to that in the general population (10%) (we also tested the classifier on data from a small study of macrophages²⁰, Supplementary Fig. S17). Because the two studies are substantially different, we took multiple steps to make them comparable. First, we subset to include only those 9,450 genes which were assayed in both. Second, because the dynamic range obtained from RNA-seq (current study) and microarrays (previous study²⁴) were different, we normalized the gene expression levels to a standard normal ($\mu$=0, $\sigma$=1) distribution (Supplementary Fig. S11; note however that this strategy is unable to correct for the inability of microarrays to accurately quantify genes with expression levels that result in fluorescence levels below the background level or above the saturation limit). Third, we corrected for the large, expected batch effect between the two studies by regressing out the first PC of the combined expression data using the limma function removeBatchEffect⁵⁷ (Supplementary Fig. S12).

To identify genes to use in the classifier, we performed a differential expression analysis on the normalized, batch-corrected data from the current study using the same approach described above (with the exception that we no longer used voom⁵⁶ since the data were no longer counts). Specifically, we tested for differential expression between susceptible and resistant individuals in the non-infected state and identified sets of genes to use in the classifier by varying the q-value cutoff. Cutoffs of 5%, 10%, 15%, 20%, and 25% corresponded to gene set sizes of 99, 385, 947, 1,934, and 3,697, respectively. We used the R package caret⁶³ to train 3 different machine learning models: elastic net⁶⁴, support vector machine⁶⁵, and random forest⁶⁶ (the parameters for each individual model were selected using the Kappa statistic). To assess the results of the model on the training data, we performed leave-one-out-cross-validation (LOOCV). In order to choose the model with the best performance, we calculated the difference between the mean of the LOOCV-estimated probabilities of being TB resistant for the samples known to be TB resistant and the corresponding mean for the samples known to be TB susceptible. This metric emphasized the ability to separate the susceptible and resistant individuals into two separate groups. Using this metric, the best performing model was the support vector machine with the 99 genes that are significantly differentially expressed at a q-value of 5% (Fig. 3a, Supplementary Fig. S13, Supplementary Data S5); however, both the elastic net (Supplementary Fig. S14) and random forest (Supplementary Fig. S15) had similar performance. Lastly, we tested the classifier by predicting the probability of being TB susceptible in the 65 healthy samples (Fig. 3b). For evaluating the predictions on the test set of individuals with unknown susceptibility status, we used a relaxed cutoff of the probability of being TB susceptible of 0.25, which was based on the ability of the model at this cutoff to classify all TB susceptible individuals in the training set as susceptible with only 2 false positives. As expected, the 99 genes used in the classifier had similar normalized, batch-corrected median expression levels in the non-infected state across both studies (Supplementary Fig. S16).

**Software implementation.** We automated our analysis using Python (https://www.python.org/) and Snakemake⁶⁷. Our processing pipeline used the general informatics software FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), MultiQC⁶⁸, samtools⁶⁹, and bioawk (https://github.com/lh3/bioawk). We used R⁷⁰ for all statistics and data visualization. We obtained gene annotation information from the Ensembl⁵⁰
and Lynx⁷¹ databases. The computational resources were provided by the University of Chicago Research Computing Center. All code is available for viewing and reuse at https://github.com/jdblischak/tb-suscept.

**Data availability.** The raw fastq files have been deposited in NCBI's Gene Expression Omnibus⁷² and are accessible through GEO Series accession number GSE94116 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE94116). The RNA-seq gene counts and other summary data sets are included as Supplementary Data and are also available for download at https://github.com/jdblischak/tb-suscept.

While a small subset of the GWAS summary statistics required for partially reproducing our results were included in Supplementary Data S4, we do not have permission to share the full set of summary statistics from the previously published TB susceptibility GWAS that would be required for fully reproducing our results. To access the summary statistics, contact the authors directly: Russia - Sergey Nejentsev (sn262@cam.ac.uk), Ghana and The Gambia - Thorsten Thye (thye@bnitm.de), Uganda and Tanzania - Scott M. Williams (smw154@case.edu). The summary statistics for the height GWAS can be downloaded from the GIANT Consortium's website (http://portals.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium_data_files).

## References

1. World Health Organization. Global tuberculosis report 2015. URL http://www.who.int/tb/publications/global_report (2015).

2. World Health Organization. Global TB facts 2015. URL http://www.who.int/tb/Global_TB_Facts.pdf (2015).

3. Glaziou, P., Sismanidis, C., Floyd, K. & Raviglione, M. Global epidemiology of tuberculosis. *Cold Spring Harbor perspectives in medicine* 5, 445–461 (2015).

4. Sotgiu, G., Centis, R., D'ambrosio, L. & Migliori, G. B. Tuberculosis treatment and drug regimens. *Cold Spring Harbor perspectives in medicine* 5, 505–516 (2015).

5. Seung, K. J., Keshavjee, S. & Rich, M. L. Multidrug-resistant tuberculosis and extensively drug-resistant tuberculosis. *Cold Spring Harbor perspectives in medicine* 5, 579–598 (2015).

6. Muñoz, L., Stagg, H. R. & Abubakar, I. Diagnosis and management of latent tuberculosis infection. *Cold Spring Harbor perspectives in medicine* 5, 517–529 (2015).

7. North, R. J. & Jung, Y.-J. Immunity to tuberculosis. *Annual review of immunology* 22, 599–623 (2004).

8. O'Garra, A. et al. The immune response in tuberculosis. *Annual review of immunology* 31, 475–527 (2013).

9. Kallmann, F. J. & Reisner, D. Twin studies on genetic variations in resistance to tuberculosis. *Journal of Heredity* 34, 269–276 (1943).

10. Comstock, G. W. Tuberculosis in twins: a re-analysis of the Prophit survey. *The American review of respiratory disease* 117, 621–624 (1978).

11. Cobat, A. et al. High heritability of antimycobacterial immunity in an area of hyperendemicity for tuberculosis disease. *The Journal of infectious diseases* 201, 15–9 (2010).

12. Möller, M. & Hoal, E. G. Current findings, challenges and novel approaches in human genetic susceptibility to tuberculosis. *Tuberculosis (Edinburgh, Scotland)* 90, 71–83 (2010).

13. Thye, T. et al. Genome-wide association analyses identifies a susceptibility locus for tuberculosis on chromosome 18q11.2. *Nature genetics* 42, 739–41 (2010).

14. Mahasirimongkol, S. et al. Genome-wide association studies of tuberculosis in Asians identify distinct at-risk locus for young tuberculosis. *Journal of human genetics* 57, 363–7 (2012).

15. Thye, T. et al. Common variants at 11p13 are associated with susceptibility to tuberculosis. *Nature genetics* 44, 257–9 (2012).

16. Png, E. et al. A genome wide association study of pulmonary tuberculosis susceptibility in Indonesians. *BMC medical genetics* 13, 1–9 (2012).

17. Chimusa, E. R. et al. Genome-wide association study of ancestry-specific TB risk in the South African Coloured population. *Human molecular genetics* 23, 796–809 (2014).

18. Curtis, J. et al. Susceptibility to tuberculosis is associated with variants in the ASAP1 gene encoding a regulator of dendritic cell migration. *Nature genetics* 47, 523–7 (2015).

19. Sobota, R. S. et al. A locus at 5q33.3 confers resistance to tuberculosis in highly susceptible individuals. *American journal of human genetics* 98, 514–24 (2016).

20. Thuong, N. T. T. et al. Identification of tuberculosis susceptibility genes with human macrophage gene expression profiles. *PLoS pathogens* 4, e1000229 (2008).

21. Bryant, P. A. et al. Susceptibility to acute rheumatic fever based on differential expression of genes involved in cytotoxicity, chemotaxis, and apoptosis. *Infection and immunity* 82, 753–61 (2014).

22. Loddenkemper, R., Lipman, M. & Zumla, A. Clinical aspects of adult tuberculosis. *Cold Spring Harbor perspectives in medicine* 6, a017848 (2016).

23. Cooper, A. M. Cell-mediated immune responses in tuberculosis. *Annual review of immunology* 27, 393–422 (2009).

24. Barreiro, L. B. et al. Deciphering the genetic architecture of variation in the immune response to Mycobacterium tuberculosis infection. *Proceedings of the National Academy of Sciences of the United States of America* 109, 1204–9 (2012).

25. Lango Allen, H. et al. Hundreds of variants clustered in genomic loci and biological pathways affect human height. *Nature* 467, 832–8 (2010).

26. Tailleux, L. et al. Probing host pathogen cross-talk by transcriptional profiling of both Mycobacterium tuberculosis and infected human dendritic cells and macrophages. *PLoS one* 3, e1403 (2008).

27. Deretic, V. Autophagy in tuberculosis. *Cold Spring Harbor Perspectives in Medicine* 4, 53–67 (2014).

28. Castrejón-Jiménez, N. S., Leyva-Paredes, K., Hernández-González, J. C., Luna-Herrera, J. & García-Pérez, B. E. The role of autophagy in bacterial infections. *Bioscience trends* 9, 149–59 (2015).

29. Spang, N. et al. RAB3GAP1 and RAB3GAP2 modulate basal and rapamycin-induced autophagy. *Autophagy* 10, 2297–309 (2014).

30. Sturgill-Koszycki, S. et al. Lack of acidification in mycobacterium phagosomes produced by exclusion of the vesicular proton-atpase. *Science (New York, NY)* 263, 678–81 (1994).

31. Hornef, M. W., Wick, M. J., Rhen, M. & Normark, S. Bacterial strategies for overcoming host innate and adaptive immune responses. *Nature immunology* 3, 1033–40 (2002).

32. Hestvik, A. L. K., Hmama, Z. & Av-Gay, Y. Mycobacterial manipulation of the host cell. *FEMS microbiology reviews* 29, 1041–50 (2005).

33. Flynn, J. L., Goldstein, M. M., Triebold, K. J., Koller, B. & Bloom, B. R. Major histocompatibility complex class I-restricted T cells are required for resistance to Mycobacterium tuberculosis infection. *Proceedings of the National Academy of Sciences of the United States of America* 89, 12013–7 (1992).

34. Grotzke, J. E. et al. The Mycobacterium tuberculosis phagosome is a HLA-I processing competent organelle. *PLoS pathogens* 5, e1000374 (2009).

35. Grotzke, J. E., Siler, A. C., Lewinsohn, D. A. & Lewinsohn, D. M. Secreted immunodominant Mycobacterium tuberculosis antigens are processed by the cytosolic pathway. *Journal of immunology* (Baltimore, Md.: 1950) 185, 4336–43 (2010).
36. Lindestam Arlehamn, C. S., Lewinsohn, D., Sette, A. & Lewinsohn, D. Antigens for CD4 and CD8 T cells in tuberculosis. *Cold Spring Harbor perspectives in medicine* **4**, 89–103 (2014).

37. Miller, M. D. & Krangel, M. S. The human cytokine 1-309 is a monocyte chemoattractant. *Proceedings of the National Academy of Sciences of the United States of America* **89**, 2950–4 (1992).

38. Blischak, J. D., Tailleux, L., Mitrano, A., Barreiro, L. B. & Gilad, Y. Mycobacterial infection induces a specific human innate immune response. *Scientific reports* **5**, 16882 (2015).

39. Sudhof, T. C. The synaptic vesicle cycle. *Annual review of neuroscience* **27**, 509–47 (2004).

40. Berry, M. P. R. et al. An interferon-inducible neutrophil-driven blood transcriptional signature in human tuberculosis. *Nature* **466**, 973–7 (2010).

41. Blankley, S. et al. The application of transcriptional blood signatures to enhance our understanding of the host response to infection: the example of tuberculosis. *Philosophical transactions of the Royal Society of London. Series B, Biological sciences* **369**, 20130427 (2014).

42. Maertzdorf, J., Kaufmann, S. H. & Weiner, J. Toward a unified biosignature for tuberculosis. *Cold Spring Harbor Perspectives in Medicine* **5**, 183–95 (2015).

43. Coscola, M. & Gagneux, S. Does M. tuberculosis genomic diversity explain disease diversity? *Drug discovery today. Disease mechanisms* **7**, e43–e59 (2010).

44. Ostuni, R. et al. Latent enhancers activated by stimulation in differentiated cells. *Cell* **152**, 157–71 (2013).

45. Netea, M. G. et al. Trained immunity: A program of innate immune memory in health and disease. *Science* (New York, NY) **352**, aaf1098 (2016).

46. Lavin, Y. et al. Tissue-resident macrophage enhancer landscapes are shaped by the local microenvironment. *Cell* **159**, 1312–26 (2014).

47. Bloom, C. I. et al. Detectable changes in the blood transcriptome are present after two weeks of antituberculosis therapy. *PloS one* **7**, e46191 (2012).

48. Chaussabel, D. et al. Unique gene expression profiles of human macrophages and dendritic cells to phylogenetically distinct parasites. *Blood* **102**, 672–81 (2003).

49. Liao, Y., Smyth, G. K. & Shi, W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. *Nucleic acids research* **41**, e108 (2013).

50. Yates, A. et al. Ensembl 2016. *Nucleic acids research* **44**, D710–6 (2016).

51. Huber, W. et al. Orchestrating high-throughput genomic analysis with Bioconductor. *Nature methods* **12**, 115–21 (2015).

52. Durinck, S. et al. BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. *Bioinformatics* (Oxford, England) **21**, 3439–40 (2005).

53. Durinck, S., Spellman, P. T., Birney, E. & Huber, W. Mapping identifiers for the integration of genomic datasets with the R/ Bioconductor package biomaRt. *Nature protocols* **4**, 1184–91 (2009).

54. Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. *Bioinformatics* (Oxford, England) **30**, 923–30 (2014).

55. Smyth, G. K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. *Statistical applications in genetics and molecular biology* **3**, Article3 (2004).

56. Law, C. W., Chen, Y., Shi, W. & Smyth, G. K. Voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. *Genome biology* **15**, R29 (2014).

57. Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. *Nucleic acids research* **43**, e47 (2015).

58. Smyth, G. K., Michaud, J. & Scott, H. S. Use of within-array replicate spots for assessing differential expression in microarray experiments. *Bioinformatics* (Oxford, England) **21**, 2067–75 (2005).

59. Liu, R. et al. Why weight? modelling sample and observational level variability improves power in RNA-seq analyses. *Nucleic acids research* **43**, e97 (2015).

60. Stephens, M. False discovery rates: a new deal. *Biostatistics* (Oxford, England) **1–20** (2016).

61. Lawrence, M. et al. Software for computing and annotating genomic ranges. *PLoS computational biology* **9**, e1003118 (2013).

62. Jurasinski, G., Koebsch, F., Guenther, A. & Beetz, S. *Flux: Flux rate calculation from dynamic closed chamber measurements*. URL http://CRAN.R-project.org/package=flux. R package version 0.3-0 (2014).

63. Kuhn, M. Building predictive models in R using the caret package. *Journal of Statistical Software* **28**, 1–26 (2008).

64. Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. *Journal of Statistical Software* **33**, 2008–2010 (2010).

65. Karatzoglou, A., Smola, A., Hornik, K. & Zeileis, A. kernlab - an S4 package for kernel methods in R. *Journal of Statistical Software* **11**, 1–20 (2004).

66. Liaw, A. & Wiener, M. Classification and regression by randomForest. *R News* **2**, 18–22 (2002).

67. Köster, J. & Rahmann, S. Snakemake—a scalable bioinformatics workflow engine. *Bioinformatics* (Oxford, England) **28**, 2520–2 (2012).

68. Ewels, P., Magnusson, M., Lundin, S. & Käller, M. MultiQC: summarize analysis results for multiple tools and samples in a single report. *Bioinformatics* (Oxford, England) **btw354** (2016).

69. Li, H. et al. The sequence alignment/map format and samtools. *Bioinformatics* (Oxford, England) **25**, 2078–9 (2009).

70. R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org (2015).

71. Sulakhe, D. et al. Lynx: a knowledge base and an analytical workbench for integrative medicine. *Nucleic acids research* **44**, D882–7 (2016).

72. Edgar, R., Domrachev, M. & Lash, A. E. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. *Nucleic acids research* **30**, 207–10 (2002).

## Acknowledgements

We thank Matthew Stephens and John Novembre for providing feedback and Gilad lab members for helpful discussion. We thank Marie-Noëlle Ungeheuer and Nathalie Jolly for help recruiting subjects. We thank Sergey Nejentsev for sharing data from the GWAS in Russia. We thank Thorsten Thye for sharing data from the GWAS in Ghana and The Gambia. We thank Rafal S. Sobota, Catherine M. Stein, Giorgio Sirugo, and Scott M. Williams for sharing data from the GWAS in Uganda and Tanzania. This study was funded by National Institutes of Health (NIH) Grant AI087658 to Y.G. and L.T. J.D.B. was supported by NIH Grant T32GM007197. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

## Author Contributions

Y.G., L.T., and L.B.B. conceived of the study and designed the experiments. C.C., E.B., A.D., G.M., O.C., C.V.P., J.H., and R.B. collected samples. L.T. coordinated sample collection and performed the infection experiments. M.M. extracted the RNA and prepared the sequencing libraries. J.D.B. analyzed the results. L.B.B. and Y.G. supervised the project. J.D.B. wrote the paper with input from Y.G., L.T., L.B.B., and R.B.
**Additional Information**

**Supplementary information** accompanies this paper at doi:10.1038/s41598-017-05878-w

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017