Lung function, airway and peripheral basophils and eosinophils are associated with molecular pharmacogenomic endotypes of steroid response in severe asthma.

INTRODUCTION: Asthma is a complex disease with heterogeneous expression/severity. There is growing interest in defining asthma endotypes consistently associated with different responses to therapy, focusing on type 2 inflammation (Th2) as a key pathological mechanism. Current asthma endotypes are defined primarily by clinical/laboratory criteria. Each endotype is likely characterised by distinct molecular mechanisms that identify optimal therapies.
METHODS: We applied unsupervised (without a priori clinical criteria) principal component analysis on sputum airway cells RNA-sequencing transcriptomic data from 19 asthmatics from the Severe Asthma Research Program at baseline and 6-8 weeks follow-up after a 40 mg dose of intramuscular corticosteroids. We investigated principal components PC1, PC3 for association with 55 clinical variables.
RESULTS: PC3 was associated with baseline Th2 clinical features including blood (rank-sum p=0.0082) and airway (rank-sum p=0.0024) eosinophilia, FEV1 change (Kendall tau-b R=-0.333 (-0.592 to -0.012)) and follow-up FEV1 albuterol response (Kendall tau-b R=0.392 (0.079 to 0.634)). PC1 with blood basophlia (rank-sum p=0.0191). The top 5% genes contributing to PC1, PC3 were enriched for distinct immune system/inflammation ontologies suggesting distinct subject-specific clusters of transcriptomic response to corticosteroids. PC3 association with FEV1 change was reproduced in silico in a comparable independent 14-subject (baseline, 8 weeks after daily inhaled corticosteroids (ICS)) airway epithelial cells microRNAome dataset.
CONCLUSIONS: Transcriptomic PCs from this unsupervised methodology define molecular pharmacogenomic endotypes that may yield novel biology underlying different subject-specific responses to corticosteroid therapy in asthma, and optimal personalised asthma care. Top contributing genes to these PCs may suggest new therapeutic targets. © Author(s) (or their employer(s)) 2022. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

Can we discover asthma endotypes using transcriptomic data alone before/after corticosteroid treatment in asthmatic airway cells that are clinically relevant to disease typing and therapeutic response?

In a population of 19 mainly severe asthmatic adults, we show that their sputum airway cell transcriptomic changes 6–8 weeks after a single systemic corticosteroid treatment can be used (without any a priori clinical/laboratory criteria) to define molecular pharmacogenomic endotypes that are associated with type 2 inflammation asthma clinical features and changes in lung function.

This methodology is readily adaptable to similar study designs of corticosteroid treatment response in asthma. The top contributing genes to these endotypes may yield novel biology of corticosteroid treatment response and new therapeutic targets.

Introduction

Asthma is a common complex disease with heterogeneity in both expression and severity. Multiple studies have addressed clinically based endotypes of asthma.1–8 Despite heterogeneity in asthma study populations and their defining clinical features/criteria, some universally recognised common clinical subgroups have emerged. These include mild atopic asthma originating in childhood, atopic asthma with lower lung function, and medication-resistant exacerbation-prone asthma. Clinical endotypes have also been consistently associated with different responses to therapy.9 With the advent of next generation biologics for the treatment of severe asthma, endotyping efforts have focused on identifying type-2 airway inflammation (Th2) asthma.10–14 Currently, asthma endotypes have been defined primarily by clinical laboratory or biomarkers criteria.9 For example, the Th2 clinical features include fractional exhaled nitric oxide (FeNO) >25 ppb, peripheral blood eosinophil count >300 (or >150, 2 cut-offs considered separately) cells/µL, and sputum eosinophils >2%. There has been no reported molecular clustering in asthma without these a priori clinical criteria.15 16 A common theme is that each endotype is likely characterised by distinct molecular mechanisms that may identify optimal therapies for targeting the endotype.1 2 5 6 17 Indeed, several studies have reported distinct transcriptomic differences that appear to correlate with clinical subgroups.7 18 19

The aim of this study is to extend the concept of molecular endotype identification in asthma in novel directions through RNA-sequencing (RNA-seq) of sputum airway cells derived from a subset of 19 asthmatic subjects in the Severe Asthma Research Program (SARP, Wake Forest University centre) at baseline before (baseline), and 6–8 weeks after (follow-up) administration of a 40 mg dose of intramuscular corticosteroid (CS).15 We used the unsupervised principal components (PCs)20–22 of transcriptomic changes following systemic CS treatment to define distinct molecular pharmacogenomic endotypes, derived independently of any a priori clinical criteria and evaluated these endotypes for their association with clinical features (focusing on Th2 ones) and clinical response outcomes.

Methods

Study population and biological samples

The study population consisted of 19 subjects from the National Heart, Lung and Blood Institute Severe Asthma Research Program 3 (SARP 3, Wake Forest University centre Institutional Review Board IRB00021507) ClinicalTrials.gov identifier NCT01750411. Twelve subjects in this longitudinal cohort have severe asthma and seven have non-severe (mild to moderate) asthma. All subjects underwent extensive phenotyping23 including an induced phenotype where each subject had systemic CS responsiveness characterised before (‘baseline’), and 6–8 weeks after (‘follow-up’) the administration of a 40 mg dose of intramuscular CS triamcinolone acetonide.15 Sputum was obtained at both baseline and follow-up using standardised SARP protocols.24–26 Thirty-eight transcriptome profile samples of subject sputum cells total were obtained using RNA-seq: 19 subjects×2 timepoints.

RNA-seq and data analysis

Expanded in online supplemental methods M1. Briefly, unsupervised principal component analysis (PCA) was applied on the subject-wise sample difference (follow-up minus baseline) column vectors of expression of 16 944 robustly detected genes to identify the dominant directions of sample variation in transcriptome space called PCs20 21 focusing on PC1–3. Sample PC coordinates were tested for association with the 55 core clinical variables, online supplemental table S1. Gene ontology (GO) enrichment analysis was performed on the top 5% contributing genes of each PC. All data used in this study will be publicly available on NCBI GEO GSE184433.

Results

The 19 SARP subjects in this study are listed on online supplemental table S2 with select clinical variables at baseline, and their change (follow-up minus baseline) 6–8 weeks after a 40 mg dose of intramuscular CS. Sputum cell samples were obtained at baseline and follow-up, and the 19×2 sputum sample transcriptome profile were obtained using RNA-seq and analysed. Figure 1 outlines the analysis workflow and key results on these transcriptomic data with select clinical variables, and on a comparable independent set of airway epithelial cells (N=14 subjects) microRNAome dataset GSE3446627 on NCBI GEO to assess generalisability of our strategy and results in silico.

Figure 1

Outline of the analysis workflow and key results from the primary study population of 19-subject sputum transcriptome dataset with select clinical variables from the Severe Asthma Research Program, and on a comparable independent 14-subject airway epithelial cells (HBEC) microRNAome dataset GSE34466 on NCBI GEO for assessing the generalisability of the strategy and results in silico. CS, corticosteroid; BASO, basophils; EOS, eosinophils; GCGS, glucocorticoid gene set; GO, Gene ontology; HBEC, human bronchial epithelial cells; NEU, neutrophils; PC#, principal component #; PCA, principal component analysis; SARP, Severe Asthma Research Program.

Sputum cell transcriptomic PC1 and PC3 sample coordinate dichotomies are associated with Th2 clinical features

We applied PCA on the CS treatment sample difference/change (follow-up minus baseline, paired) column vectors of 16 944 genes for all 19 subjects. Without using any clinical variables, PCA identifies the directions of dominant global variance among the 19 subjects in this 16 944-gene dimensional sputum cell transcriptome space. PC1–3 accounted for 28.93%, 19.95% and 11.93% of the global transcriptomic variance, respectively. We plotted the 19 subjects in PC1–3 space, and visual inspection suggests transcriptomic clusters within quadrants of PC1/PC3 space, figure 2A–B.

Figure 2

(A, B) Principal component (PC) analysis of corticosteroid induced transcriptomic change (follow-up minus baseline) sample coordinates in the PC1–3 planes. (C) PC1 sample coordinate versus baseline blood basophil count. (D) PC3 sample coordinate versus baseline blood eosinophil count. Each subject is indicated by their unique 2-digit identifier, 7/19 non-severe asthmatics are in cyan. The PC1 and PC3 dichotomy line at 0 is draw in (A, B). The linear regression line is drawn in (C, D) and the Kendall tau-b correlation and 95% CI are indicated.

We dichotomised PC1 and PC3 coordinates at 0, respectively, to define transcriptomic endotypes, and we checked whether these dichotomies/endotypes were associated with CS treatment response in particular Th2 clinical features (eg, FeNO, sputum eosinophil %) and lung function change, table 1.

Table 1

Select clinical variables (asthma severity, sex, age, body mass index (BMI), pre-bronchodilator (BD) FEV1% predicted) and Th2 clinical feature associations of two molecular endotypes of corticosteroid treatment response defined by corticosteroid induced sputum cell transcriptomic change (=follow up minus baseline) in PC1 and PC3 sample coordinates, respectively, dichotomised at 0 from the 19-subject Severe Asthma Research Program study

Clinical variables	Ratio // median (IQR)		OR (95% CI)//	Ratio // median (IQR)		OR(95% CI)//
	PC1 <0	PC1 ≥0	ranksum p	PC3 <0	PC3 ≥0	ranksum p
# Not severe: severe asthma	2:6	5:6	0.400 (0.055 to 2.933)	4:6	3:6	1.333 (0.204 to 8.708)
# Female: male	6:2	10:1	0.300 (0.022 to 4.060)	8:2	8:1	0.500 (0.037 to 6.684)
Age (years, baseline)	43.43 (41.64–50.40)	45.81 (38.80–52.29)	0.8404	50.40 (45.08–54.47)	41.77 (39.22–43.18)	0.0350*
BMI (baseline)	32.28 (25.76–44.97)	33.41 (32.32–39.98)	0.5448	33.38 (31.32–35.42)	32.66 (30.88–42.61)	0.9048
Baseline blood basophil count (cells/ul)	4 (3–5)	69 (16–95)	0.0191*	5 (4–89)	50 (4–74)	0.9506
Baseline blood eosinophil count (cells/µl)	162 (90–244)	200 (148–394)	0.4073	244 (180–456)	140 (68–188)	0.0082*
Baseline blood neutrophil count (cells/µl)	4025 (2921–5963)	4554 (3968–5752)	0.3950	4295 (3640–5874)	4554 (3725–6057)	0.8421
FeNO (ppb) baseline	17.5 (7.0–37.5)	29.0 (16.3–39.8)	0.3942	32.5 (20.0–38.0)	14.0 (7.0–36.5)	0.3143
Follow-up	14.5 (9.0–27.5)	29.0 (24.3–39.8)	0.0786	27.0 (20.0–32.0)	24.0 (8.5–36.3)	0.6448
Change	−3.0 (−10.0 to 1.5)	−2.0 (−6.3 to 10.5)	0.4785	−4.5 (−7.0 to 2.0)	−1.0 (−6.0 to 4.8)	0.5350
Change, #<0 : #≥0	5 : 2	7 : 4	0.952 (0.144 to 6.281)	7 : 3	5 : 4	1.867 (0.283 to 12.310)
Sputum eosinophil % baseline	0.70 (0.30–1.85)	0.20 (0.00–0.65)	0.1604	1.20 (0.40–2.00)	0.00 (0.00–0.20)	0.0024*
Follow-up	0.50 (0.30–1.10)	0.20 (0.00–0.43)	0.1171	0.40 (0.20–1.00)	0.20 (0.00–0.73)	0.2252
Change	−0.30 (−0.55 to 0.10)	−0.20 (−0.45 to 0.15)	0.6132	−0.45 (−1.50 to 0.20)	0.00 (−0.20 to 0.40)	0.0181*
Change, #<0 : #≥0	6:2	8:3	1.125 (0.141 to 8.995)	9:1	5:4	7.200 (0.622 to 83.342)
Pre-BD FEV1% predicted baseline	70.39 (57.55–86.41)	78.07 (70.71–89.94)	0.3950	78.45 (59.77–86.87)	72.33 (67.42–85.47)	0.9048
Follow-up	70.67 (58.44–84.36)	80.32 (59.27–87.41)	0.5999	82.20 (57.08–87.51)	75.25 (61.80–82.90)	0.7802
Change	−0.17 (−2.05 to 1.89)	−3.43 (−8.38 to 2.03)	0.3101	0.51 (−2.80 to 2.25)	−3.43 (−7.39 to 0.28)	0.2428
Change, #<0 : #≥0	4:4	7:4	0.571 (0.090 to 3.641)	4:6	7:2	0.190 (0.025 to 1.432)

For continuous-valued clinical variables, the median (IQR) and 2-sided Wilcoxon rank-sum p value are shown. For binary-valued clinical variables (eg, asthma severity, change in pre-BD FEV1% predicted dichotomised at 0), the ratios, ORs and Fisher exact test (95% CI) are shown.

*indicates significance at p<0.05.

PC1 dichotomy was significantly associated with baseline blood basophil count (p=0.0191, Wilcoxon rank-sum). PC3 dichotomy was significantly associated with baseline blood eosinophil count (p=0.0082), baseline (p=0.0024) and change (p=0.0181) in sputum eosinophil %. Thus, the molecular endotypes of CS treatment response defined by PC1 and PC3 dichotomies are directly associated with baseline blood basophilia (PC1), baseline blood eosinophilia (PC3), and baseline/change in sputum/airway eosinophilia (PC3). We also noted there were no significant associations with FeNO or pre-bronchodilator (BD) FEV1 percent predicted (%p). PC2 dichotomised at 0 was not associated with any clinical variables here.

Sputum cell transcriptomic PC1 and PC3 sample coordinates associated with baseline blood basophil and eosinophil counts respectively

Next, to increase the power of our transcriptomic endotype analysis, we treated the subject PC1–3 coordinates as continuous variables (not dichotomised as above), and we checked whether these coordinates were associated with 55 core clinical variables, online supplemental table S1). We highlight select clinical variable associations for investigation in table 2 focusing on Th2 clinical features. In addition, we investigated these associations within the 12/19 severe asthma sub-population.

Table 2

Select clinical variable associations of sputum cell transcriptomic change PC1 and PC3 sample coordinates and severity (binary) variable in the Severe Asthma Research Program study using Kendall tau-b correlation where* indicates significance where 0 lies outside the respective 95% CI

Variables tested for association	Kendall tau-b correl R (95% CI)
	PC1 28.93%	PC3 11.93%	Severe	PC1 28.93%(in 12/19 severe)	PC3 11.93%(in 12/19 severe)
PC1 28.93%	–	0.018 (−0.307 to 0.338)	−0.184 (−0.478 to 0.148)	–	0.212 (−0.238 to 0.587)
PC3 11.93%	0.018 (−0.307 to 0.338)	–	0.000 (−0.323 to 0.323)	0.212 (−0.238 to 0.587)	–
Severe	−0.184 (−0.478 to 0.148)	0.000(−0.323 to 0.323)	–	–	–
Sex	−0.022 (−0.342 to 0.303)	0.066 (−0.262 to 0.381)	−0.031 (−0.351 to 0.294)	0.275 (−0.174 to 0.630)	0.110 (−0.334 to 0.514)
Age baseline	0.088 (−0.242 to 0.399)	−0.404 (−0.642 to 0.093)*	−0.367 (−0.617 to 0.051)*	−0.273 (−0.628 to 0.176)	−0.394 (−0.704 to 0.042)
BMI baseline	0.146 (−0.185 to 0.448)	−0.018 (−0.338 to 0.307)	0.033 (−0.292 to 0.352)	0.152 (−0.296 to 0.545)	0.091 (−0.351 to 0.500)
Baseline blood BASO count	0.340 (0.020 to 0.597)*	0.000 (−0.323 to 0.323)	0.025 (−0.300 to 0.345)	0.152 (−0.296 to 0.545)	0.152 (−0.296 to 0.545)
Baseline blood EOS count	0.199 (−0.132 to 0.490)	−0.457 (−0.680 to 0.158)*	−0.050 (−0.367 to 0.277)	0.198 (−0.251 to 0.578)	−0.290 (−0.639 to 0.158)
Baseline blood neutrophil count	0.193 (−0.138 to 0.485)	0.053 (−0.275 to 0.369)	0.367 (0.051 to 0.617)*	0.061 (−0.378 to 0.477)	0.242 (−0.208 to 0.608)
FeNO baseline	0.047 (−0.280 to 0.364)	−0.177 (−0.473 to 0.154)	−0.370 (−0.619 to 0.054)*	0.109 (−0.336 to 0.513)	−0.233 (−0.601 to 0.218)
Follow-up	0.219 (−0.112 to 0.506)	−0.041 (−0.359 to 0.285)	−0.211 (−0.500 to 0.120)	0.264 (−0.186 to 0.622)	−0.047 (−0.466 to 0.390)
Change	0.212 (−0.118 to 0.501)	0.083 (−0.247 to 0.395)	0.438 (0.134 to 0.666)*	0.321 (−0.125 to 0.659)	0.229 (−0.221 to 0.599)
Sputum EOS % baseline	−0.068 (−0.383 to 0.260)	−0.367 (−0.617 to 0.051)*	0.098 (−0.232 to 0.407)	−0.127 (−0.527 to 0.319)	−0.509 (−0.769 to 0.102)*
Follow-up	−0.237 (−0.520 to 0.092)	−0.275 (−0.549 to 0.053)	−0.125 (−0.430 to 0.206)	−0.259 (−0.619 to 0.191)	−0.323 (−0.660 to 0.122)
Change	−0.042 (−0.360 to 0.285)	0.294 (−0.032 to 0.563)	−0.094 (−0.404 to 0.236)	0.016 (−0.416 to 0.441)	0.362 (−0.078 to 0.685)
Sputum macrophage % baseline	−0.029 (−0.349 to 0.296)	0.228 (−0.102 to 0.513)	−0.067 (−0.381 to 0.261)	−0.030 (−0.453 to 0.403)	0.333 (−0.111 to 0.667)
Follow-up	−0.129 (−0.434 to 0.202)	−0.199 (−0.490 to 0.132)	−0.385 (−0.629 to 0.071)*	−0.290 (−0.639 to 0.158)	−0.168 (−0.556 to 0.281)
Change	−0.064 (−0.379 to 0.264)	−0.485 (−0.699 to 0.193)*	−0.050 (−0.367 to 0.277)	−0.061 (−0.477 to 0.378)	−0.545 (−0.789 to 0.153)*
Sputum neutrophil % baseline	0.018 (−0.307 to 0.338)	−0.333 (−0.592 to 0.012)*	−0.150 (−0.451 to 0.181)	−0.121 (−0.523 to 0.324)	−0.303 (−0.647 to 0.144)
Follow-up	0.053 (−0.275 to 0.369)	0.333 (0.012 to 0.592)*	0.217 (−0.114 to 0.504)	0.182 (−0.268 to 0.566)	0.303 (−0.144 to 0.647)
Change	0.181 (−0.150 to 0.476)	0.556 (0.284 to 0.745)*	0.067 (−0.261 to 0.381)	0.273 (−0.176 to 0.628)	0.455 (0.032 to 0.739)*
Pre-BD FEV1%p baseline	−0.064 (−0.379 to 0.264)	0.053 (−0.275 to 0.369)	−0.100 (−0.409 to 0.230)	−0.091 (−0.500 to 0.351)	0.152 (−0.296 to 0.545)
Follow-up	−0.135 (−0.438 to 0.197)	−0.088 (−0.399 to 0.242)	−0.284 (−0.555 to 0.043)	−0.242 (−0.608 to 0.208)	−0.061 (−0.477 to 0.378)
Change	−0.146 (−0.448 to 0.185)	−0.333 (−0.592 to 0.012)*	−0.334 (−0.593 to 0.013)*	−0.273 (−0.628 to 0.176)	−0.455 (−0.739 to 0.032)*
BD response baseline	−0.018 (−0.338 to 0.307)	0.029 (−0.296 to 0.349)	0.100 (−0.230 to 0.409)	0.485 (0.071 to 0.756)*	0.000 (−0.429 to 0.429)
Follow-up	0.135 (−0.197 to 0.438)	0.392 (0.079 to 0.634)*	0.000 (−0.323 to 0.323)	0.576 (0.195 to 0.806)*	0.515 (0.111 to 0.773)*
Change	0.146 (−0.185 to 0.448)	0.146 (−0.185 to 0.448)	−0.067 (−0.381 to 0.261)	0.121 (−0.324 to 0.523)	0.364 (−0.077 to 0.685)

There associations were also studied within the 12/19 severe asthma subpopulation.

*Indicates statistical significance.

BD, bronchodilator; BMI, body mass index.

PC1 was significantly associated with baseline blood basophil count (Kendall tau-b correlation R=0.340 (0.020 to 0.597) 95% CI). PC2 was not associated with any clinical variables here. PC3 with baseline blood eosinophil count (R=−0.457 (−0.680 to –0.158)) and baseline sputum eosinophil % (R=−0.367 (−0.627 to –0.051)), figure 2C–D. PC1 and PC3 were not associated with baseline/follow-up/change in FeNO. Within the 12/19 severe asthma sub-population, the baseline blood basophil (PC1) and eosinophil (PC3) associations were not significant, but the baseline sputum eosinophil % (PC3) became greater (R=−0.509 (−0.769 to –0.102)).

Focusing on the PC1 versus baseline blood basophil count, and PC3 versus baseline blood eosinophil count associations, we used ROC curve analyses to assess the diagnostic performance of PC1 and PC3 as binary classifiers of high versus not high baseline blood basophil and eosinophil counts respectively, as the PC discrimination threshold is varied, figure 3A–B. We consider a subject to have a high blood basophil count if their blood basophil count >10 cells/μL. For blood eosinophil count, we considered two cut-offs: >300 cells/μL yielding just 4/19 (21%) high blood eosinophil subjects, and >150 cells/μL yielding 12/19 (63%) high blood eosinophil subjects. The areas under receiver operator characteristic curve (AUCCH) values of 0.883 (PC1 vs blood basophil), and 0.900, 0.815 (PC3 vs blood eosinophil at >300, >150 cells/μL cut-offs) indicate the CS treatment response PCs to be robust predictors of high baseline blood basophil and eosinophil counts, respectively.

Figure 3

Areas under receiver operator characteristic curve (AUCCH) analyses for: (A) PC1 sample coordinate and baseline blood basophil count >10 cells/μL. (B) PC3 sample coordinate and baseline blood eosinophil (EOS) count >300 cells/μL. (C) PC3 sample coordinate and baseline blood EOS count >150 cells/μL. (D) PC3 sample coordinate and change in pre-BD FEV1% predicted >0.

Sputum cell transcriptomic PC3 sample coordinates sputum macrophage and neutrophil % following CS treatment

While the primary focus in our work here is Th2 clinical feature associations of PC1/3, we noted that PC3 was significantly associated with baseline/follow-up/change in sputum neutrophil % (Kendall tau-b R=0.556 (0.284 to 0.745)), and change in sputum macrophage % (R=−0.485 (−0.699 to –0.193)), table 2, figure 4C–D. Change in sputum neutrophil and macrophage % associations remain significant within the 12/19 severe asthma subpopulation. Blood neutrophil count is not associated with PC1–3.

Figure 4

(A) PC3 sample coordinate versus follow-up bronchodilator response. (B) PC3 sample coordinate versus change in pre-bronchodilator FEV1% predicted. (C) PC3 sample coordinate versus change in sputum macrophage %. (D) PC3 sample coordinate versus change in sputum neutrophil %. Each subject is indicated by their unique 2-digit identifier, 7/19 non severe asthmatics are in cyan. The linear regression line is drawn in each and the Kendall tau-b correlation and 95% CI are indicated.

Sputum cell transcriptomic PC3 sample coordinates associated with change in pre-BD FEV1 following CS treatment

A notable result when the subject PC3 coordinates is treated as a continuous variable (not dichotomised at 0) is that PC3 is significantly associated with two lung function measurements following CS, table 2, figure 4A–B, namely change in pre-BD FEV1%p (Kendall tau-b R=−0.333 (−0.592 to –0.012)), follow-up BD response (R=0.392 (0.079 to 0.634)). Within the 12/19 severe asthma subpopulation both associations became greater, and interestingly PC1 is now significantly associated with baseline (R=0.485 (0.071 to 0.756)) and follow-up (R=0.576 (0.195 to 0.806)) BD response.

Focusing on the change in pre-BD FEV1%p with PC3, we used a ROC curve analysis to assess the diagnostic performance of PC3 for discriminating lung function improvement. We consider a subject to have improved lung function if their change in pre-BD FEV1%p following CS is >0. In this mainly severe asthmatic population (12/19 severe), just 2/19 had >10% increase in pre-BD FEV1%p and 8/19 had >0. The AUCCH value of 0.864 indicate PC3 to be a robust predictor of lung function improvement following CS treatment, figure 3D. That is, PC-defined molecular endotypes resulting from CS treatment response can be used to as a pharmacogenomics predictor of clinical response.

In contrast to a supervised analysis, no baseline/follow-up Th2 clinical features were significantly associated with BD response or with change in pre-BD FEV1%p, table 2, online supplemental table S3. We did note that change in FeNO was associated with change in pre-BD FEV1%p.

Top contributing genes to sputum cell transcriptomic PC1 and PC3 enriched for immune system and inflammation GOs

We regard the top 5% PC loading coefficient magnitude genes as the top contributing genes driving variation along each respective PC coordinate axis. There are ~847/16 944 (~5%) such genes ranked by loading coefficient magnitude per PC after adjusting for ties and repeat genes. We investigated their GO enrichments separately. In contrast to a supervised analysis, we similarly investigated the top 5% genes rank correlated (magnitude) with Th2 clinical features: baseline blood basophil/eosinophil counts and FeNO, and change in FeNO.

First, we found strong overlap between top 5% PC1 genes (OR 4.07 (3.22 to 5.13) Fisher exact 95% CI) and top 5% PC3 genes (OR 3.73 (2.94 to 4.73)) with our previously reported glucocorticoid gene set (GCGS),28 table 3. In contrast, the supervised top 5% Th2 clinical features correlated gene overlaps with GCGS were not significant or had smaller ORs. Also notable is that the top PC1 and PC3 gene sets strongly overlapped (OR 70.28 (59.02 to 83.70)), whereas the supervised top Th2 clinical features correlated gene sets overlaps were less strong (ORs range 0.22–6.29), table 4.

Table 3

Overlaps between previously reported glucocorticoid gene set (GCGS) with the top 5% PC1 and PC3 contributing genes (unsupervised), the top 5% Th2 clinical feature rank correlated genes (supervised) from Severe Asthma Research Program, and target genes of the top 5% PC3 contributing microRNAs (unsupervised) from GSE34466

GCGS overlap with each test variable	PC1 28.93%	PC3 11.93%	Baseline blood baso	Baseline blood eos	Baseline FeNO	Change FeNO	GSE34466PC3	Test variable	GCGS
# genes in overlap	15 598	15 592	15 531	15 542	15 544	15 556	16 040	No	No
	484	490	553	545	545	526	522	No	Yes
	753	759	820	809	807	795	311	Yes	No
	95	89	26	34	34	53	57	Yes	Yes
OR(95% CI)	4.07 (3.22 to 5.13)*	3.73 (2.94 to 4.73)*	0.89(0.60 to 1.33)	1.20(0.84 to 1.71)	1.20(0.84 to 1.71)	1.97 (1.47 to 2.64)*	5.63 (4.19 to 7.57)*

The number of overlapping genes, ORs and Fisher exact test 95% CI are indicated.

*Indicates statistical significance.

Table 4

Overlaps between gene sets from the 19-subject Severe Asthma Research Program study: the top 5% PC1 and PC3 contributing genes (unsupervised), and the top 5% Th2 phenotypes rank correlated genes (supervised)

Top 5% gene overlaps	PC3 11.93%	Baseline blood baso	Baseline blood eos	Baseline FeNO	Change FeNO
PC1 28.93%	510	62	36	23	109	# gene overlap
	70.28*	1.54*	0.84	0.52*	3.06*	OR
	(59.02 to 83.70)*	(1.18 to 2.01)*	(0.60 to 1.18)	(0.34 to 0.79)*	(2.47 to 3.79)*	(95% CI)
PC3 11.93%		47	74	36	82
		1.12	1.90*	0.84	2.14*
		(0.83 to 1.52)	(1.48 to 2.44)*	(0.60 to 1.18)	(1.69 to 2.72)*
Baselineblood baso			73	53	9
			1.88*	1.30	0.20*
			(1.46 to 2.41)*	(0.97 to 1.73)	(0.10 to 0.38)*
Baseline blood eos				84	10
				2.24*	0.22*
				(1.77 to 2.84)*	(0.12 to 0.41)*
ChangeFeNO					180
					6.29*
					(5.24 to 7.54)*

The number of overlapping genes, ORs and Fisher exact test 95% CI are indicated.

*indicates statistical significance.

We separated the top 5% PC1 and PC3 contributing genes by loading coefficient sign prior to GO enrichment analysis, focused on the 10% highest enrichment scoring significant GO clusters, and summarised their keywords, table 5. For PC1 and PC3, we found significant and repeated enrichments of immune system and inflammation related ontologies among the top 10% scoring clusters. Interestingly, steroid metabolism ontologies are found in PC1 positive GO cluster 23, PC3 positive GO cluster 12 and PC3 negative GO cluster 100, data not shown. In contrast, the supervised top 5% Th2 clinical features correlated genes (except change in FeNO) had lesser enrichments for immune system and inflammation related GO clusters, online supplemental table S4.

Table 5

Top 10% significant gene ontology (GO) clusters and their keywords for the top 5% PC1 and PC3 negative and positive contributing genes from the 19-subject Severe Asthma Research Program study

	PC1	PC3
# genes top 5% negative (N)	378	371
# GO clusters (N)	93	70
# genes top 5% positive (P)	470	477
# GO clusters (P)	90	117

Airway epithelial cells microRNAome PC3 sample coordinates associated with change in pre-BD FEV1 following 8-week inhaled CS regimen

To assess the generalisability of our findings in silico, we applied the above unsupervised analysis on an comparable independent microRNA dataset GSE3446627 on NCBI GEO of 862-microRNA profiles of airway epithelial cells from bronchoscopic brushing of 14 steroid-naive (6 weeks before enrolment) asthma subjects at baseline before, and after an 8-week inhaled CS 200 µg two times a day, with pre-BD FEV1 measurements at both timepoints, figure 1. In contrast to the mainly severe asthmatic SARP study population where just 2/19 had >10% and 8/19 had >0% increase in pre-BD FEV1%p, in this population 4/14 had >10% and 8/14 had >0% increase in pre-BD FEV1%p.

First, we applied PCA on the 14-subject 8 weeks minus baseline sample difference column vectors of 862 microRNAs with PC1–3 accounting for 38.64%, 14.71% and 12.13% of the global microRNAome variance respectively, as we had done with the SARP transcriptome data above. Second, we found PC3 to be associated with change (8 weeks minus baseline) in pre-BD FEV1%p (Kendall tau-b R=0.389 (0.001 to 0.675)). We highlight the positive correlation (R=0.389) of pre-BD FEV1%p with the GSE34466 microRNAome PC3, in contrast to the negative correlation (R=−0.333) of pre-BD FEV1%p with the SARP transcriptome PC3 above. Finally, we identified 405 functionally validated target genes (miRTarBase July 15 201929) of the top 5% PC3 contributing microRNAs by their PC3 loading coefficient magnitude and found significant overlaps with our earlier top 5% PC3 contributing genes (OR 3.65 (2.73 to 4.89), data not shown), and with our previously reported GCGS (OR 5.63 (4.19 to 7.57)), table 3.

Discussion

We identified molecular pharmacogenomic endotypes to an induced phenotype of response to systemic CSs in mainly severe asthmatic subjects using unsupervised PCA of a subject’s sputum airway cell transcriptomic changes.15 16 We used resulting the PC1 and PC3 coordinates (figure 2) to define endotypes in two ways: dichotomised at the 0 coordinate point (binary variable) and as is (continuous variable). Respectively, we investigated PC1 and PC3 separately for associations with a set of 55 clinical variables, focusing on Th2 clinical features. We used the top 5% genes contributing to these PCs to uncover biologically and clinically meaningful treatment response genes and ontologies. Crucially, the molecular endotypes defined in this fashion were associated with clinical response to CS. Namely, PC3 is associated with treatment-related change in pre-BD FEV1 and follow-up BD response (table 2). While the molecular endotypes were also related to baseline Th2 status, with PC1 associated with baseline blood basophils, and PC3 with baseline blood and sputum eosinophils (tables 1–2, figures 2–3), none of these baseline Th2 clinical features were themselves associated clinical response to CS (online supplemental table S3). Within the 12/19 severe asthma subpopulation, the baseline blood basophil (PC1) and eosinophil (PC3) associations were diminished, PC3 associations with sputum eosinophil % and change in pre-BD FEV1 and follow-up BD response became greater. The differences between the molecular treatment response endotypes and Th2 status is further emphasised through different biologies represented by ontology enrichment analysis of top loading (contributing) genes from the PC-defined endotypes (table 5) versus the top correlated genes with any Th2 clinical feature (online supplemental table S4). Finally, we demonstrated that the molecular endotypes are generalisable at the level of the airway epithelial cell, in an independent inhaled CS treatment cohort analysis of microRNA data. Overall, molecular treatment endotypes appear to be more representative of CS induced changes in biology and more specific predictors of clinical CS treatment response.

Underscoring the Th2 biologic significance of our findings, the endotype defined by positive PC1 coordinates mainly has subjects with baseline peripheral blood basophilia, figure 3A. Basophils are underappreciated mediators of Th2 immune response that are increasingly recognised as mediators of thymic stromal lymphopoietin induced Th2 inflammation.30 31 In contrast, the endotype defined by negative PC3 coordinates has subjects with the classic Th2 phenotype of blood and sputum eosinophilia,32 33 figure 3B–C. This emphasises the role that Th2 status may play in response to CS34 and supports Th2 status as a key underlying feature of our molecular endotypes. From a clinical response perspective, as the PC3 coordinates increased, there was a significant decrease in lung function change following CS administration, figure 4B. Given that blood/sputum eosinophils were not found to be associated with lung function change, the PC3 molecular endotype contains features that represent a more comprehensive biological picture of treatment response than blood/sputum eosinophils alone. This is consistent with previous work suggesting that molecular features might be superior to clinical characteristics in predicting response to CS in severe asthma from children in SARP.35 The top 5% PC1 and PC3 contributing genes were enriched for CS signalling indicated by their significant enrichment for our previously reported GCGS28 table 3. Thus, these mainly severe asthmatic subjects appear to cluster into at least two distinct types of clinically and biologically relevant molecular responses to CS treatment.

We previously applied unsupervised PCA to identify top 5% PC contributing genes to be associated with human fetal lung development in the context of embryonic gestational age, vitamin D signalling and CS response.20 21 36 While those analyses were aimed at expanding basic understanding of lung development biology, our current study suggests that molecular subgroups resulting from such unsupervised analysis may also reflect distinct biologic mechanisms in response to CS treatment. The method may be generally applied to classifying and predicting CS treatment response. In contrast to supervised analyses using clinical variables, unsupervised PCA (and more generally surrogate variable analysis SVA22) has two distinct advantages. First, a gene’s contribution is not limited to a single PC as shown by the strong overlap between top PC1 and PC3 contributing genes, table 4. Indeed this might reflect the interactome being the basis of a systems biology model, or pleiotropy in general. Second, we can systematically identify and prioritise the top N contributing genes per PC to clarify both individual genes and their ontologies that underlie the endotypic subgroups. We highlight the effectiveness of our method by contrasting it with supervised analysis using the top 5% Th2 clinical features (ie, blood and sputum eosinophilia, high FeNO) correlated genes which ultimately had fewer overlapping genes with each other (table 4) and lesser enrichment for GCGS (table 3), immune system and inflammation related ontologies (online supplemental table S4) suggesting less biologic certainty for prominent involvement of these processes.

Our work complements current clinical cluster analyses. A recent SARP analysis focused on CS response in severe asthmatics18 suggested four potential CS response clusters, with the cluster associated with the greatest CS response demonstrating the highest levels of blood eosinophils. Another current SARP-3 analysis26 identified differences in pre-BD and post-BD lung function longitudinal trajectories between subjects stratified into four groups by baseline sputum eosinophil % high/low+neutrophil% high/low. While our subset analysis does not have sufficient sample size to assess the presence of such clinical clusters, the association of blood eosinophils and sputum neutrophils with PC3 supports a possible overlap of a portion of the clinical endotype with molecular features. Nonetheless, given the deficiencies of supervised analysis to fully capture meaningful ontological associations, it is also likely that molecular endotyping can independently yield additional novel biology when used as an independent metric, as supported by our ontology enrichment analysis. Others and we previously cited the advantage of molecular endotyping for predicting therapeutic response in asthma.37 Fitzpatrick et al 35 noted that clinical features were of limited use for predicting triamcinolone response, whereas systemic mRNA expression of inflammatory cytokines and chemokines related to interleukin (IL)-2, IL-10 and tumour necrosis factor signalling pathways, including, AIMP1, CCR2, IL10RB, and IL5, strongly differentiated children who failed to achieve control with triamcinolone treatment. Our analysis also includes prediction of response using molecular features that were not noted using clinical biomarkers alone,15 16 but in the context of molecular endotypes generated by genome-wide RNA-seq CS response data.

Several studies have described gene expression signatures that may predict CS responsiveness in asthma.35 38 However, these have generally focused on circulating blood cells rather than airway cells, and on baseline differential gene expression differences between predefined responding and non-responding subject populations. Moreover, their prior expression analyses focused on panels of preselected candidate genes. In contrast, our molecular endotyping of each subject is based on whole transcriptomic change in sputum cells following CS treatment. Our model has the unique advantage of capturing dynamic changes inherent in transcriptional response to CS treatment that is agnostic to clinical outcomes/features a priori and gene/molecular preselection biases, and yet our resulting molecular endotypes were found to be associated with clinical variables at baseline and follow-up. Interestingly, 3 genes* from a 6-gene baseline sputum gene expression signature (CLC, CPA3, DNASE1L3, IL1B*, ALPL*, CXCR2*)39 40 reported to discriminate inflammatory phenotypes of asthma and predict ICS treatment response are among the top 5% genes contributing to PC1 and PC3 with IL1B ranked 2, online supplemental table S5, Berthon et al 40 noted the 3 genes to be neutrophil markers and PC3 was significantly associated with baseline/follow-up/change in sputum neutrophils, table 2, figure 4D.

The small sample size of our study is an important multifaceted issue. Clinically, severe asthma is heterogeneous and sample size limits our ability to fully characterise this heterogeneity. With regards to generalising the pre-BD FEV1%p association to other asthma populations, we addressed this technically in silico in a comparable independent microRNA dataset GSE3446627 (cf. Result 6) where we found the resulting microRNAome change PC3 to be associated with change in pre-BD FEV1%p with an opposite signed correlation. This is intriguing considering that expression of microRNAs and their target genes are generally understood to be inversely correlated. Also, MALAT1 (ranked 6 in PC3, and a target of a top 5% GSE34466 microRNAome PC3 hsa-miR-429, online supplemental table S5) is one of five key long non-coding RNAs in an asthma-associated competing endogenous network targeted by eight potential new drugs.41 Despite sample size, we achieved potentially meaningful clinical endotypes. It is notable that neutrophils, eosinophils and basophils (peripheral blood/sputum inflammatory cells clinically known to be associated with asthma severity) are relevant to transcriptomic endotyping for steroid response thus suggesting that a simple Th2 versus non-Th2 dichotomy may be insufficient to describe response to inhaled steroids in severe asthma.

Conclusions

Our analysis of sputum airway cell RNA-seq data in well-characterised mainly severe asthmatics identified distinct molecular pharmacogenomic endotypes in response to an induced phenotype with systemic CSs. Since CSs remain a mainstay of therapy for most asthmatics with persistent symptoms, these endotypes may clarify biological differences underlying the heterogeneity in CS treatment response, and provide an objective useful model for predicting treatment response and developing novel therapies. Our analysis suggest that unsupervised PCA/SVA can be used to objectively define testable (across multiple datasets and populations) molecular endotypes that may yield novel biology underlying different responses to CS therapy in asthma, identify potential new therapeutic targets, and improve the identification of subjects most likely to respond to CSs alone versus those who may benefit from a biologic added to their CS regimen leading to optimal personalised asthma care. Potentially such molecular endotyping could be applied as a diagnostic test using sputum induction and rapid gene expression measurements designed to delineate therapy specific responses in asthma.