A network-based meta-analysis approach identifies a biologically meaningful common gene signature specific to active tuberculosis

First, we aimed to identify a set of common genes consistently differentially expressed in ATB cases relative to other conditions, including HC, LTBI, OLD and Tx, while accounting for cohort/population heterogeneity among the studies. To accomplish this, we developed a network-based meta-analysis approach where we established a gene covariation network to consider genes that were both differentially expressed, and co-varied among multiple studies based on the meta-analysis (Fig 1A). We identified a TB-specific common gene set and used it as the training dataset to build the ML predictive model for TB risk estimation and treatment monitoring. Specifically, we collected whole blood transcriptome data from 27 published studies (the discovery dataset) in which each dataset had ATB and at least one other contrasting condition for a given cohort of individuals (Tables 1 and S1). For each contrasting condition comparison (e.g., ATB vs. HC) we conducted a differential gene expression analysis within a cohort resulting in a list of differentially expressed genes with corresponding log fold-changes (logFC) (Fig 1A1). After repeating the analyses across cohorts, we created an N×M matrix X with logFC values for M cohorts and N genes (Fig 1A2). We further simplified X by converting the matrix values to +1, 0, -1 to account for uncertainty of the data, resulting in a matrix that captured differentially expressed genes that were either upregulated or downregulated between ATB and the other conditions (HC, LTBI, OLD or Tx).

Download:

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

Fig 1. An end–to–end gene signature model development scheme.

(A) The network–based meta–analysis approach. (A1) Schematic of differential gene expression calculations for a single cohort given a disease condition. (A2) Steps for integrating differential gene expression results for multiple cohorts into a gene covariation network. (A3) Four different networks were constructed corresponding to four different disease conditions: ATB v HC, ATB v LTBI, ATB v OLD, ATB v Tx. From these four networks, the top weighted nodes across networks were selected to form the common ATB–specific gene signature. (B) ML predictive model training and testing based on the common gene signature. (C) Steps for the predictive model validation.

https://doi.org/10.1371/journal.pcbi.1010770.g001

Download:

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

Table 1. Cohort lists used for model training and validation.

TB–related datasets and used for training and validation analyses, along with accession IDs, platforms, and which disease comparisons each dataset was used in (see S1 Table for full details including which datasets were used in each network, and details for the respiratory viral infection datasets for model validation). Abbreviations: Net/ML = Network model & ML training/testing; Val = ML validation.

https://doi.org/10.1371/journal.pcbi.1010770.t001

To capture co-varying genes across cohorts, we constructed a gene covariation network from the matrix X where each node represented a gene, and computed the edge weight between nodes by taking the dot product between each pair of genes in X. We only retained edges in the network with weights ≥ 3, which we determined to represent significant covariation from permutation tests (Figs 1A2 and S1). We then calculated the node centrality for each node, or the weighted degree of each node as the sum of the weights of the edges assigned to the node (S2 Fig). We initially observed that the closer to the center (higher weighted degrees) the nodes were placed in their respective networks, the more consistently the genes responded to active TB among various studies (Fig 2H–2K). To define a robust gene signature discriminating ATB from different conditions, we built four covariation networks for each condition comparison and rank ordered genes based on their centrality (Fig 1A3). Finally, we identified a set of 45 genes that were in the top 5% of central genes from the weighted degree distribution of the network, were central in at least two of four networks, and had average logFC ≥ 0.5 or ≤ -0.5 across the studies (Fig 2A–2E).

Download:

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

Fig 2. Network analysis results and gene signature.

(A–D) Network degree distribution fitting. The distribution of weighted degrees for each network was fit with a probability density function to help determine top–ranked genes to retain for the candidate gene set (Methods). (E) Overlap of top 5% of genes by weighted degree in each network. (F) Gene signature mapped on protein–protein interaction network. Protein–protein association network of 45 candidate genes constructed with STRING. Width of edge corresponds to edge confidence. Associations correspond to physical interactions or represent proteins that act functionally together. (G) A visualization of S2 Table, showcasing which pathways (and corresponding genes) were enriched for our 45–gene set (red: gene is in pathway, blue: gene is not part of pathway). (H–K) Volcano plots displaying the mean log2(Fold Change) of gene expression across cohorts against the weighted degree for each network for all genes within each network. Forty–five candidate genes highlighted in orange within each plot. (L–O) Heatmaps of log2(Fold Change) for each of the 45 candidate genes (rows) between disease conditions displayed for each cohort (columns) that were used in constructing each network. The corresponding GSE ID of each cohort was labeled on columns. If a cohort contains multiple clinically defined populations, a specific comparison is also highlighted on the column. Bolded gene names are genes that are also included in the reduced model, genes with an asterisk are not included in either the full or reduced model.

https://doi.org/10.1371/journal.pcbi.1010770.g002

We validated this gene set with a volcano plot showing the average logFC across datasets against the weighted degree of each gene for each network constructed and highlighted the 45 gene set selected (Fig 2H–2K). Consistent with our selection strategy, we found that most genes in this gene set were highly differentiated when averaged across cohorts for each disease comparison and heavily weighted (center) in most of the networks. In addition, we plotted the differential expression of the 45 genes across all datasets (Fig 2L–2O) and observed that gene patterns of upregulation and downregulation were consistent across cohorts within each disease condition and between disease conditions. Genes upregulated in ATB vs. HC for example were also upregulated in ATB vs. OLD, demonstrating a consistent signal for this gene set across cohorts within and between disease conditions. To validate this pattern for the 45-gene set, we averaged the logFC values for each gene across datasets within each comparison (ATB vs HC, ATB vs LTBI, ATB vs OLD, and ATB vs Tx) and computed the correlation between these averaged expression patterns between each pair of comparisons (e.g. ATB vs HC & ATB vs LTBI) to construct a correlation matrix (S3 Fig). We observed that these differential expression patterns correlated strongly across disease conditions with r ranging between 0.61–0.96 and all p-values ≤ 1e-05. However, a few genes (such as HP, AEACAM1) showed slightly inconsistent expression patterns in ATB vs. OLD.

We next sought to identify a biological interpretation for this 45-gene set that contained several genes used in prior TB gene signature studies (S4 Fig). First, we performed a gene set enrichment analysis and found an enrichment of genes involved in immune function such as interferon gamma (IFN-γ) and alpha/beta (α/β) signaling, IL-6 signaling and Toll-like receptor cascades (Fig 2G and S2 Table). Next, we mapped all 45 genes in the protein-protein association network using the STRING database [39]. Strikingly, we found 71% (32/45) of the genes were interconnected (Fig 2F). The connections between genes may indicate physical interaction between some of the encoded proteins, or that the proteins jointly contribute to a shared function [39]. Proteins IFIT2 and IFIT3 centered in the network are IFN-α/β induced proteins and highly expressed during cell-intrinsic immune response to Mtb infection [40]. Further, we observed that STAT1 forms associations with many proteins in the network (n = 17); STAT1 is known to be an immunoregulatory factor that regulates both IFN-α/β and IFN-γ -mediated responses [7,41,42]. Taken together, this report using a network-based meta-analysis approach on whole blood transcriptome datasets provided the first systematic evidence to identify a set of 45 genes that commonly appeared in all 27 studies and was enriched for genes actively involved in interferon activation responding to infection.

Two novel machine-learning (ML) models predict TB disease states

The scope of the ML work was to generate a TB score, a continuous variable, with a fixed range between 0 and 1, representing a dynamic continuum from Mtb infection to disease as well as following disease treatment back to cured. We employed multivariate regression to establish the relationship between outcome (TB disease states) and the expression patterns of the ATB-specific gene signature identified from our network-based meta-analysis (Fig 1B, Materials and Methods). To capture heterogenous responses of individuals in different TB states among the cohorts, we pooled together the whole discovery dataset (27 studies) (Tables 1 and S1), and the ML model regressed to two states: 1 (ATB) and 0 (HC, LTBI, Tx or OLD). The ML framework began with feature selection where the features fed into ML models were derived from an expression ratio between any pair of the common genes (Materials and Methods). Feature selection underwent two steps: univariate and multivariate feature selection. Combining these two steps allowed removal of irrelevant or redundant features from the original 990 gene-paired features down to 41 features which included 42/45 genes (S4 Table) and mitigated the risk of model overfitting while improving predictive performance. We then assessed the best models among 7 ML regression algorithms using the discovery dataset based on 5-repeated 5-fold nested CV (cross-validation) framework. The random forest (RF) model consistently outperformed the other 6 ML models based on the performance metrics (average outer cross-validation R² = 0.51, MSE = 0.11, AUROC = 0.91), while multilayer perceptron model, an artificial neural network approach, was second-best (average outer cross-validation R² = 0.48, MSE = 0.11, AUROC = 0.90) (S3 Table and S5 Fig). We generated a final RF model trained using the full discovery dataset; we refer to this model as the full model, its parameters are listed in S4 Table.

Given that blood transcriptomic signatures hold great promise for development of screening and triage tests but they impose a limit on the number of genes they can probe [30,32,38], we sought to construct a model with fewer features but without reducing predictive performance. We further adapted a stability selection approach to perform a sensitivity analysis and down-selected the most robust features while preserving predictive performance (S6 Fig). Using a forward stepwise selection within a 5-fold nested CV framework, we evaluated the RF algorithm by adding one feature at a time starting from the top ranked features. This final model with 12 features (paired genes) was determined with slightly reduced performance (average outer cross-validation AUROC = 0.85) (S7 Fig) and refer to this model as the reduced model.

As a result, we generated two models, labeled as the full and reduced model (S4 Table). Importantly, among the 18 genes selected by the reduced model, 8/18 genes (ANKRD22, APOL6, BATF2, DUSP3, FCGR1B, GBP4, GBP5, and VAMP5) were repeatedly identified by the best performing gene signature models for incipient TB diagnosis reported previously [43] (S4 Fig), in which DUSP3 and GBP5 were included in the Sweeney 3 gene signature and have been tested in a PoC setting for TB diagnosis [5,38,44]. Seven additional genes (ADM, CD274, FBOX6, LMNB1, IFIT2, NELL2, and ZNF438) have been identified as part of several gene signatures associated with active TB [7,28,45,46]. For the remaining 3 new genes (CD5, LRRK2 and SPOCK2), CD5-expressing regulatory B cells were found to negatively regulate Th17 and to associate with active TB [47]. SPOCK2 expression was shown to be associated with lung injury induced by viral infection [48] and the mouse LRRK2 knockout model suggested that LRRK2 regulates innate immune responses and lung inflammation during Mtb infection [49]. Taken together, we demonstrate that our gene signature discovery approach is robust enough to produce comparable results to other models, and to expand on previous models with better performance.

Validation of the gene signature models in TB disease risk prediction

To systematically validate the models, we collected 10 independent longitudinal studies (7 TB progression and 3 TB treatment cohorts) from multiple countries and different target populations (Fig 1C) (Tables 1 and S1). We first pooled all 7 TB progression datasets [17–19,27,28,31] and stratified the data using different time intervals to disease and LTBI without progression. TB scores generated from the full model showed a graded increase along time intervals beginning >2 years to disease, and the score distributions from all progression groups differed from the no progression group (Fig 3A). Similar observations were found in the reduced model (S8 Fig). Furthermore, one of the progression studies (the Leicester household contact study [31]) categorized subjects into incipient, subclinical, or clinical TB, according to their clinical phenotypes at the time of sampling. The full model was able to statistically separate the subjects’ scores between clinically defined subgroups and healthy controls as well as in-between clinically defined subgroups (with the exception of subclinical TB & clinical TB) (Fig 3B), suggesting that the TB scores derived by the RF model recapitulated host responses to clinical TB pathogenesis. With the same pooled dataset, we compared TB scores generated by four previously published models. These models demonstrated accurate performance in short-term TB risk estimation [3,5,18,27]. We observed similar trends between the scores and time intervals to disease among the models, except for the Sweeney 3 signature; this model may be sensitive to different cohort populations, resulting in the bimodal distributions we observed in several groups (S9 Fig). For most of the published models the score distributions in the groups distant to disease progression (> 12 months) were less discriminant to the no progression group, which poses a challenge on incipient TB diagnosis at early timepoints (S9 Fig).

Download:

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

Fig 3. Systematic model validation using TB progression cohorts.

(A) The distributions of TB scores generated by the full model, stratified by categorical interval to disease (datapoints n = 1281), and (B) clinically defined TB states (datapoints n = 137) are shown in the violin plot. Individual datapoints are plotted as a point in the violin plot. P–values were calculated using the Mann–Whitney U test with Bonferroni correction (ns, p < 1; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001). (C–D) Receiver operating characteristic curves depict diagnostic performance of the model for incipient TB, (C) stratified by time intervals to disease (< 3, < 6, <12, <18, < 24, < 30 months) and (D) mutually exclusive time intervals to disease (0–3, 3–6, 6–12, 12–18, 18–24, 24–30 months). (E) The distributions of TB scores generated by the full model, stratified by different TB disease stages (HC, LTBI, and ATB), other lung disease (OLD) and viral infection (VI), are visualized in the violin. (F) Area under the curve and 95% confidence intervals for each interval to disease are also shown. Diagnostic performance of the model in differentiating between ATB versus HC, LTBI, OLD, and/or viral infection (VI) using a pooled dataset of all 57 collected cohort studies (S1 Table) (datapoints n = 6290) are showed in ROC curves.

https://doi.org/10.1371/journal.pcbi.1010770.g003

We next examined prognostic performance of the models by computing AUROCs, stratified by multiple time intervals to disease (<3, <6, <12, <18, <24, <30 months). Alternatively, we examined mutually exclusive time intervals to disease (0–3, 3–6, 6–12, 12–18, 18–24, 24–30 months) as a sensitivity analysis. For the full model, AUROC for identification of incipient TB within 3 months was able to achieve 0.94 (95% CI 0.91–0.97) and declined along with increasing interval to disease. Despite this, the AUROC for the group over a 2-year period (< 30 months) was still 0.85 (95% CI 0.82–0.88) (Fig 3C and Table 2). Strikingly, we observed similar AUROCs in the reduced model even though this model incorporated fewer features (S8 Fig and Table 2). Both models outperformed the previously published models in discriminating between progressors and non-progressors in all time periods, except for the RISK6 signature which provided a comparable but moderately reduced AUROC (S9 Fig and Table 2). Furthermore, we investigated prognostic performance in each mutually exclusive time-period across the models. AUROC values for the full model ranged from 0.94 (0–3 months, 95% CI 0.91 to 0.97) to 0.74 (18–24 months, 95% CI 0.62 to 0.86); values for the reduced model ranged from 0.94 (0–3 months, 95% CI 0.91 to 0.97) to 0.73 (12–18 months, 95% CI 0.63 to 0.82). We observed no statistical difference between the two models (Table 3). Compared to the previously published models, both of our models had better discrimination between progressors and non-progressors in all mutually exclusive time periods except RISK6, which performed slightly better than our full model 6–12 months prior to disease (Table 3). Remarkably, when assessing prognostic performance, our models outperformed other models, particularly for the groups distant to disease (12–18, 18–24, and 24–30 months) with AUROCs ranging from 0.74 to 0.77 for the full model and 0.73–0.89 for the reduced model (compared to range 0.67–0.70 for RISK6) (Table 3).

Download:

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

Table 2. Prognostic performance of models stratified by inclusive time intervals.

Prognostic performance (AUROC with 95% confidence intervals) of the models developed in this report and published previously on the combined validation datasets for identification of incipient TB within a 2.5–year period, stratified by inclusive time interval to disease.

https://doi.org/10.1371/journal.pcbi.1010770.t002

Download:

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

Table 3. Prognostic performance of models stratified by exclusive time intervals.

Prognostic performance (AUROC with 95% confidence intervals) of the models developed in this report and published previously on the combined validation datasets for identification of incipient TB within a 2.5–year period, stratified by mutually exclusive time interval to disease.

https://doi.org/10.1371/journal.pcbi.1010770.t003

Next, we tested the sensitivity and specificity of each model using cut-offs identified by the maximal Youden Index, which is based on the best tradeoff between sensitivity and specificity from each ROC of the time interval. Both of our models met or approximated the minimum criteria (>75% sensitivity and >75% specificity) in the WHO target product profile (TPP) for prediction of progression to TB disease [50] over a 30-month period (full model sensitivity 74.2% and specificity 78.3%, reduced model sensitivity 77.3% and specificity 74.9%). RISK 6 and Suliman 4 only met the criteria over 0–12 months and 0–6 months, respectively. Sweeney 3 and BATF2 did not meet the criteria in any of the time periods (S5 Table). Based on a pre-test probability of 2% with Youden Index cutoff, we estimated the PPV (positive predictive value) and NPV (negative predictive value) of each model. Both of our models achieved a PPV ranging from 16.8% to 17.7% for a 0–3 month and marginally surpassed the WHO target PPV of > 5.8% over 2 years (full model PPV 6.6% and reduced model PPV 5.9%). None of the other models except RISK 6 met the WHO PPV criteria over 0–2 years (RISK 6 PPV 5.9%) (S5 Table). We then used a pre-specified cutoff of 2 standard deviations (SDs) above the mean of the non-progressing group to prioritize PPV and specificity, which was proposed in previous reports [27,43]. At this threshold, the PPVs estimated in our models significantly improved to 27.9% and 30% over 0–3 months, and 19.5% and 20.3% over 0–24 months for the full and reduced model, respectively, which substantially surpassed the WHO target PPV over 0–2 years. We observed similar improvements in the PPV for the RISK 6 and Suliman 4 models (S6 Table). However, overall sensitivities using this cutoff were significantly reduced compared to the one using Youden Index although specificities remained high (> 96% in all models). Sensitivities over 0–3 months were 66.7% for the full model and 69% for the reduced model. Over 0–30 months, sensitivities were 41.1% for both of our models. Many cases at distant time points to disease progression (>12 months) could not be detected by using this cutoff.

Gene signature models discriminate active TB from non-active TB

Two published gene signature models (Sweeney 3 and RISK 6) demonstrated their utility as a screening or triage test for TB diagnosis in multicenter prospective studies [32,38]. Here we examined diagnostic performance of our models to discriminate active from non-active TB subjects in comparison with the published models using the pooled datasets of all 37 studies. Using this pooled dataset allowed us to recapitulate various confounding variables to mimic a real-world setting. After pooling the datasets, subjects were stratified by different TB states, and the scores for every subject were generated by individual models. Additionally, as Mulenga et al. highlighted, respiratory viral infections induce interferon-stimulated genes, which can be a potential confounding variable for gene signature predictions [37]. We therefore also included an additional 20 datasets containing 15 different respiratory viral infections (S1 Table). Note, since our RF models were generated by using part of these pooled datasets (discovery dataset), the assessment of our model performance in this section was subjected to 5-fold nested CV to ensure a blinded model validation. The full model demonstrated score distributions able to distinguish between the groups and provided a strong discrimination between ATB vs. HC, LTBI or OLD (Fig 3E and 3F). Additionally, the full model was able to discriminate viral infection (VI) from ATB (AUROC 0.93, 95% CI 0.91–0.94), the overall AUROC between ATB and non-ATB (HC, LTBI, OLD and VI) achieved 0.92 (95% CI 0.91–0.93). As expected, the reduced model using fewer features had reduced diagnostic performance but maintained an AUROC at 0.87 (95% CI 0.86–0.89) when separating ATB from non-ATB (S10 Fig). Its reduced performance was primarily attributed to insufficient differentiation between ATB and OLD (AUROC 0.75, 95% CI 0.73–0.77) and between ATB and VI (AUROC 0.83, 95% CI 0.81–0.85) (S10 Fig). In contrast, the pooled TB scores generated by the published models exhibited multimodal or dispersed distributions in multiple groups (S11 Fig). This resulted in less overall discrimination between ATB and non-ATB (Sweeney 3 AUROC 0.81, 95% CI 0.79–0.83; RISK 6 AUROC 0.75, 95% CI 0.73–0.77; BATF2 AUROC 0.70, 95% CI 0.68–0.72; Suliman 4 AUROC 0.79, 95% CI 0.75–0.82) (S11 Fig). The results suggest that published models assessing the expression of fewer genes, especially the BATF2 model testing a single gene, may be more sensitive to confounding variables.

A probabilistic model to assess TB risk

Identifying a positive case from a screening or triage test for disease diagnosis or prognosis usually requires a pre-defined cutoff, which needs to be validated by post-hoc analyses from independent studies. Given heterogeneous blood transcriptomic responses to TB progression due to confounding variables as previously described, this pre-defined cutoff varies in different settings. Similar evidence has been discussed in other reports [3,16,51]. In addition to using the cutoff approach for test readout, here we also propose a new probabilistic model to estimate TB risk using a cutoff-free method while considering individual and population heterogeneity in the model. To build the probability model, we used our full model to compute all TB scores and stratified the datapoints into different TB states (Fig 4A). The score distributions increased according to TB state, beginning with the healthy control and latent infection without progression groups, followed by incipient TB groups ranging from distant to proximal to disease. The active TB group had the highest scores. Using this stratified score data, we generated a probability model to estimate TB risk as a function of TB scores ranging between 0 to 1. The model composed six best-fitted probability functions, each of which estimated the probability associated with the observed outcome (healthy control, latent infection without progression, < 24m, < 12m, < 3m to disease, and active TB), for any given TB score (Fig 4B). For example, a subject with a TB score = 0.2 is estimated to be 82% in LTBI, 58% in HC and only <5% in either incipient or active TB. In contrast, a subject with a TB score = 0.6 is estimated to be 78% in < 24m, 72% in < 12m, 58% in < 3m to disease, 42% in active TB, and only < 5% in either HC or LTBI. As the score increases and approaches 0.8, the probabilities for incipient TB < 3m to disease dramatically increase and approach 100%, as is the case with active TB in which the distribution is shifted further to the right. This probabilistic model allows us to generate an additional readout for TB risk stratification or TB screening while taking into account uncertainty in the data due to biological or technical variation.

Download:

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

Fig 4. Building a probabilistic TB risk model.

(A) The distributions of the full model–generated TB scores, stratified by categorical interval to disease, and active TB, are shown in the violin plot, based on the pooled dataset collected from all the 37 studies. Individual datapoints of the dataset are plotted as a point in the violin plot (datapoints n = 3869). (B) The 6 probability functions (cumulative distribution function) that depict the probabilities of the outcomes (HC, LTBI, < 24m, < 12m, < 3m to disease, and active TB) as a function of TB scores were showed in lines with different colors. The shaded area represented the 95% confidence interval area.

https://doi.org/10.1371/journal.pcbi.1010770.g004

Validation of reduced gene signature model in treatment monitoring and treatment outcome prediction

We have shown that these whole blood gene signatures serve as a correlate of immune status responding to Mtb bacterial load during TB progression. Accumulated evidence has also demonstrated the utility of gene signatures in treatment monitoring along with Mtb elimination [3,16,30–35]. Here we interrogated whether our gene signature model, in particular the reduced model, could be used to infer treatment responses and predict clinical outcomes (cure, failure, and recurrence). We used two independent treatment studies which provided clinical characterization with publicly available transcriptome and bioassay data to assess the consistency of our model prediction with the study results [31,35,52].

The Catalysis treatment response dataset contained 96 HIV-uninfected adults with newly diagnosed pulmonary TB undergoing standard 6-month treatment [52]. Of the 96 participants, 8 did not achieve cure and failed treatment when the sputum culture remained positive at month 6, and 12 of 88 cured participants had recurrent TB within 2 years after treatment completion. We first examined the dynamics of TB scores generated by our model in response to treatment. Since the reduced RF model’s performance in the prognostic assessment is non-inferior to the full model (Tables 2 and 3), and holds greater potential for clinical development, we used only the reduced model for the validation in treatment monitoring. The predicted scores significantly declined during the first week after treatment initiation and continuously decreased with treatment duration. At the end of treatment (EOT) (Day 168), the scores from most of the participants were at a similar level as healthy controls (Fig 5A). The model was able to discriminate between healthy controls and before treatment (Day 0) (AUC 0.92, 95% CI 0.87–0.97), 1 week (AUC 0.76, 95% CI 0.67–0.85), and 4 weeks after treatment initiation (AUC 0.72, 95% CI 0.63–0.81). Poor discrimination between healthy controls and EOT was observed which agreed with the clinical outcome where most of the participants (88 out of 96) were cured at EOT (Fig 5B).

Download:

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

Fig 5. Model validation using the Catalysis TB treatment study.

(A) Comparison of TB scores generated by the reduced model at Day 0, 7, 28, and 168 after treatment initiation and from healthy controls are shown in a violin plot. Individual datapoints are plotted as a point in the violin plot (datapoints n = 395). (B) Receiver operating characteristic curves depict model discrimination between healthy controls and the patients from different time points after treatment initiation. Area under the curve and 95% confidence intervals for each interval to disease are also shown. (C) At each timepoint the TB scores were stratified by the time of sputum culture conversion to negative (negativity at day 28, 56, 84 and 168, and no conversion at day 168 [failed]) and plotted in the violin plot (datapoints n = 360). (D) ROC curves, stratified by different timepoints after treatment initiation, depict predictive performance of the models for discrimination between the patients with bacteriological cure and those with treatment failure at EOT. (E) Comparison of TB scores against treatment outcome from poorly compliant subjects at day 0. P–values were calculated using the Mann–Whitney U test (**, p < 0.01).

https://doi.org/10.1371/journal.pcbi.1010770.g005

We further stratified the participants by time to sputum culture conversion to negative (negativity at day 28, 56, 84 and 168, and no conversion at day 168 [treatment failed]). The participants who failed treatment retained high TB scores throughout the course of treatment, which were statistically different from cured groups. In contrast, at 7 days after treatment initiation the TB scores from the group with the earliest culture conversion to negative (day 28) had approached a similar level as at EOT (Fig 5C and S12). Importantly, the reduced model statistically discriminated between participants with bacteriological cure and those with treatment failure not only at EOT (AUC 0.95, 95% CI 0.83–1.00) but at baseline (AUC 0.73, 95% CI 0.51–0.95) (Fig 5D). RISK 6 and BATF2 also provided a comparable result in treatment outcome prediction (S14A Fig). Strikingly, when we compared TB scores at baseline from poorly adherent subjects against treatment outcome (patients who missed more than 20% of treatment were stratified in the group of poor adherence [52]), the subjects with treatment failure showed statistically higher levels of TB scores than those who were cured (p < 0.01) (Fig 5E).

While treatment failure was statistically associated with poor adherence in this study (p = 0.00032, ANOVA test), we demonstrate that TB scores at baseline serve as another confounding factor that impacts treatment outcome. Moreover, with a portion of cured participants in this study developing recurrent disease within 2 years after treatment completion, we investigated whether the gene signature could differentiate participants with and without recurrence before EOT. The model did not provide enough discriminative evidence to separate the two groups at any timepoint during the course of treatment (S13 Fig) and neither did the published models (S14B Fig). This could be because relapses were not distinguished from reinfections for these recurrent cases, the models should be retested on validated relapses. Moreover, we investigated whether our model contained properties linked to lung inflammation based on [18F] FDG PET-CT imaging responses in this study, which were reported in RISK 6 and Sweeney 3 [3,44]. The scores generated from our model statistically correlated with lung lesion activity as measured by total glycolytic ratio activity (TGRA) at three time points (Day 0 Spearman r = 0.6, p < 0.0001, Day 28 Spearman r = 0.66, p < 0.0001, and Day 168 Spearman r = 0.31, p = 0.0022) (S15C–S15E Fig). Additionally, the TB scores correlated with both MGIT culture time to positivity (Spearman r = -0.62, p < 0.0001) and Xpert Ct values (Spearman r = -0.52, p < 0.0001); both which were measured before treatment (S15A and S15B Fig). Further, the scores at baseline were significantly higher in the participants with radiologically persistent lung inflammation at EOT than those with radiologically cleared lung inflammation (persistent or cleared lung inflammation was defined by the cutoff of TGRA = 400 from a previous report [52]) (S15F Fig). Taken together, these findings suggest that the scores from our model are associated to both metabolic activity of TB lung lesion and active Mtb infection, both of which contribute to EOT outcome prediction.

We further sought to validate our reduced model on data from another independent treatment study; the Leicester treatment cohort is composed of 74 participants with pulmonary TB. Of these participants, 16 received a 6-month standard regimen treatment (clinical cure < 200 days), 39 received extended treatment (requiring >200 days of standard regimen treatment due to clinical suspicion of TB at EOT), 7 were defined as difficult TB cases (requiring extended treatment due to treatment intolerance and/or adherence issues), 4 were TB drug resistance patients, and 8 were infected with a chronic local outbreak TB strain [31]. Consistent with the Catalysis analysis, we observed declining TB scores after treatment initiation. However, the scores did not reach the same level as healthy control until > 1 year after treatment initiation, which might be due to 39/74 participants requiring extended treatment (Fig 6A). The model consistently discriminated between healthy controls and participants in treatment until month 4. The performance of the model began to wane during month 4–6 (AUC 0.70, 95% CI 0.61–0.79), and month 7–12 (AUC 0.71, 95% CI 0.63–0.80), at which point a portion of the patients had been cured based on clinical assessment (Fig 6B).

Download:

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

Fig 6. Model validation using the Leicester TB treatment study.

(A) The distributions of the reduced model–generated TB scores, stratified by different time intervals after treatment initiation and healthy control, are shown in the violin plot. Individual datapoints are plotted as a point in the violin plot (datapoints n = 728). (B) Receiver operating characteristic curves depict model discrimination between healthy controls and the patients from different time intervals after treatment initiation. Area under the curve and 95% confidence intervals for each interval to disease are also shown. Comparison of TB scores between smear positive and negative patients at TB diagnosis (C), between the patients requiring standard and extended anti–TB treatment (ATT) at week 3–4 (D), month 2–3 (E) and month 4–6 (F) after treatment initiation are shown in the violin plots. P–values were calculated using the Mann–Whitney U test (ns, p < 1; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001). (G) Scatterplots show TB scores throughout treatment course, color–coded by the patients requiring standard or extended ATT. The line represents the median of the stratified group with 95% confidence interval around the median shown in the shaded area. (F) ROC curves, stratified by different timepoints after treatment initiation, depict predictive performance of the models for discrimination between the patients requiring standard and extended ATT.

https://doi.org/10.1371/journal.pcbi.1010770.g006

The average TB score was highest at baseline when compared to the different time intervals, however the score distribution was scattered with a long tail (Fig 6A). Strikingly, when we stratified the patients based on smear positive and negative, we found that smear positive patients were statistically more likely to have higher TB scores (p < 0.0001). Most of the samples in the long tail of the score distribution were from smear negative patients (Fig 6C). Furthermore, the scores from patients requiring extended treatment were consistently higher than those on standard treatment at week 3–4 (p > 0.05), month 2–3 (p < 0.01) and month 4–6 (p < 0.0001) (Fig 6D–6F). A previous report identified that smear positive patients mostly fell within the extended treatment patient group [31]. Our model also showed a consistent prediction where patients requiring extended treatment displayed higher scores throughout the course of treatment (Fig 6G). These observations suggest that the duration of treatment required in individual patients could be predicted at early timepoints. Indeed, TB scores at month 2–3 were a stronger predictor of treatment duration and significantly discriminated between the patients requiring standard treatment versus those requiring extended treatment after a 6-month treatment (Fig 6H). In summary, model validation in two independent treatment studies suggests that our model holds potential to monitor TB treatment success and to predict the duration of treatment required for cure at the end of treatment in TB patients.

Transcriptome datasets

We collected 37 publicly available transcriptome datasets from NCBI GEO and EMBL ArrayExpress from previously published studies (Tables 1 and S1). Among the datasets, there were 27 studies containing one or multiple clinical TB states including HC, LTBI, ATB, and treatment (Tx) or other lung diseases (OLD) (datapoints n = 2914), and 10 progression or treatment studies with longitudinal timepoints (datapoints n = 1281). Of these, 26/37 datasets were derived from microarray experiments and 11/37 were derived via RNAseq. We used 27 out of 37 studies (most are the cross-sectional studies) as the discovery datasets to build the networks for our common gene signature discovery, and to train/test the ML model. The remaining 10 longitudinal cohorts as the validation datasets were used to independently validate prediction performance of the ML model on TB prognosis, diagnosis, and responses to TB drug treatment. To further test model specificity to TB disease, we collected additional 20 microarray datasets from GEO (datapoints n = 2095) from patients across 12 countries with multiple respiratory viral infections (S1 Table).

Processing of microarray datasets

Microarray data were downloaded from GEO using the function getGEOData from the MetaIntegrator package v2.1.3 with option qNorm=FALSE [56]. Missing values in each expression matrix (probes x individuals) were replaced with 0 before performing quantile normalization with the package preprocessCore v1.55.2. Rows were reduced for each dataset by retaining only a single probe for each gene (the probe with the highest gene expression sum across all samples in the matrix). This transformed and filtered expression matrix was used for downstream analysis.

Processing of RNAseq datasets

RNAseq raw data were downloaded from NCBI Sequence Read Archive (SRA) and processed through our customized bulk RNAseq pipeline. Fastp v0.20.1 [57] with option qualified_quality_phred = 20 was used to perform quality control, adapter/tail trimming and read filtering on FASTQ data, followed by Salmon v1.4.0 [58] with options recoverOrphans = true, gcBias = true, and seqBias = true for quantification of gene transcript expression from the sequence reads. The trimmed and filtered reads were mapped to the human genome Ensembl GRCh38 (release 102). A read counts matrix was generated from the pipeline for downstream analysis.

Differential gene expression analysis

We performed pairwise differential gene expression analysis on the 27 discovery datasets used in the network analysis. To run differential gene expression between two groups within a given dataset, we required that at least three individuals belong to each group. Four pairwise comparisons were considered for each dataset: (1) ATB vs. HC, (2) ATB vs. LTBI, (3) ATB vs. OLD, and (4) ATB vs. Tx.

Building gene covariation networks

We constructed four gene covariation networks corresponding to four pairwise comparisons as follows: (1) ATB vs. HC, (2) ATB vs. LTBI, (3) ATB vs. OLD, and (4) ATB vs. Tx.

For each comparison, we collected M datasets for each of which we ran a differential gene expression analysis between the groups. From the results, we retained only genes that had both, (1) a corresponding false discovery rate of 5% or less (BH adjusted p-value < 0.05) and (2) an effect size > 1.5-fold (logFC > log₂(1.5)). For each dataset, we then constructed a 1×N vector that presented N genes and contained the logFC for each gene that passed the above filters or 0 otherwise. We combined vectors from M datasets to generate an N×M matrix and convert the matrix to +1 if the corresponding logFC > 0, 0 if the corresponding logFC = 0, or −1 if the corresponding logFC < 0 to create matrix X which held information corresponding to which genes were significantly upregulated and downregulated across all datasets (Fig 1A).

We constructed a covariation network using NetworkX v2.6.3 (https://networkx.org/) with N nodes (each node represents a gene). Edges between nodes were weighted by considering the covariation between each pair of genes in X. We calculated the dot product between each pair of genes (a pair of rows in X) to generate a distribution (D) of edge weights (range—M to M, M = number of datasets). To filter out low-value non-zero edges that may have occurred by chance, we performed a permutation test on X. We constructed a null distribution of edge weights by randomly shuffling the columns of X within each row (S1 Fig). We repeated this process twenty-five times to create the null distribution (D’). For each distribution, D and D’, we computed the proportion of edge weights to the left of each edge weight for t∈{−M, 0} (normalized by all negative edge weights) and the proportion of edge weights to the right of each edge weight for t∈ {0, M} (normalized by all positive edge weights). Let t be the proportion of edge weights computed from D and t’ the proportion of edge weights computed from D’. A threshold T_l for edge weights < 0 and T_r for edge weights > 0 was determined by the exact p-value . The edges within T_l and T_r were then excluded in the network.

Identifying common genes among the networks

First, the weighted degree of a node, the sum of the edge weights of a node, was calculated to represent the degree centrality of a gene in the network (S2 Fig). A gene with higher degree centrality indicates its response to ATB more consistently across the cohorts. To look for the most common genes in the network, we modeled the weighted degree distribution based on the assumption that the covariation network is a scale-free network [62] where the degree distribution follows a power law (Fig 2A–2D). Therefore, a power-law distribution (a generalized beta distribution) is used to estimate the weighted degree distribution W and the cumulative distribution function (CDF) of W was also estimated. For each network, we identified the genes within the top 5% of W. Then we collected the top-ranked genes that occurred in at least two of the networks (Fig 2E). Next, we excluded genes with abs(mean(logFC))<0.5. The final common gene list was composed of 45 genes, where 41 genes were up-regulated and 4 genes were down-regulated in ATB relative to other disease conditions.

Pathway enrichment analysis

To conduct pathway enrichment, we used REACTOME pathway modules taken from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/). We ran our 45-candidate gene set through Enrichr [63] as implemented in Python through the package GSEApy v0.10.1 (https://github.com/zqfang/GSEApy).

Establishing optimized predictive models

Model validation in TB prognosis and responses to treatment

For the progression datasets (7/10 validation cohorts), the datasets were directly compiled without further batch correction to mimic the real-world situation where the prognostic triage tool is performed at an individual subject level. AUC ROC was assessed to evaluate predictive performance that discriminated between the progression group in different prespecified intervals to disease (< 3 months, < 6 months, < 12 months, < 18 months, < 24 months, and < 30 months) and the group including LTBI without progression during the 2-year follow-up. Additionally, sensitivity, specificity, and positive/negative predictive values (PPVs/NPVs) were assessed for each of these time intervals, when assuming 2% pre-test probability, based on two parameters—1) the predetermined cutoffs defined by 2 standard deviations (SDs) above the mean of the control group [27] to prioritize PPVs and specificity (S6 Table), and 2) the maximal Youden Index from AUC ROC which provides the best tradeoff between sensitivity and specificity (S5 Table). Furthermore, to determine the uncertainty of the estimates of model performance metrics, 95% confidence intervals (CIs) were calculated as described in Ying et al [68]. To investigate the changes of the model-predicted TB scores along with disease progression among the cohorts, the data points were binned using mutually exclusive time intervals to disease of 0–3 months, 3–6 months, 6–12 months, 12–18 months, 18–24 months, 24–30 months, and >30 months. For the treatment cohorts (3/10 validation cohorts), AUC ROC was assessed to evaluate the performance of treatment monitoring that distinguished active TB at baseline and at timepoints after treatment initiation. Furthermore, we also examined the prediction of treatment outcomes (cure, failure, and recurrence) using the data at early timepoints (baseline, week 1, week 2, or month 1), assessed by AUC ROC. To study TB score dynamics over the treatment course across cohorts with different sample collection schedules, the data points were bucketed into several representative time points.

Developing a probabilistic TB risk model

To accommodate the heterogeneity of immune responses to TB progression in different demographic groups, we developed a probabilistic model to incorporate the uncertainty of the responses associated with the clinical outcomes. We first compiled all the transcriptomic data (37 cohorts) and extracted the data that was on the trajectory of TB progression as well as from HC, LTBI, and ATB. TB scores for each sample were generated using our predictive model. Next, we focused on six groups–HC, LTBI, < 24 month to disease, < 12 month to disease, < 3 month to disease, and ATB. For each group, we modeled TB scores using a best-fitted probability distribution in the following steps. Initially, the data was to fit to various distributions in which SciPy v1.7.3, a Python library, has a list of predefined continuous distributions that can be tested against the data, and the best-fitted distribution with appropriate parameters estimated was identified by using sum of square error between the actual and fitted distributions. Next, given the best-fitted distribution, its probability density function (PDF) and cumulative density function (CDF) were estimated. Therefore, consider TB score as a continuous random variable X with PDF f_X(x), the CDF can be obtained from:

Since TB score is ranged between 0 and 1, the area under PDF between 0 and 1 must be one.

For the groups of different time intervals to disease and ATB, we calculated the cumulative probability that TB score takes a value less than or equal to t:

Conversely, for the groups of HC and LTBI, we calculated the cumulative probability that TB score takes a value large than or equal to t:

We used the cumulative probabilities to estimate the probability associated to each of six groups, given a TB score from individuals. To determine 95% CIs of CDF, the Dvoretzky-Kiefer-Wolfowitz inequality was used to generate a CDF-based confidence band [69].

Statistics

A two-sided Wilcoxon Rank Sum Test was used for pairwise analyses comparing TB disease scores across any clinically defined groups. The BH method was used to calculate false discovery rate (FDR)- adjusted p-values for multiple-comparison post-hoc correction. An adjusted p-value of less than 0.05 was considered significant. Spearman’s rank analysis was used to test for correlations between variables.

Code

For RNAseq processing pipeline, we used Nextflow (https://www.nextflow.io/) for analytics workflow development and workflow submission into SLURM, a workload manger, to distribute computing jobs on AWS HPC parallel cluster. Microarray/RNAseq data process and differential gene expression analysis were conducted in R with customized scripts. The remaining works, including network analysis, supervised learning model building and validation, and probabilistic modeling were then established in Python. The code for machine learning was adapted to the multi-processing framework to allow parallel computing by leveraging AWS EC2 instance with 64 CPUs.