Transcriptomics of bronchoalveolar lavage cells identifies new molecular endotypes of sarcoidosis

Abstract

Background: Sarcoidosis is a multisystem granulomatous disease of unknown origin with a variable and often unpredictable course and pattern of organ involvement. In this study we sought to identify specific bronchoalveolar lavage (BAL) cell gene expression patterns indicative of distinct disease phenotypic traits.

Methods: RNA sequencing by Ion Torrent Proton was performed on BAL cells obtained from 215 well-characterised patients with pulmonary sarcoidosis enrolled in the multicentre Genomic Research in Alpha-1 Antitrypsin Deficiency and Sarcoidosis (GRADS) study. Weighted gene co-expression network analysis and nonparametric statistics were used to analyse genome-wide BAL transcriptome. Validation of results was performed using a microarray expression dataset of an independent sarcoidosis cohort (Freiburg, Germany; n=50).

Results: Our supervised analysis found associations between distinct transcriptional programmes and major pulmonary phenotypic manifestations of sarcoidosis including T-helper type 1 (Th1) and Th17 pathways associated with hilar lymphadenopathy, transforming growth factor-β1 (TGFB1) and mechanistic target of rapamycin (MTOR) signalling with parenchymal involvement, and interleukin (IL)-7 and IL-2 with airway involvement. Our unsupervised analysis revealed gene modules that uncovered four potential sarcoidosis endotypes including hilar lymphadenopathy with increased acute T-cell immune response; extraocular organ involvement with PI3K activation pathways; chronic and multiorgan disease with increased immune response pathways; and multiorgan involvement, with increased IL-1 and IL-18 immune and inflammatory responses. We validated the occurrence of these endotypes using gene expression, pulmonary function tests and cell differentials from Freiburg.

Conclusion: Taken together, our results identify BAL gene expression programmes that characterise major pulmonary sarcoidosis phenotypes and suggest the presence of distinct disease molecular endotypes.

FIGURE 1

Study workflow. Summary of all steps in the supervised and unsupervised analyses of bronchoalveolar lavage (BAL) from the Genomic Research in Alpha-1 Antitrypsin Deficiency and Sarcoidosis (GRADS) cohort (n=209) together with validation using BAL from the Freiburg cohort (n=50) is presented. WGCNA: weighted gene co-expression network analysis.

FIGURE 2

Association of gene expression with clinical traits and bronchoalveolar lavage (BAL) cell differentials identified using supervised analysis. Scadding staging, pulmonary function tests (PFTs), age, computed tomography (CT) scan features (with severity measurements) and BAL cell differentials were considered continuous variables and Spearman’s ρ test was used. Race, gender and CT scan features (without severity measurements) were considered as categorical and the Wilcoxon rank sum and Kruskal–Wallis tests were used. a) A matrix showing the number of genes associated with each clinical trait on diagonal in yellow and the overlap in assicated genes bewteen any two traits off the diagonal in grey to yellow. Each entry has three rows. For entries on the diagonal, the first row describes the total number of genes assocated with the trait. The second and third rows represent the number of positively and negatively associated genes marked by (+) and (−), respectively. For entries off the diagonal, the first row represents the number of overlapping genes associated with both traits. The second and third rows describe the number of genes positively and negatively associated with both traits marked by (++) and (−−), respectively. Entries coloured grey represent nonsignificant overlap between two traits (p>0.05), while white-to yellow-coloured entries represent a significant (p<0.05) overlap in genes associated with the two traits assessed with the significance assessed by Chi-squared test. b) The top five significant (false discovery rate (FDR) <0.05) enriched pathways for genes associated with each trait as well as overlap between mutliple traits are shown in a Sankey plot. Bars on the left represent genes significantly associated with each clinical trait, with their height representing the number of genes. In the middle of the plot, positively and negatively correlated genes for each clinical trait are represented by yellow and purple coloured bars, respectively, again with heights representing the number of genes. Each set of negatively or positively correlated genes was connected to its top five significant (FDR<0.05) enriched pathways (with at least three genes in the pathway), listed on the right, by grey lines.

FIGURE 3

Idenification of gene modules that associate significantly with clinical traits and bronchoalveolar lavage (BAL) cell differentials using weighted gene co-expression network analysis. a) Heatmap showing the correlation (positive in red and negative in blue) between the Eigen genes of five chosen gene modules (1, 47, 4, 33, 18) and clinical traits (demographics, pulmonary function tests, computed tomography scan features, BAL cell differentials, Scadding stage and treatment within stage II–III and IV). b) Sankey plot visualising the significantly correlated clinical traits as well as the top five significantly enriched pathways (with at least three genes) for each gene module and their overlap. Bars on the left represent clinical traits and are connected to gene modules with which they are significantly correlated (p<0.05 in panel a). Bars in the middle represent the five chosen gene modules with bar heights describing module size. Each gene module was split into genes positively and negatively correlated with its eigen gene represented by yellow (+) and purple (−) bars, respectively. Finally, each gene set was connected to its top five significnatly enriched pathways represented by bars on the right, similarly with bar heights describing the number of module member genes from each pathway. Details of the gene modules can be found at https://yale-p2med.github.io/SARC_BAL. FEV₁: forced expiratory volume in 1 s; FVC: forced vital capacity; D_LCO: diffusing capacity of the lung for carbon monoxide.

FIGURE 4

Identification of novel endotypes of sarcoidosis using weighted gene co-expression network analysis gene modules. Heatmaps show the gene expression pattern of five selected gene modules and their corresponding results of K-means clustering of patients in the Genomic Research in Alpha-1 Antitrypsin Deficiency and Sarcoidosis (GRADS) study as a discovery cohort (left panels) and Freiburg as a validation cohort (right panels). Each row represents a gene and each column represents a patient. Panels a) to e) are heatmaps of the modules 1, 47, 4, 33 and 18, respectively. The patient clusters (endotypes) identified by K-means clustering are labelled as clusters A, B, C and D at the top of the heatmaps.