Baseline characteristics of TBM, PTB, and healthy controls

Association between baseline clinical characteristics with TBM mortality in RNA-seq cohorts

Objectives and cohorts flow

TBM: TB meningitis, HIV: human immunodeficiency virus, PTB: pulmonary TB, HC; healthy controls

Blood transcriptomic profiles of four cohorts: healthy controls (n=30), PTB (n=295), HIV-negative TBM (n=207) and HIV-positive TBM (n=74)

(A) Principle component analysis (PCA) of whole transcriptomic profile of HC, PTB PTB and TBM with and without HIV. Each symbol represents one individual with color coding different cohorts. The x-axis represents principle component (PC) 1, while y-axis represents PC2. (B-G) Enrichment scores of known innate immunity pathways associated with TBM pathogenesis. Gene list of these pathway were depicted in additional file 1, table S5. Pathway enrichment scores were calculated using single sample GSEA algorithm (ssGSEA) (Barbie DA, et al 2009). Each dot represents one participant. The box presents median, 25th to 75th percentile and the whiskers present the minimum to the maximum points in the data.

Blood transcriptomic profiles of three-month mortality at baseline in all TBM and TBM stratified by HIV status

Volcano plot showed differential expression (DE) genes by fold change (FC) between death and survival in all TBM (A), HIV-negative (B) and HIV-positive TBM (C). Each dot represents one gene. The x-axis represents log2 FC. The y-axis showed –log10 FDR of genes. DE genes were colored with red indicating up-regulated, blue indicating down-regulated genes which having fold discovery rate (FDR) <0.05 and absolute FC > 1.5.

Blood transcriptional modules associated with mortality in TBM

(A) Associations between WGCNA modules with two clinical phenotypes TBM disease severity (MRC grade) and three-month mortality in discovery and validation cohorts, and their associated biological processes. The heatmap showed the association between principle component 1 (PC1) of each module and the phenotypes, particularly Spearman correlation r for MRC grade and hazard ratio per increase 1/10 unit of PC1 (HR) for mortality. The HRs were estimated using a Cox regression model adjusted for age, HIV status and dexamethasone treatment. False discovery rate (FDR) corrected based on Benjamini & Yekutieli procedure, with significant level denoted as * < .05, ** < .01 and *** <.001. Gradient colors were used to fill the cell with green indicating negative r or HR < 1, red color indicating positive r or HR > 1. The order of modules was based on hierarchical clustering using Pearson correlation distance for module eigengene. On the left, biological processes, corresponding to modules, were identified using Gene Ontology and KEGG database. (B) Validation of the association between WGCNA modules and mortality in discovery and validation cohorts. X-axis represents –log10 FDR in discovery cohort and Y-axis represents –log10 FDR in validation cohort. Red dash lines indicate FDR = 0.05 as the threshold for statistically significant in both cohorts. Five modules (blue, brown, red, black and cyan) with FDR < 0.05 were validated.

Biological processes, pathways and hub genes of validated modules associated with mortality

(A-D) showed biological processes and pathways identified in four mortality associated modules: blue, brown, red and black module, by over representation analysis (ORA). Bar plots show the top representative GO biological processes or KEGG pathways. The bars indicates biological processes or pathways having ORA FDR < 0.05 and size corresponding to fold enrichment calculated as the ratio of gene number of pathway in the input list divided by the ratio of gene number of the pathway in reference. (E-H) showed gene co-expression networks and hub genes of blue, brown, red and black module, respectively. Each node represents one gene. Each edge represents the link between two genes. Hub genes were shown by bigger nodes and bold text. The gradient color of node corresponds to its HR to mortality, with red indicating HR > 1, and blue HR < 1.

Validation of hub genes in PCR validation cohort

Gene expression of representative hub genes in healthy controls (n=30), PTB (n=295), HIV-negative TBM (n=207) and HIV-positive TBM (n=74)

Each dot represents gene expression from one participant. (A, B) expression of FCAR and MCEMP1 hub genes from the blue module. (C, D) expression of NELL2 and TRABD2A hub genes from the brown modules. (E, F) expression of PLCG1 and NLRC3 hub genes from the red module. (G, H) expression of CD247 and MATK hub genes from the black module. The box presents median, 25th to 75th percentile and the whiskers present the minimum to the maximum points in the data. Comparisons were made between dead (red) with survival (blue) or between HIV-negative and HIV-positive TBM by Wilcoxon rank sum test with p-values displayed as significance level above the boxes and the horizontal bars (* < .05, ** < .01, *** <.001).

Relationship between known pathways associated with TBM pathogenesis and mortality

Enrichment scores of known immune pathways associated with TBM pathogenesis. Gene list of these pathway were depicted in additional file 1, table S5. Pathway enrichment scores were calculated using single sample GSEA algorithm (ssGSEA) (Barbie DA, et al 2009). Each dot represents one participant. The box presents median, 25th to 75th percentile and the whiskers present the minimum to the maximum points in the data. The comparisons were made between survival and death using Wilcoxon rank sum test. Only significant results are presented with * < .05, ** < .01, *** <.001.

Consensus transcriptional modules associated with TBM mortality stratified by HIV-infection

(A) Associations between 16 consensus WGCNA modules with two clinical phenotypes TBM severity (MRC grade) and mortality in HIV-negative (n=207) and HIV-positive (n=74) TBM participants, and their associated BP Gene ontology or KEEG database. The heatmap showed the association between modules and the phenotypes, with Spearman correlation r for MRC grade and hazard ratio per increase 1 unit of PC1 of module (HR) for mortality in HIV-positive and HIV-negative cohorts. The consensus sub-panel presented associations of the consensus modules and clinical phenotypes with same trend detected in both HIV cohorts, otherwise were annotated with missing (NA) values. False discovery rate (FDR) corrected using Benjamini & Yekutieli procedure, with significant level denoted as * < .05, ** < .01 and *** <.001. Gradient colors were used to fill the cell with green indicating negative r or HR < 1, red color indicating positive r or HR > 1. The order of modules was based on hierarchical clustering using Pearson correlation distance for module eigengene. It is noted that these consensus modules were not identical to the identified modules in the primary analysis in Figure 1. (B-C) Functional enrichment analysis of HIV-positive pathway (blue module) and HIV-negative pathway (yellow module), respectively. (D-E) Gene co-expression network of blue and yellow modules. Each node represents one gene. Each edge represents the link between two genes. Hub genes were shown by bigger nodes with bold text. The gradient color of node corresponds to its HR to mortality, with red indicating HR>1, and blue HR<1.

Enrichment score of immunity pathways in healthy controls (n=30), PTB (n=295), HIV-negative TBM (n=207) and HIV-positive TBM (n=74)

Pathway enrichment scores were calculated using single sample GSEA algorithm (ssGSEA) (Barbie DA, et al 2009). Each dot represents one participant. (A-C) showed box-plots depicting enrichment scores of the innate immunity pathways from the blue module. (E-H) enrichment scores of the adaptive immunity pathways from the red and brown modules. (D) normalized expression of TNF. Gene list of these pathway were depicted in additional file 1, table S5. The box presents median, 25th to 75th percentile and the whiskers present the minimum to the maximum points in the data. Comparisons were made between dead (red) with survival (blue) or between HIV-negative and HIV-positive TBM by Wilcoxon rank sum test with p-values displayed as significance level above the boxes and the horizontal bars, respectively (* < .05, ** < .01, *** <.001).

Comparison of gene signatures in distinguishing survival and death in TBM prognostic models

Construction of WGCNA in discovery cohort

Analysis of network topology for various soft-thresholding powers of the top 5,000 most variant genes based on the scale-free network model. (A) showed the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The red horizontal line corresponds to R2 = 0.85 and soft-thresholding power β=8, which was chosen for the construction gene-expression network. (B) displayed the mean connectivity on the y-axis as a function of the soft-thresholding power on the x-axis. The adjacency matrix of the scale-free network between genes was determined as A = (aij), where aij = |cor(genei, genej|β. (C) displays the dendrogram resulted from the hierarchical clustering analysis using topological overlap of the adjacency matrix A served as a dissimilarity metric. Each cluster was referred as a module and assigned with a color. (D) showed bar-plots indicate the number of genes contained in each module.

Preservation of discovery modules in validation cohort

(A) The median rank preservation statistics of the modules. Each module was represented by a point, labeled with the corresponding color and name. The Y-axis represents the median rank of observed preservation statistics per module, while the X-axis indicates the number of genes within each module. A low preservation median rank value indicates a high level of preservation. The gold module was an artificial module comprised of 500 randomly selected genes. The grey module consisted of non-connected genes identified in the WGCNA analysis of the discovery cohort. These two modules exhibited the highest median rank, indicating low preservation, which ensured our preservation analysis controlling well the background noise signal. (B) The Zsummary preservation statistics of the modules. Each module was represented by a point, labeled with the corresponding color and name. The Y-axis represents the Zsummary statistic of each module based on 1,000 permutations of module labels, while the X-axis indicates the number of genes within each module. The horizontal green dash line corresponds to the threshold of Zsummary>10 indicating strongly persevered modules. The presence of all modules in discovery were validated in validation cohort with Zsummary>10.

Correlations between gene module membership and gene significance with mortality in 4 associated modules in discovery and validation cohorts

(A, C, E and G) are scatter plots for blue, brown, red and black modules in discovery cohort. (B, D, F and H) are scatter plots for blue, brown, red and black modules in validation cohort. For each module, each dot represents one gene in the module. The X-axis represents module membership calculated by Pearson correlation between gene expression level and its corresponding PC1 of that module. The Y-axis represents –log10(p-value) of the association between gene expression and three-month mortality based on the Cox regression model adjusted for age, HIV, and dexamethasone treatment. Hub-genes were those on the right corner of the plot. The displayed R2 and P value were taken from linear regression models.

Construction of consensus WGCNA in HIV-negative (n=207) and HIV-positive (n=74)

Consensus analysis of network topology for various soft-thresholding powers of the top 5,000 most variant genes based on the scale-free network model. (A, B) showed the summary network indices (y-axes) as functions of the soft thresholding power (x-axes). Numbers in the plots indicate the corresponding soft thresholding powers. The plots indicated that approximate scale-free topology is attained around the soft-thresholding power of 8 for both sets, which is the lowest power that satisfies the approximate scale-free topology criterion for both cohorts, HIV-positive (black dots) and HIV-negative (red dots). The adjacency matrix of the scale-free network between genes of each cohort was determined and scaled. The consensus topological overlap of two adjacency matrix were identified and served as dissimilarity metric for the hierarchical clustering. (C) displayed the dendrogram of the consensus gene co-expression network. Each cluster was referred as a module and assigned with a color. (D) the bar-plot indicates the number of genes contained in each module.

Gene expression of representative hub genes in healthy controls (n=30), PTB (n=295), HIV-negative TBM (n=207) and HIV-positive TBM (n=74)

Each dot represents gene expression from one participant. (A, B) expression of RNF24 and ZNF721 hub genes from the specific blue module. (C, D) expression of BCL9L and IL11RA hub genes from the specific yellow modules. The box presents median, 25th to 75th percentile and the whiskers present the minimum to the maximum points in the data. Comparisons were made between dead (red) with survival (blue) or between HIV-negative and HIV-positive TBM by Wilcoxon rank sum test with p-values displayed as significance level above the boxes and the horizontal bars, respectively (* < .05, ** < .01, *** <.001).

Pathway fold change enrichment to mortality of top hit pathways in blue, brown and red modules from the primary analysis

Pathway fold change enrichment to mortality of each biological pathway defined as the distribution of the total of mean differential expression between dead and survival TBM of all genes involved in the pathway based on Qusage method (Yaari G at el, 2013). Median and 95% confidence interval of pathway fold change of each top-hit biological pathway in the blue (A), brown (B), and red module (C) in the primary analysis. Blue lines are from HIV-negative cohort and red lines from HIV-positive cohort. Pathway fold change enrichment above 0 indicates up-regulation while below 0 indicates down-regulation in dead. Vertical dashed lines indicate no difference in dead compared to survival TBM.

Gene expression and mortality of selected hub genes in TBM, PTB and healthy participants

(A-E) Boxplots visualized distribution of gene expression of MCEMP1, TRABD2A, CD4, NELL2 and ZNF354C in healthy (n=30), PTB (n=295), TBM without HIV (n=207) and with HIV (n=74). Boxes indicate median and inter-quantile range. Dots indicate data in individuals. The comparisons were made between dead (red) and survival (blue) by Wilcoxon rank sum test with p-value displayed as * < .05, ** < .01 and *** <.001. (F-J) Associations between mortality and gene expression of MCEMP1, TRABD2A, PKD1, NELL2 and ZNF354C. In each figure, the upper panel corresponds to gene distribution in HIV-negative (blue) and HIV-positive (red), the lower panel presents the approximation of the association between that gene and mortality. Their gene expression were divided into 15 groups using equal-distant quantiles (1/15,…, 14/15) intervals and the proportion of mortality within those groups of patients were computed. Each point (error bar) presented the proportion with its confidence interval of mortality per group. The line represents the logistic curve that illustrates the mortality trend corresponding to a two-fold increase in gene expression.

Performance of optimal gene set for TBM mortality prediction

Receiver operating characteristic (ROC) curves for three-month mortality for our developed model using expression level of 4 genes (MCEMP1, NELL2, CD4 and ZNF354C) with 2 clinical predictors (age, MRC grade) (red line) and the model develop by Thao et al. (PMID: 29029055) (blue line) for HIV-negative cohort (A), for HIV-positive cohort (B) and for HIV-negative PCR validation cohort (E). The optimism corrected area under the curve (AUC) was calculated based on 1,000 times bootstrap subsampling as described in the methods. Calibration plots were presented for the corresponding prediction models for HIV-negative (C), HIV-positive cohort (D) and HIV-negative PCR validation cohort (E).

Batch correction for RNA-seq normalized count data

Principle component analysis (PCA) of transcriptomic data before (A) and after batch correction by combat function in SVA r package. Each symbol represents one individual with color coding different cohorts. The x-axis represents principle component (PC) 1, while y-axis represents PC2

Analysis workflow diagram

TBM participants were recruited from two randomized control trials for adjunctive dexamethasone treatment of HIV-negative and HIV-positive adults with TBM (*LAST ACT: NCT03100786; ACT HIV: NCT03092817). Samples from 207 HIV-negative and 74 HIV-positive patients were used for RNA sequencing and randomly divided into discovery and validation cohorts to identify biological pathways and hub genes associated with three-month mortality.