Validation of a 52-gene risk profile for outcome prediction in patients with idiopathic pulmonary fibrosis: an international, multicentre, cohort study.

BACKGROUND: The clinical course of idiopathic pulmonary fibrosis (IPF) is unpredictable. Clinical prediction tools are not accurate enough to predict disease outcomes.
METHODS: We enrolled patients with IPF diagnosis in a six-cohort study at Yale University (New Haven, CT, USA), Imperial College London (London, UK), University of Chicago (Chicago, IL, USA), University of Pittsburgh (Pittsburgh, PA, USA), University of Freiburg (Freiburg im Breisgau, Germany), and Brigham and Women's Hospital-Harvard Medical School (Boston, MA, USA). Peripheral blood mononuclear cells or whole blood were collected at baseline from 425 participants and from 98 patients (23%) during 4-6 years' follow-up. A 52-gene signature was measured by the nCounter analysis system in four cohorts and extracted from microarray data (GeneChip) in the other two. We used the Scoring Algorithm for Molecular Subphenotypes (SAMS) to classify patients into low-risk or high-risk groups based on the 52-gene signature. We studied mortality with a competing risk model and transplant-free survival with a Cox proportional hazards model. We analysed timecourse data and response to antifibrotic drugs with linear mixed effect models.
FINDINGS: The application of SAMS to the 52-gene signature identified two groups of patients with IPF (low-risk and high-risk), with significant differences in mortality or transplant-free survival in each of the six cohorts (hazard ratio [HR] range 2·03-4·37). Pooled data showed similar results for mortality (HR 2·18, 95% CI 1·53-3·09; p<0·0001) or transplant-free survival (2·04, 1·52-2·74; p<0·0001). Adding 52-gene risk profiles to the Gender, Age, and Physiology index significantly improved its mortality predictive accuracy. Temporal changes in SAMS scores were associated with changes in forced vital capacity (FVC) in two cohorts. Untreated patients did not shift their risk profile over time. A simultaneous increase in up score and decrease in down score was predictive of decreased transplant-free survival (3·18, 1·16-8·76; p=0·025) in the Pittsburgh cohort. A simultaneous decrease in up score and increase in down score after initiation of antifibrotic drugs was associated with a significant (p=0·0050) improvement in FVC in the Yale cohort.
INTERPRETATION: The peripheral blood 52-gene expression signature is predictive of outcome in patients with IPF. The potential value of the 52-gene signature in predicting response to therapy should be determined in prospective studies. FUNDING: The Pulmonary Fibrosis Foundation, the Harold Amos Medical Faculty Development Program of the Robert Wood Johnson Foundation, and the National Heart, Lung, and Blood Institute of the US National Institutes of Health.

Introduction

Idiopathic pulmonary fibrosis (IPF) is a progressive and highly lethal interstitial lung disease of unknown etiology. The median survival without transplant is approximately three to four years[1]. The natural history of the disease is highly variable and unpredictable with some patients demonstrating long term clinical stability and others experiencing a more rapid disease course[2]. While clinical parameters allow staging of patients, they do not predict outcome accurately[3]. In recent years, evidence emerged that blood molecular and genetic markers may be indicative of disease outcome and potentially improve the accuracy of clinical predictions[4-9]. However, the majority of these studies were limited in scope and replication.

We previously identified a 52 gene expression signature in peripheral blood mononuclear cells (PBMC) that predicted transplant-free survival (TFS) in IPF and validated four of these genes (CD28, ICOS, LCK, ITK) by qRT-PCR[6]. In this study, we hypothesized that genomic risk profiles based on the peripheral blood, 52-gene expression signature, would accurately predict outcome in IPF. Our objectives were to determine the outcome prediction accuracy of 52-gene risk profiles in multiple cohorts, to determine whether adding 52-gene risk profiles to currently accepted clinical staging tools improved their outcome prediction accuracy, and to identify whether genomic risk profiles change with disease progression or in response to anti-fibrotic therapy.

Methods

The following sections summarize our methods. The online supplement includes more details:

Design, settings and participants

Study design is summarized in figure 1. For time to event analyses, patients were recruited from the Universities of Yale (n=48), Imperial College London (n=55), Chicago (n=45), Pittsburgh (n=120), Freiburg (n=38) and Brigham and Women’s Hospital-Harvard Medical School (BWH-HMS) (n=119) (Table 1). Recruitment started on 07/2004 and ended on 8/2015. For time course analyses, samples were available from Pittsburgh and Yale cohort patients (Figure 1b). IPF diagnosis was established by a multidisciplinary group at each institution following ATS/ERS guidelines[10]. Studies were approved by Institutional Review Boards at each institution and informed consent was obtained from all patients. Demographic, clinical information, pulmonary function test and diffusion capacity of the lung for carbon monoxide (DLCO) were collected at the time of blood draw. The Gender, Age and Lung Physiology (GAP) index was calculated as reported by Ley and colleagues[3].

Figure 1

Study design

The outline summarizes the (a) time to event and (b) time course analysis design for this study including the cohorts, blood compartments, experiments and statistical methods used in each independent cohort and in the pooled data analysis. PBMC: Peripheral blood mononuclear cells. BWH-HMS: Brigham and Women’s Hospital-Harvard Medical School. Dates of enrollment for each cohort are included in figure 1a. For figure 1b, time is presented in years (average and range, in parenthesis).

Table 1

Clinicopathological characteristics of the IPF patients in the cohorts for time to event analysis

FVC%: Forced vital capacity percent predicted, DLCO%: Carbon monoxide diffusing capacity, percent predicted. FEV1% Forced expiratory volume in 1 second, percent predicted. HRCT, high resolution computed tomography. UIP: Usual Interstitial Pneumonia. BWH-HMS: Brigham and Women’s Hospital-Harvard Medical School

Characteristic	Yale University (N=48)	Imperial College London (n=55)	University of Chicago (n=45)	University of Pittsburgh (n=120)	Freiburg University (n=38)	BWH-HMS (n=119)
Age at enrollmentMean ± SD	70·8 ± 6·6	67·3 ± 8·1	67 ± 8·1	68·2 ± 8·5	66·8 ± 8·8	67 ± 8
Gender, n (%)
Males	39 (81·3)	36 (65·5)	40 (88·9)	87 (72·5)	34 (89·5)	82 (68·9)
Females	9 (18·7)	19 (34·5)	5 (11·1)	33 (27·5)	4 (10·5)	37 (31·1)
Race, n (%)
Caucasian	47 (97·9)	52 (94·5)	37 (82·3)	118 (98·34)	38 (100)	108 (90)
Black	0 (0)	0	3 (6·6)	1 (0·83)	0 (0)	6 (5)
Hispanic	1 (2·1)	0	5 (11·1)	0 (0)	0 (0)	2 (1·5)
Other	0 (0)	3 (5·5)	0 (0)	1 (0·83)	0 (0)	3 (2·5)
Smoking status, n (%)
Ever smoker	37 (77·1)	39 (70·9)	27 (60)	80 (66·7)	27 (71·1)	82 (68·9)
Never smokers	11 (22·9)	16 (29·1)	18 (40)	40 (33·3)	11 (28·9)	37 (31·1)
Immunosuppression use, n (%)
No	45 (93·8)	46 (83·6)	43 (95·6)	106 (88·3)	25 (65·8)	90 (75·6)
Yes	3 (6·2)	9 (16·4)	2 (4·4)	14 (11·7)	13 (34·2)	29 (24·4)
Spirometry, mean ± SD
FVC (%)	73·6 ± 15·1	72·8 ± 20·4	61 ± 14·7	66·4 ± 18·6	65 ± 18	65·3 ± 18·5
DLCO (%)	39·6 ± 12·5	39·5 ± 14	43·3 ± 17·7	50·1 ± 18·9	46·8 ± 17·7	42·2 ± 16·2
FEV1 (%)	80·7 ± 19·2	73·5 ± 19	73·9 ± 17·3	78·3 ± 21·2	64·6 ± 16·2	70·8 ± 18·4
GAP Index
Mean ± SD	4·3 ± 1·4	3·9 ± 1·6	4·3 ± 1·6	3·8 ± 1·5	4·4 ± 1·5	3·9 ± 1·3
Diagnosis, n (%)
HRCT + UIP	16 (33·3)	0	24 (53·3)	64 (53·3)	16 (42·1)	66 (55·5)
HRCT Biopsy	32 (66·7)	55 (100)	21 (46·7)	56 (46·7)	22 (57·9)	53 (44·5)

Sample collection, RNA extraction and quality assessment

Yale, Chicago, Pittsburgh and Freiburg cohorts

PBMC collection, total RNA extraction and quality assessment methods have been previously described[6]. BWH-HMS and Imperial College Cohorts: Whole blood was collected using PAXgene blood RNA tubes (PreAnalytiX) and total RNA was extracted using the PAXgene Blood RNA Kit, following the manufacturer’s protocol.

52-gene signature measurement

Yale, Pittsburgh, Freiburg and BWH-HMS cohorts

The nCounter® analysis system (Nanostring)[11] was used to validate the 52-gene signature. Imperial College London cohort: The 52-gene signature was analyzed from a previously published gene expression dataset of whole blood[12] (GEO accession number: GSE93606). University of Chicago cohort: The expression of the 52-gene signature was analyzed from a previously published gene expression dataset of PBMC from IPF patients[6] (GEO accession number: GSE27957). Gene expression microarrays were performed in accordance to MIAME guidelines. Gene normalization was performed by cohort (see online supplement for more details). transformed Log2 gene expression values were used for statistical analyses.

MMP7 measurement

Serum samples were obtained from Pittsburgh cohort patients who had PBMC collected simultaneously (N=114) in the time to event analysis. The MMP7 Elisa assay (R&D Systems) has been previously validated by us [13,14]

Statistical methods and analysis

Development of the Scoring Algorithm of Molecular Subphenotypes (SAMS)

SAMS is a classification algorithm of gene expression data generated from the calculation of two scores (up and down scores). The following steps summarize the calculation of SAMS Up and Down scores: Step 1: Geometric mean normalization - We subtract the log2 value of the gene from the geometric mean of the same gene in all the samples in the cohort. A gene with a positive value is considered increased, and a gene with a negative value is considered decreased. Step 2: Determination of increased and decreased ratios – This ratio is calculated by dividing the number of genes changed in a certain direction (increased or decreased) in a sample divided by the number of genes expected to change in the same direction. The 52-gene signature contains 7 increased and 45 decreased genes. Thus the increased ratio is calculated by dividing the number of actually increased genes by 7 and the decreased ratio is the number of actually decreased genes divided by 45 (see example in online supplement). Step 3: Sums of the values of increased or decreased genes are calculated per sample. Step 4: Calculation of the scores - the up score is derived by multiplying the sum of the values of the increased genes by the increased ratio and the down score by multiplying the sum of the decreased genes by the decreased ratio. Because the gene expression values are log2 – the up score will be positive and the down score will be negative

To determine 52-gene risk profiles in each independent cohort, patients with up scores above the median value and down scores below the median value in each cohort were classified as “high risk”. Patients without this pattern of expression were classified as “low risk”. Analysis of variance (ANOVA) was used to identify significant differences in SAMS scores between cohorts. The SAMS calculator is publicly available at http://gem.med.yale.edu/SAMSWeb3/index.jsp

Time to event analysis

Patients were followed from study entry until death, loss of follow up, or transplant. Because two cohorts (Yale and Imperial college) did not contain transplants, we used different outcome definitions, transplant-free survival (TFS) in Chicago, Pittsburgh, Freiburg and BWH-HMS and mortality in Yale and Imperial College. The association between genomic risk profiles and outcomes was determined by univariate Cox Proportional-Hazard models. For TFS, both transplants and deaths were considered events. To determine whether genomic risk profiles were predictive of the predetermined clinical outcomes, we pooled the data from all cohorts and adjusted for age, gender, percent predicted forced vital capacity (FVC%) and immunosuppressive therapy, defined as the use of prednisone, azathioprine, or a combination of both at the time of blood draw. Multivariate competing risk[15] and Cox proportional-hazard[16] models were applied to the pooled data to determine association with mortality or TFS respectively. For mortality analyses in the pooled data, transplants were considered a competing risk (Figure 1a). Differences in mortality and TFS between patients with high and low risk genomic profiles were evaluated using cumulative incidence and Kaplan Meier curves, respectively. To test whether 52-gene risk profiles could improve outcome prediction when used in combination with the GAP index[3], we fit competing risk[15] and Cox proportional-hazard[16] models as follows: GAP only, genomic only, GAP and genomic or the G-GAP index. The G-GAP index was calculated by adding three points (the maximum score in the GAP index) to the GAP index if a patient had a high risk genomic profile and no points if they had a low risk profile. To determine the prediction accuracy of these models and to compare their predictive performance, we used time-dependent Receiver Operating Characteristic (ROC) for censored data[17] and Area Under the Curve (AUC) using a 10-fold cross validation procedure. Pooled data analysis results were adjusted by patient’s age, gender, immunosuppression use and percent predicted forced vital capacity (FVC%). MMP7 and 52-gene risk profiles where compared head-to-head using the Concordance index (C-index), an equivalent of the area under the curve (AUC) in a receiver operator curve (ROC), a well-accepted measure of the probability that predicting the outcome is better than chance. [18]. GAP index was not included in the comparisons between 52-gene risk profiles and MMP7.

Time course analysis

Time course analyses were performed in the Pittsburgh and Yale, time course cohorts (Figures 1b and S1). Trends in SAMS scores and forced vital capacity volumes (FVC) were plotted to identify shifts in genomic risk profiles over time. To identify statistically significant differences in up and down scores, and FVC across time between high and low risk patients, we used a linear mixed-effect (LME) model[19] with random intercepts. A linear mixed-effect model with random intercepts was also used to study the associations between changes in up and down scores and changes in FVC in patients with simultaneous measurements. LME models were adjusted by patient’s age, gender and therapy (immunosuppression therapy in the Pittsburgh cohort and anti-fibrotic therapy in Yale). Anti-fibrotic therapy was initiated after baseline sample was collected and defined as the use of Pirfenidone or Nintedanib. To determine the association between changes in SAMS scores and survival, we calculated the relative changes in up and down scores for each IPF patient based on their first two visits, adjusting to the duration of time intervals. Thus, for each one of the Pittsburgh cohort patients with at least two visits (N=66), we calculated the relative changes in both scores per month, from one to six months (see supplementary methods). Patients were classified as high risk if relative changes in up and down scores between two subsequent visits occurred simultaneously and were ≥10% (bidirectional changes). A Cox proportional-hazards model was used to determine the association between bidirectional changes in SAMS scores and TFS. Results were adjusted by patient’s age, gender, FVC and immunosuppressive therapy. Finally, we used a LME model to compare the rate of FVC decline per year in IPF patients from the Yale cohort who had a simultaneous decrease in up score and increase in up score (N=6) from those with other time course changes in SAMS scores (N=16), after initiation of anti-fibrotic therapy. Statistical significance was defined as two-sided P<0·05. Analyses were performed using R. Details on the R packages are provided in the online supplement.

Results

52-gene risk profiles are predictive of outcome in IPF and non-inferior to serum MMP7 levels

We measured the expression of a peripheral blood 52-gene signature[6] in IPF patients from six independent academic centers (Table 2). Gene expression levels, clinical and demographic data were collected at baseline in all patients and across time in patients from the Pittsburgh and Yale cohorts to perform time to event (Figure 1a) and time course analysis (Figure 1b), respectively. To classify patients as high or low risk we calculated up and down scores using SAMS. Up or down scores were not significantly different between cohorts suggesting a similar distribution of patients with 52-gene, high risk profiles in each cohort (Figure S2). SAMS scores separated patients into high and low risk groups with impressive similarity in gene expression patterns within risk groups across the various cohorts (Figure 2a). Univariate Cox Proportional-hazard models demonstrated that patients in the high risk group had significantly (P<0·05) higher mortality (Yale and Imperial College London cohort) or lower TFS (Chicago, Pittsburgh, Freiburg and BWH-HMS cohorts), respectively, when compared to patients in the low risk group (Figure 2b). The hazard ratios (HR) for mortality and TFS ranged from 2·03 to 4·37 indicating that patients with a 52-gene, high risk profile had at least a two-fold increased risk of dying or having a lung transplant during follow-up in each independent cohort.

Table 2

Clinicopathological characteristics of the IPF patients in the two risk groups (pooled data) for time to event analysis

P-values were calculated using the Fisher’s exact test except for age, pulmonary function tests and GAP index where an unpaired, two tailed, t-test was used. FVC%, forced vital capacity, percent predicted, DLCO%, carbon monoxide diffusing capacity, percent predicted. FEV1%, forced expiratory volume in 1 second, percent predicted. HRCT, high-resolution computed tomography. UIP, usual interstitial pneumonia.

Characteristics	Low risk (n=278)	High risk (n=147)	P-value†
Age (yr)
Mean ± SD	67·4 ± 7·9	68·4 ± 8·7	0·24
Gender, n (%)
Males	198 (71·2)	120 (81·6)	0·019
Females	80 (28·8)	27 (18·4)
Race, n (%)			0·077
Caucasian	257 (92·4)	143 (97·7)
Black	10 (3·6)	0 (0)
Hispanic	5 (1·8)	3 (2)
Other	6 (2·2)	1 (0·7)
Smoking status, n (%)			0·27
Ever smoker	185 (66·5)	106 (72·1)
Never smoker	93 (33·5)	41 (27·9)
Immunosuppression use, n (%)
No	252 (90·6)	103 (70·1)	<0·0001
Yes	26 (9·4)	44 (29·9)
Spirometry (mean ± SD)
FVC%	69·3 ± 18·4	62·7 ± 17·3	0·0004
DLCO%	46 ± 17·3	40·9 ± 16·2	0·005
FEV1%	76 ± 19·8	70·6 ± 18·4	0·007
GAP Index			0·002
Mean ± SD	3·9 ± 1·4	4·3 ± 1·5
Diagnosis, n (%)			0·41
HRCT+ UIP biopsy	126 (45·3)	60 (40·8)
HRCT	152 (54·7)	87 (59·2)

Figure 2

52-gene risk profiles are predictive of outcome in IPF

(a) Clustering of IPF patients based on52-gene risk profiles (high vs low) derived using SAMS in each one of the six cohorts studied. Every row represents a gene and every column a patient. Color scale is shown adjacent to heat maps in log-based two scale; yellow denotes increase over the geometric mean of samples and purple, decrease. (b) Mortality and Transplant-free survival (TFS) differs between high vs low risk profiles based on the 52-gene signature in each independent cohort.

To determine how outcome prediction using 52-gene risk profiles compared to serum MMP7, we measured MMP7 concentrations by ELISA in Pittsburgh cohort patients with simultaneous PBMC and serum collections (N=114) and compared their TFS prediction performance using the C-index. Our analysis demonstrated that the C-index for TFS prediction in the Pittsburgh cohort was significantly higher (P=0·011) when using 52-gene, genomic risk profiles (C-index=0·72 95%CI 0.659, 0.779) versus MMP7 levels in serum (C-Index=0·61, 95%CI 0.535, 0.683).

52-gene, genomic risk profiles are predictive of outcome independent of demographic and clinical variables

To identify demographic and clinical characteristic differences between 52-gene risk profiles, a pooled data analysis was performed using data from all 425 IPF patients (Figure 3a). High risk patients were predominantly Caucasian males with lower FVC% and DLCO% at presentation. There were more high-risk patients under immunosuppressants (Table 3). A high risk, 52-gene profile was independently predictive of mortality (HR 2·18, 95% CI 1·53, 3·09, P<0·0001) or TFS (HR 2·04, 95%CI 1·52, 2·74, P<0·0001) (Figure 3b and c) after adjusting for age, gender FVC% and immunosuppressive therapy in the pooled dataset. To account for possible cohort heterogeneity, we also performed multivariate competing risk and Cox PH models stratified by cohort in the pooled data and the results did not differ significantly (HR 2·36, 95% CI 1·67, 3·35, P= 1.3e-6 for mortality and HR 2·08, 95% CI 1·54, 2·80, P= 1.6e-6 for TFS). Because of the known adverse effects of immunosuppressive therapy on survival of patients with IPF[20], we repeated the analysis only on patients that did not receive immunosuppression. A 52-gene, high risk genomic profile was also independently predictive of mortality (HR 2·27, 95% CI 1·54, 3·35, P<0·0001) or TFS (HR 2·13, 95%CI 1·54, 2·96 P<0·0001) in this dataset, after excluding patients under immunosuppressants (Figure S3). A prediction model based on the calculated G-GAP index outperformed all other prediction models studied (Supplementary Tables 1 and 2) and significantly improved accuracy prediction of mortality or TFS (Figure 3d and e). The maximal Area Under the Curve changed by 13% (69% to 82%) or 10·6% (70% to 80·6%) for a 30-day mortality and TFS prediction, respectively.

Figure 3

52-gene risk profiles are predictive of outcome independent of demographic and clinical variables

(a) Pooled data analysis comparing high vs low risk profile patients from all cohorts. Color scale is shown adjacent to heat maps in log-based two scale. (b) Mortality and (c) Transplant-free survival (TFS) differs between high vs low risk patients from all cohorts after adjusting for age, gender, FVC% and immunosuppressive therapy. (d) Area Under the Curve (AUC) of time-dependent ROC analysis for (d) mortality and (e)TFS based on the GAP index alone or the G-GAP index in all patients.

Association of 52-gene expression trends over time with disease progression and survival

For time course analyses, we measured the expression of the 52-gene signature in RNA isolated from PBMC using the nCounter system, calculated up and down scores at each time point and collected FVC values over time in two cohorts (Pittsburgh and Yale, Figure 1b). Details about number of visits and follow up duration can be seen in Figure 1.

To determine the association between changes in up and down scores over time with FVC, we performed a LME model adjusted for age and gender in Pittsburgh and Yale cohorts. In both cohorts, up scores were negatively associated with FVC and down scores were positively associated with FVC. The association of up scores with FVC was −0·025 (95%CI −0·039, −0·011, P=0·0004) in the Pittsburgh cohort and −0·010 (95%CI −0·017, −0·004, P=0·004) in the Yale cohort. Similarly, the association of down scores with FVC was 0·008 (95%CI 0·005, 0·011, P<0·0001) in the Pittsburgh cohort, and 0·027 (95%CI 0·004, 0·051, P=0·029) in the Yale cohort.

To determine whether 52-gene, high or low risk patients, not on anti-fibrotic drugs (Pittsburgh cohort) shifted their risk profile, we plotted up and down scores and FVC trends and compared their values across time in high versus low risk groups using a LME model. Our results indicate no shift in risk profiles or FVC trends (Figure 4a, b and c), results confirmed by the LME model. This model demonstrated a significant difference for up scores (high risk: 4·05 vs low risk: 0·99, P<0·0001), down scores (high risk: −14·9 versus low risk: −4·57, P<0·0001) and FVC (high risk: 2·28 liters versus low risk: 2·60 liters, P=0·04) between high and low risk groups across time in this cohort.

Figure 4

52-gene signature trends over time demonstrate association with disease progression and survival

(a) up and down (b) scores from SAMS, and (c) FVC volumes do not shift their trends over time in high (continuous red line) vs low (continuous black line) risk groups (Pittsburgh cohort). Pointwise confidence intervals are represented in purple. (d) Bidirectional changes in SAMS scores (Simultaneous increase in up score and decrease in down score) can be observed during disease course in IPF and are more prominent in high risk individuals (example shown in dotted black line box). (e) Bidirectional changes in SAMS scores are predictive of Transplant-free survival (TFS). Dotted blue line (high risk) –Pittsburgh cohort patients with 30-day bidirectional changes in SAMS scores ≥10%. Continuous red line (low risk) – Pittsburgh cohort patients with 30-day bidirectional changes in SAMS scores <10%. Results adjusted by age, gender, FVC and immunosuppressive therapy.

We also assessed whether substantial changes in SAMS scores over time were predictive of IPF survival in patients not on anti-fibrotic drugs (Pittsburgh cohort). Since relative changes in FVC ≥10% have been associated with decreased IPF survival [21,22], we hypothesized that a relative increase in up score and a simultaneous decline in down score ≥10%, was also predictive of IPF survival. Univariate and multivariate Cox models (Supplementary Table 3) demonstrated that a simultaneous ≥10% increase in up score and decrease in down score (Bidirectional changes), between two measurements obtained 30-days apart (Figure 4d), was significantly predictive of future transplant-free survival (HR: 3·18· 95%CI 1·16, 8·76, P=0·025) (Figure 4e). Only three out of 32 patients in the Yale time course cohort had ≥10% bidirectional changes across time thus we could not assess the relationship between bidirectional changes and survival in this cohort.

Changes in 52-gene expression trends over time are associated with clinical response to anti-fibrotic agents

To determine the effect of anti-fibrotic drugs on 52-gene risk profiles, we first plotted up and down score trends over time in the Yale time course cohort. Low risk profile patients exhibited the same patterns as observed in the Pittsburgh cohort, but high risk profile patients exhibited shifts in up scores (Figure S4a) and down scores (Figure S4b). Because a higher proportion of high risk patients were initiated on anti-fibrotic therapy (90%) compared to low risk patients (59%) (Supplementary table 4), we analyzed the interaction between changes in scores and response to therapy. Impressively, in patients who exhibited a simultaneous decrease in up score and increase in down score, we observed an average increase in FVC (0·06 liters per year), while in patients that did not exhibit these changes in scores, we observed an average decrease in FVC (−0.21 liters per year). The difference was statistically significant (P=0·005) (Figure S4c).

Discussion

We have previously identified a 52-gene signature predictive of TFS in two IPF patient cohorts by using microarray analysis of PBMC[6]. Here, we analyzed the 52-gene signature in the peripheral blood from 425 IPF patients from six independent cohorts. Using the novel Scoring Algorithm of Molecular Subphenotypes (SAMS), we derived risk profiles from the 52-gene signature that identified two classes of IPF patients with significant differences in outcome in all six cohorts. The prediction accuracy of 52-gene risk profiles was better than serum concentrations of MMP7 and adding 52-gene risk profile information to the clinical GAP index significantly increased its prediction accuracy. Temporal analysis revealed that untreated patients generally did not change their risk profiles; however, simultaneous increase in up score and decrease in down scores was predictive of subsequent transplant-free survival. In patients initiated on anti-fibrotic therapy, a simultaneous decrease in up score and increase in down score was associated with stabilization of FVC.

The recognition of the variable clinical course in IPF led to a substantial effort to identify clinical tools and reliable peripheral blood biomarkers for risk stratification. Changes in peripheral blood proteins such as MMP7[4,13], ICAM and IL8[4], SP-A and SP-D[5], KL-6[23], CCL-18[24], YKL40[25], CXCL13[26], POSTN[27], anti-hsp70 IgG antibodies[28] and protease degradation products[29], have been found to be predictive of poor IPF outcomes. Changes in circulating cells (CD4+CD28+ T cells[6], fibrocytes[30] and Semaphorin 7a+ regulatory T cells[31]), gene polymorphisms (TOLLIP[32], TLR3[33] and MUC5B[7]) and aging biomarkers (Telomere length[8], free mitochondrial DNA[34]) have also been associated with mortality in IPF. While these studies strongly suggested the value of peripheral blood biomarkers for risk stratification in IPF, no marker is currently used in clinical practice. This is, in part, because the majority of the studies did not have truly independent replication cohorts, nor did they demonstrate added value over clinical staging tools. In contrast to previous studies, our study provides validation of our 52-gene expression signature in six independent IPF cohorts and demonstrates a substantial improved accuracy when incorporated with currently used clinical tools. This is important, because accurate outcome prediction has very practical implications for IPF patients. Based on the current lung allocation score, and on their clinical characteristics, nearly all of the patients in our study would be referred for transplant evaluation, and many would be eligible for lung transplantation. However, our data suggests that only patients with a high risk genomic profile could require this evaluation urgently, and many may not require lung transplantation even three to five years after diagnosis. Thus, incorporating 52-gene risk profiles in the evaluation of IPF patients, may enhance the precision of lung transplantation referral – avoiding delays in transplants to those who need it early, and delaying those who may not need it. Similarly, when lung transplantation is not an option, this test could also help physicians deciding when to refer IPF patients to palliative care, a currently significant unmet need[35] or distinguish patients who respond to drug therapy from those who do not. Similarly, the majority of previous studies did not assess the change of markers over time. This is important, as it is unknown whether IPF patients shift their risk profiles. We demonstrate that a patient’s 52-gene, genomic risk profile rarely changes in the absence of anti-fibrotic therapy. However, when the profile does change it is important. In untreated patients, a simultaneous increase in up score and decrease in down score reflects subsequent increased mortality.

In patients treated with anti-fibrotic agents, a simultaneous decrease in up score and increase in down score, reflects stabilization or even increase in FVC. Thus, our study demonstrates that 52-gene risk profiles at presentation are predictive of outcome and changes in a patient’s genomic risk profile are informative of clinical deterioration as well as potential response to anti-fibrotic therapies.

While our study focuses on the biomarker applications of the 52-gene signature for risk stratification in IPF, it could also serve to generate hypotheses for follow up studies. We have previously shown that four genes of this signature (CD28, ICOS, LCK and ITK), that belong to the T-cell co-stimulatory signaling pathway, were correlated with the percentage of CD4+CD28+ T cells in the circulation of these patients[6]. Similarly, previous reports have demonstrated that changes in circulating CD4+ T cells with CD28 down-regulation[36] of IPF patients are also associated with poor disease outcomes. These reports suggest a potential link between changes in the expression of genes in the 52-gene signature with phenotypic shifts in circulating immune cells. Similarly, a recent report suggested that down-regulation of T cell co-stimulation markers is associated with T cell exhaustion and poor outcomes in inflammatory and autoimmune diseases[37]. While IPF is not generally considered an autoimmune disease, T cell exhaustion is a mechanism that should be explored as a potential explanation of our findings. Additionally, other members of the 52-gene expression signature may have also some clues about the role of immune aberrations in IPF. As an example, MCEMP1 (mast cell-expressed membrane protein 1) one of the outcome predictive genes when overexpressed, encodes a transmembrane protein isolated from human mast cells[38], known to work in concert with fibroblasts to aggravate pulmonary fibrosis[39] or FLT3 (Fms-related tyrosine kinase 3) a strong Nintedanib-responsive tyrosine kinase with unknown roles in pulmonary fibrosis. While such studies were beyond the scope of this paper, they could potentially shed light on the role of immune aberrations in IPF.

Despite the impressive reproducibility of our findings, we need to recognize some of the limitations of our study. First, SAMS scores were calculated for each individual after normalization within each cohort. The normalization within cohort was required because the data was obtained by different technologies using RNA extracted from whole blood or PBMC (Figure 1). This of course limits the clinical applicability of our results because the expressions of the 52 genes of an entire IPF cohort need to be available for the calculation of the genomic risk profile of an individual patient. For our results to be implemented in the clinic, we would need to generate a set of reference values for the 52 genes in IPF patients. Such reference values could be used to calculate SAMS up and down scores for every new sample and determine the 52-gene risk profile of patients, independently of a specific cohort. The significant reproducible performance of the 52-gene signature, should encourage the development of this reference set and the standardization as an essay for clinical use. Second, we did not determine the specificity of the 52-gene signature to IPF. To assess the effect of aging, we analyzed the 52-gene signature in control individuals older than 90 years of age[40,41] (Figure S5), and found that it was not predictive of mortality in the aged, but we did not study other chronic lung disease. Third, treatment guidelines have changed in IPF in some cohorts, patients were at least initially on immunosuppressive therapy, which it is well known, affect outcome. However, the 52-gene signature was originally discovered in a cohort (Chicago cohort) where only two out of 45 patients were on immunosuppressive therapy at study entry. Such small number of patients under immunosuppression should not account for the transplant-free survival and mortality prediction accuracy of the signature. To further address this, we performed a separate analysis in which we excluded all patients on immunosuppressive therapy at the time of blood draw. The 52-gene signature was predictive of outcome in this population indicating that immunosuppression did not confound our results. Fourth, our initial predictive model was not adjusted to DLCO because we had missing data especially among high risk patients who did not have DLCO measurements performed at the time of blood draw. However, we did address the effect of DLCO indirectly, through the comparison to the GAP index. DLCO is a component of the GAP index, and adding the 52-gene risk profile to the GAP index significantly improved its outcome predictive accuracy. Finally, our longitudinal analysis was limited by the size of cohorts and the difference between them, however we have demonstrated significant reproducibility on two observations, that untreated IPF patients do not generally shift their genomic risk profile and that 52-gene SAMS scores are significantly associated with FVC. The observation that in treated high risk patients, a simultaneous decrease in up score and increase in down score is associated with a significant stabilization of FVC is intriguing, but will require replication, as it is based on a very small number of patients.

In conclusion, our study demonstrates that the 52-gene risk profiles are reproducible predictors of outcome in IPF patients. The enhanced outcome prediction accuracy when 52-gene risk profiles are added to the GAP index (G-GAP index) and the association of changes in genomic risk profiles with changes in FVC, survival and potential response to anti-fibrotic therapy, indicate the potential value of the 52-gene signature as a blood test to risk stratify and monitor disease in IPF. To develop this blood test, we would need prospective studies that specifically address some of the limitations of our study including, the establishment of universal reference values for the 52 genes, a prospective comparison to other molecular markers, and determination whether the 52 gene signature is predictive or associates with acute exacerbations.