Abstract

Introduction:

Idiopathic pulmonary fibrosis (IPF) is a progressive and fatal lung disease with increasing incidence, prevalence, and mortality. It was reported that IPF prevalence in 65 year or older subjects has doubled in the past decade, from 202.2 cases per 100,000 people in 2001 to 494.5 cases per 100,000 people in 2011 (1). The pathological feature of IPF is characterized by temporally heterogeneous interstitial fibrosis. A variety of risk factors have been identified, including older age, male sex and cigarette smoking (2) as well as genetic mutations in a number of genes, such as surfactant proteins (e.g., surfactant protein C (3)), those involved in telomere maintenance (e.g., TERT (4)) and MUC5B (5). Genetic mutations for IPF were comprehensively reviewed by Steele et al. (6).

Transcriptomic analyses of peripheral blood mononuclear cells (PBMCs) or whole blood have been performed to identify important genes/pathways as biomarkers or functional modules for IPF. For example, a study of 120 IPF patients identified CD28, ICOS, LCK and ITK genes and “The costimulatory signal during T cell activation” Biocarta pathway for association with transplant-free survival (7). Additionally, network analyses of the gene expression data from whole blood identified gene modules associated with IPF, highlighting important roles for genes such as NLRC4, PGLYRP1, MMP9, and DEFA4 (8).

These genomics studies, together with others using peripheral blood or PBMCs, have made important contributions to the understanding of genomics and pathogenesis of IPF. However, since PBMCs or whole blood contain heterogeneous populations of different types of immune cells, e.g., T cells, B cells, monocytes, NK cells, etc., and each cell type has their own unique expression signature and function, the overall expression profiles of PBMCs or whole blood represent a “mosaic” of these mixed expression patterns from different cells. Differential expression at the gene level using PBMCs or whole blood may hence be caused by differential “abundance” or “proportions” of one or several cell types. Moreover, a specific type of immune cell may be uniquely important to the IPF disease process and using PBMCs or whole blood that pools different immune cells together may “cloud the signals” of contributions to IPF pathogenesis from a specific cell type.

Motivated by the above reasons we initiated a study of the abundance of immune cells from PBMCs or whole blood in IPF subjects so as to correlate the proportion of a specific immune cell type to the IPF disease outcome. We took advantage of a powerful software, Cibersort (9), which enumerates cell composition of a complex specimen based on deconvolution of the genomic profile of the specimen using signature profiles of different cell types. Using the LM22 signature file (provided by Cibersort (9)) that contains 22 distinct types of immune cells (derived from microarray data of human immune cells), our deconvolution of immune cell types from PBMCs or whole blood has been successful as shown by the extreme low p values (0.000) for each specimen (Appendix 1).

We selected several large (sample size) datasets for genomic analyses of PBMCs or whole blood in IPF with GEO accession numbers of GSE28042 (7), GSE27957 (7), GSE93606 (8) and GSE38958 (10). The total sample size involved is ~300 subjects. After inferring the 22 immune cell proportions from each specimen, we correlated a cell’s proportion to IPF clinical outcomes,4 including the survival time (in two studies (7, 8)), DLco (diffusing capacity for carbon monoxide) (in the third study (10)) or IPF disease status (by comparing with the healthy control subjects) (in the third study (10)).

Our study has identified that the proportions of resting memory (RM) T cells and monocytes are associated with IPF outcomes. Specifically, a higher proportion of RM T cells was associated with a longer survival time in IPF patients, and these cells were at higher proportions in healthy vs. IPF subjects. In contrast, a higher proportion of monocytes was associated with a shorter survival time in IPF patients, and monocytes were at lower proportions in healthy vs. IPF subjects. In addition, a higher proportion of RM T cells was associated with a higher DLco, while a higher proportion of monocytes was associated with a lower DLco. Overall, our studies suggested a protective role for RM T cells and a detrimental role for monocytes in IPF.

Our study provides evidence of the importance of two specific immune cell subsets in the peripheral circulation, i.e., RM T cells and monocytes, for IPF pathogenesis. These data may form the basis for IPF diagnosis/prognosis based on counting specific immune cells from blood. It also provides justification for more focused analyses (e.g., RNA-seq or epigenomics analyses) of these two specific cells in terms of their mechanistic functions in IPF pathogenesis.

Results:

Datasets 1 and 2 combined:

We then combined the data (the re-scaled cell proportion data and patient characteristics data) from these 45 IPF subjects (Dataset 2) with those from the 75 IPF subjects in Dataset GSE28042 (Dataset 1), forming a combined dataset with 120 IPF patients.

These 120 subjects have a mean age of 68.3 years (43–85 years, sd (standard deviation) = 8.2 years) and a gender make-up of 28 females vs. 92 males. Among them, 58 subjects had observed survival-related outcome (death or lung transplantation) and 62 were censored at the end of the study. The YTO has a mean of 1.79 years (from 0.008 to 4.19 years, sd = 1.10 years).

The summary statistics of cell proportions (in percent) for these 120 subjects are shown in Appendix2, Table 4.

As shown in the table, 7 cell types (e.g., Dendritic.cells.resting) have a low variation (standard deviation < 1 percent) across all the subjects. This low variation is due to 0 or near 0 percent of proportion for the cell type for a large number of subjects. We therefore removed the 7 cell types and kept the remaining 14 cell types for downstream analyses.

Using the combined data, we performed survival analysis (via survival regression model, assuming exponential distribution of survival time) to correlate the YTO with a cell’s proportion, adjusting for age and gender. Adjusting for multiple testing using FDR, several cell types have significant association with IPF survival at FDR<0.05, which include monocytes, activated NK cells, RM T cells, M2 macrophages, and naïve T cells.

Table 1 shows survival regression coefficient for each cell type that achieved an FDR <0.05 in the survival model as well as the p value for the coefficient and the FDR for multiple testing correction.

Table 1:

Cells whose proportions are significantly associated with survival of IPF patients

Cell type	Regression Coefficient	P value	FDR
Monocytes	-0.04	3.51E-04	7.72E-03
NK.cells.activated	-0.37	7.42E-04	8.03E-03
T.cells.CD4.memory.resting	0.08	1.10E-03	8.03E-03
Macrophages.M2	0.51	3.99E-03	2.20E-02
T.cells.CD4.naive	0.09	8.35E-03	3.68E-02

Open in a separate window

Note: Data come from dataset 1 (GSE28042) and dataset 2 (GSE27957).

As shown in the table, proportions of monocytes and activated NK cells appear to have negative associations with IPF survival, i.e., the higher the proportions, the shorter the survival time. On the other hand, RM T cells, M2 macrophages and naïve T cells appear to have positive association with IPF survival, i.e., the higher the proportions, the longer the survival time.

According to the table, on average, a 10 percentage increase of monocyte proportion may accelerate the time to event by a factor of [exp(−0.04)]^10, which is 0.67, i.e., 0.67 times shorter survival time compared to the baseline.

Also, on average, a 10 percentage increase of RM T cell proportion may delay the time to event by a factor of [(exp(0.08)]^10, which is 2.23, i.e., 2.23 times longer survival time compared to the baseline.

We plotted Kaplan-Meier curves (Figure 1) to compare the survival times among IPF subjects whose monocyte proportions are located at 1^st, 2^nd, 3^rd and 4^th quartiles. We plotted the Kaplan-Meier curves (Figure 2) also for RM T cells. As shown in the two plots, there is a clear separation of survival curves between 1^st quartile and 4^th quartile subjects (i.e., whose monocytes or RM T cells belong to 1^st or 4^th quartile in terms of cell proportions).

Open in a separate window

Figure 1:

Kaplan-Meier curve comparing 120 IPF patients with different monocyte proportions Note:1.Data come from dataset 1(GSE28042) and dataset 2 (GSE27957). 2. Log-RANK p value is for comparison between subjects with cell proportions at the 1^st or the 4^th quartiles, respectively.

Open in a separate window

Figure 2:

Kaplan-Meier curve comparing 120 IPF patients with different RM T cell proportions Note: 1.Data come from dataset 1(GSE28042) and dataset 2 (GSE27957). 2. Log-RANK p value is for comparison between subjects with cell proportions at the 1^st or the 4^th quartiles, respectively.

The Log-Rank test to compare the survival times between subjects whose monocyte proportions are at the 1^st vs. those whose monocyte proportions are at the 4^th quartile achieved a p value of 2.8E-04.

The Log-Rank test to compare the survival times between subjects whose RM T cell proportions are at the 1^st vs. those whose monocyte proportions are at the 4^th quartile achieved a p value of 6.1E-03.

According to Figure 1, the survival rate for subjects at ~2 years was >80% for subjects having the lowest monocyte proportion (1^st quartile of monocyte proportion, with a mean proportion of 33.4% and a range of 18.7% - 39.9%) whereas the survival was <40% for those subjects having the highest monocyte proportion (4^th quartile of monocyte proportion, with a mean proportion of 64.1% and a range of 58.3% - 77.1%).

Similarly (Figure 2), the survival rate at ~2 years was ~80% for those subjects having the highest RM T cell proportions (the 4^th quartile of RM T cell proportion, with a mean proportion of 17.0% and a range of 12.6% - 26.1%), while the survival was only ~50% for those subjects having the lowest RM T cell proportions (1^st quartile of RM T cell proportion, with a proportion of 0%).

Dataset 3:

GEO accession number: GSE93606 (8).

This is a transcriptomic study of whole blood from IPF subjects. The dataset contains 57 IPF (38 male and 19 female) subjects and 20 control (12 male and 8 female) subjects. The IPF subjects had a mean age of 67.4 years (50–84 years, sd = 8 years). The control subjects had a mean age of 66.0 years (48–83 years, sd = 10.6 years).

Time to outcome for those patients with observed outcome (i.e., death or decline in predicted Percent of Forced Vital Capacity (FVC) >10% over a six-month period) ranged from 2 to 32 months, with a mean of 16.15 months and a standard deviation of 9.40 months.

Cibersort analyses were performed on each subject using the LM22 signature file. Again, each sample achieved a p value of 0.000 in the Cibersort analyses. The results are shown in Appendix 1.3.

From some IPF subjects, whole blood samples were collected at multiple time points (i.e., months 0, 1, 3, 6 and 12). Among them, 13 subjects were collected at all 5 time points, 6 subjects were collected at 4 time points, 10 subjects were collected at 3 time points, 7 subjects at 2 time points. The remaining 21 subjects were collected at only one time at baseline (month 0). We therefore used the mean proportion across multiple time points for those IPF subjects whose blood samples were collected at multiple times.

Shown in Appendix 2, Table 5 are the summary statistics for the proportions (in percent) of different immune cell types.

As shown in the table, 3 cell types (follicular helper T cells, M1 macrophages and resting dendritic cells) are absent (0%) in all subjects and ten cell types (e.g., memory B cells) have a low variation (standard deviation < 1 percent) across all the subjects. This low variation is due to 0 or near 0 percent of proportion for the cell type for a large number of subjects. We therefore remove the above cell types (a total of 13) and keep the remaining nine cell types for downstream analyses.

The subjects also have DLco and FVC (Percent predicted Forced Vital Capacity) information, with a mean DLco of 39.2% and standard deviation of 14.1% and a mean FVC of 72.2% and a standard deviation of 20.3%.

Survival analyses were performed on these 57 IPF subjects using the survival data from each subject. At nominal significance level (alpha = 0.05) RM T cells (p = 3.26e-2) and neutrophils (p=3.92e-2) achieved statistical significance for association with IPF survival.

Kaplan-Meier analyses were performed to compare subjects belonging to the 4 quartiles of a RM T cell proportion (Figure 3). Again, RM T cell proportions in the blood have a positive association with IPF survival; the higher the proportion the longer the survival. For example, at ~15 months, the survival rate is ~85% for those subjects with the highest (4^th quartile, with a range of cell proportions of 4–11%) proportion of RM T cells in the blood, which is in contrast to the survival rate of ~45% for those subjects with the lowest (1^st quartile, with a range of cell proportions of 0–0.3%) proportion of the cells.

Open in a separate window

Figure 3:

Kaplan-Meier curve comparing 57 IPF patients with different RM T cell proportions Note: 1.Data come from dataset 3 (GSE93606). 2. Log-RANK p value is for comparison between subjects with cell proportions at the 1^st or the 4^th quartiles, respectively.

In contrast, the proportion of neutrophils in the blood has a negative association with IPF survival; the higher the proportion the worse the survival (Figure 4). For example, at ~15 months, the survival rate is ~80% for those subjects with the lowest (1^st quartile, with a range of cell proportions of 24–45%) proportion of neutrophils in blood, which compares to the survival rate of ~40% for those subjects with the highest (4^th quartile, with a range of cell proportions of 58–71%) proportion of the cells.

Open in a separate window

Figure 4:

Kaplan-Meier curve comparing 57 IPF patients with different neutrophil proportions Note: 1.Data come from dataset 3 (GSE93606). 2. Log-RANK p value is for comparison between subjects with cell proportions at the 1^st or the 4^th quartiles, respectively.

Note the above proportion for RM T cells is within whole blood, not within the PBMC compartment, as in the first two datasets. However, if we re-scaled the RM T cells’ proportion to that for the PBMCs compartment by dividing its proportion in blood with (1-neutrophil’s proportion), we still achieved a significant result for the survival analysis (p = 0.047).

Datasets 1, 2 and 3 combined:

To maximize the statistical power, we combined the Dataset 3 with the combined 120 subjects (Datasets 1 and 2).

Note Datasets 1 and 2 are from PBMCs yet Dataset 3 is from whole blood. The major difference between the PBMCs datasets (Datasets 1 and 2) and the whole blood dataset (Dataset 3) is neutrophils that on average makes up ~50% of cell proportion in whole blood (see Appendix 2, Table 5). Therefore, in the combined dataset, we did not include neutrophils.

Also, the cell proportions in Datasets 1 and 2 are in the scale of PBMCs, i.e., relative to all the cells in PBMCs. Yet the cell proportions in Dataset 3 is in the scale of whole blood, i.e., relative to all the leukocytes in the whole blood. To make the cell proportions comparable in terms of scale, we re-scaled the cell proportions in Dataset 3 (so that they will be in the scale of PBMCs) by dividing each cell’s proportion with (1-neutrophil’s proportion) for a sample.

As mentioned above, in Dataset 3 we removed those cells (13 in total) whose proportion’s standard deviation is <1%. Also, in the Datasets 1 and 2 combined, we removed 7 cell types whose proportion’s standard deviation is <1%. Therefore, in the Datasets 1, 2 and 3 combined, we only have 8 common cell types for analysis, which are naïve B cells, CD8 T cells, CD4 naïve T cells, CD4 resting memory T cells, CD4 activated memory T cells, resting NK cells and monocytes.

Lastly, we transformed the survival time in the scale of YTO in Datasets 1 and 2 to survival months (by multiplying 12) in order to make the time to outcome in the same scale of months across the three datasets.

The combined dataset has 177 subjects, among whom, 130 are males and 47 are females.

Among the 177 subjects, 92 (74 males and 18 females) have the observed IPF outcome (death, lung transplantation or a >10% decline of FVC over a six-month period) and 85 (56 males and 29 females) are censored.

The subjects have a mean age of 68 years, with a standard deviation of 8.1 years. The time to IPF outcome has a mean of 15.1 months and a standard deviation of 11.5 months. The time of observation for censored subjects has a mean of 27.3 months and a standard deviation of 9.4 months.

Using the combined data, we performed survival analysis (via a survival regression model, assuming exponential distribution of survival time) to correlate the survival months with a cell’s proportion, adjusting for age and gender.

RM T cell is the only cell type that achieved significant relationship (p = 2.65e-4, FDR = 2.12e-with the IPF survival time, with a survival regression coefficient of 0.0733. Under the exponential survival model, this coefficient may translate into a prediction that a 10 percent increase of RM T cell proportion may delay the time to IPF outcome by a factor of 2.08, i.e., (exp(0.0733))^10.

Cox proportional hazards regression analysis was also performed for the relationship between survival time with RM T cell proportions, adjusting age and gender. A p value of 3.67e-4 was achieved for the RM T cell proportions, with a coefficient of −0.072, which suggests that a 1 percent increase of RM T cell proportion may be associated with a hazard ratio of 0.9305 (exp(−0.072)), i.e., a decrease of hazard of ~7% for IPF. Correspondingly, a 10 percent increase of RM T cell proportion may be associated with a decrease of IPF hazard of >50% (= 1-(exp(−0.072))^10) for IPF.

Kaplan-Meier analyses were performed to compare survival rates for subjects belonging to the 4 ranks of a RM T cell proportions (Figure 5), where 1^st rank subjects have RM T cell proportions >/=0 and < 5%, 2^nd rank: >/= 5% and < 10%, 3^rd rank: >/= 10% and <15% and 4^th rank: >/= 15%.

Open in a separate window

Figure 5:

Kaplan-Meier curve comparing 177 IPF patients with different RM T cell proportions Note: 1.Data come from dataset 1(GSE28042), dataset 2 (GSE27957) and dataset 3 (GSE93606). 2. Log-RANK p value is for comparison of subjects with cell proportions at the 1^st vs. the 4^th ranks, respectively.

As shown in the figure, RM T cell proportions have a positive association with IPF survival; the higher the proportion the longer the survival. For example, at 30 months, the survival rate is ~78.7% for those subjects with the highest proportion of RM T cells (i.e., 4th rank, RM T cell proportions ranging from 15.03% to 26.11%), which is in contrast to the survival rate of ~33.5% for those subjects with the lowest proportion of the cells (i.e., 1st rank, RM T cell proportions ranging from 0 to 4.56%).

Dataset 4:

GEO accession number: GSE38958 (10)

This dataset contains gene expression analyses of PBMCs from 70 IPF patients (58 males vs. 12 females) vs. 45 healthy controls (27 males vs. 18 females).

Among these subjects, 60 IPF patients had age information and their mean age is 68.2 years. Thirty-five healthy controls had age information and their mean age is 69.3 years. The remaining subjects do not have age information.

The majority of the subjects are Caucasians (n =85). In addition, 17 subjects are African Americans, 6 subjects Hispanic Americans, 1 subject Asian American and 6 subjects with ethnicity unknown.

60 IPF patients had DLco information, with a mean of 43.3% (range: 30.8%−97%) and a standard deviation of 18.3%. 60 IPF patients also had FVC information, with a mean of 62.4% (range: 30%−92%) and standard deviation of 15.0%.

The results of our Cibersort analyses results of these 115 subjects are shown in Appendix 1.4. Again, for each of these 115 PBMC samples, the p values achieved in the Cibersort analysis (9) is 0.000. The summary statistics (in percent) for cell proportions inferred from Cibersort analysis (9)using LM22 signature file for the total 115 PBMCs samples are listed in Appendix 2, Table 6.

As shown in the table, the mean neutrophil percentage is 2.03, suggesting a minor granulocyte (neutrophil) contamination in the specimens during the PBMC isolation process. Histogram and stem-leaf plot analyses (not shown here) indicated that the majority (n>90) of samples have a neutrophil percentage value of 0, suggesting that the contamination is limited to only a few samples.

To adjust for the effects from the neutrophil contamination for some samples, we used the same method (as above) by re-calculating (re-scaling) each cell type’s proportion in each sample using the original proportion divided by (1-neutrophil’s proportion). Appendix 2 (Table 7) shows the summary statistics for the revised (re-scaled) cell proportions. Note three cell types, M1 macrophages, resting dendritic cells and activated mast cells, with 0 percent across all the samples were removed from further analyses and hence are not shown here. Also excluded is the neutrophil itself. As shown in the re-scaled table (Appendix 2, Table 7), there is little difference in the cell proportions from the original summary statistics table (Appendix 2, Table 6).

Using the above re-scaled cell proportions, we performed logistic regression analyses, modeling IPF status as the dependent variable and the cell proportions as independent variables, while adjusting for age, gender and ethnicity as covariates. Only monocytes (p= 0.0034) and naïve B cells (p = 0.0039) achieved significance, while RM T cells achieved a p value (p = 0.13) that is close to marginal significance (p = 0.10) and is the most significant after monocytes and naïve B cells.

Here, monocytes have a higher proportion in IPF (mean = 43.9%) vs. control subjects (mean = 37.3%) (Figure 6) and RM T cells have a lower proportion in IPF (mean = 14.2%) vs. control subjects (mean = 17.5%) (Figure 7). This direction of association of monocytes and RM T cells with IPF status is consistent with the survival analyses, where a higher proportion of monocytes was associated with a worse outcome (lower survival) yet a higher proportion of RM T cells associated with a better outcome (higher survival). As shown in Figures 6 and ​and7,7, the p values are 1.2E-4 and 6.1E-3 for comparing monocytes and RM T cells, respectively, using Wilcoxon rank test.

Open in a separate window

Figure 6:

Boxplot comparing monocyte proportions in IPF vs. control subjects Note: Data come from dataset 4 (GSE38958).

Open in a separate window

Figure 7:

Boxplot comparing RM T cell proportions in IPF vs. control subjects Note: Data come from dataset 4 (GSE38958).

We also examined the association of cell proportions with DLco, adjusting for age, gender and ethnicity. Again, RM T cells and monocytes achieved significant p values (p <0.05) for association with DLco. While the RM T cell proportion has a positive association with DLco (p = 0.023) (Figure 8), the proportion of monocytes has a negative association with DLco (p = 0.044) (Figure 9). Again, these data are consistent with the previous findings suggesting a better outcome for subjects with a higher proportion of RM T cells and a worse outcome for subjects with a higher proportion of monocytes.

Open in a separate window

Figure 8:

Scatterplot showing the linear relation between RM T cell proportion and DLco Note: 1. Data come from dataset 4 (GSE38958). 2. P value is based on a full regression model adjusted for age, gender and ethnicity

Open in a separate window

Figure 9:

Scatterplot showing the linear relation between monocyte proportion and DLco Note: 1. Data come from dataset 4 (GSE38958). 2. P value is based on a full regression model adjusted for age, gender and ethnicity

Discussion:

In this study, we took advantage of the genomics data of PBMCs or whole blood of IPF patients to enumerate proportions of 22 immune cell types in their blood. Through correlating the immune cell components with IPF outcomes, including survival time, IPF status and DLco, we found a consistent pattern of differential proportion of RM T cells and monocytes, which are associated with IPF survival, IPF status and DLco. In general, a higher proportion of RM T cells was observed in healthy controls than in IPF patients (Figure 6) and was associated with a better survival (Figures 2, ​,33 and ​and5)5) and a higher level of DLco in IPF patients (Figure 8). In contrast, a higher proportion of monocytes was observed in IPF than in control subjects (Figure 6) and was associated with a worse survival (Figure 1) and a lower level of DLco in IPF patients (Figure 9). Our data suggested a protective role for RM T cells and a detrimental role for monocytes in IPF.

Our findings are supported by studies using mouse models. For example, a recent study demonstrated the protective role of tissue-resident memory T cells (Trm) for IPF (11), where vaccine-induced lung Trm was found to reverse lung collagen, fibrocytes, and histologic injury and improve physiologic function in a mouse model of bleomycin-induced lung fibrosis. The study proposed a mechanism, where Trm in the lung may promote an immune microenvironment that can arrest and reverse the chronic processes in IPF.

For monocytes, their detrimental role to IPF has been suggested by a study, where CCR2+ infiltrating monocyte-derived macrophages may play a critical role in the development of radiation-induced pulmonary fibrosis (12). In another study, it was found that monocyte-derived alveolar macrophages may drive lung fibrosis by expressing profibrotic genes, while tissue-resident alveolar macrophages did not contribute to fibrosis (13). Both of the studies pointed to the potential importance of circulating monocytes as a source for the local macrophages that mediate the pathological process of IPF, which are again consistent with our findings of the unfavorable outcome (e.g., worse survival) for subjects with a higher proportion of circulating monocytes.

Our study also suggested neutrophils as another detrimental cell type for IPF as a higher proportion of neutrophils in the blood was associated with a lower survival rate (Figure 4). The detrimental role of neutrophils in IPF was suggested by a study, where high levels of neutrophil elastase were detected in the lung parenchyma and also in both the BALF (bronchoalveolar lavage fluid) and sera (14). Considerable numbers of neutrophils were also shown to infiltrate the lung parenchyma according to immunohistochemistry (14). A more recent study also demonstrated that neutrophil elastase may promote myofibroblast differentiation in lung fibrosis (15). Also, accumulation of neutrophils in the lungs of IPF patients was mediated by chemotactic factors released by alveolar macrophages (16).

Taken together, the findings from our study are supported by earlier evidence. As a potential usage to the field, RM T cells’ protective role should be considered in IPF treatment so that patients’ survival can be improved. For example, strategies that expand the population of memory T cells may be explored. In addition, cell and fraction counts of these key immune cell types, especially RM T cells and monocytes within the PBMC compartment and neutrophils within whole blood, may also represent useful tools for IPF prognosis and hence should be

closely monitored during the disease process. Future genomics studies of IPF may be focused more on these important cells, especially RM T cells and monocytes, so that the mechanistic roles of these specific cells in the disease process can be further elucidated.

Cibersort (9) is a software for cell proportion enumeration based on gene expression data. According to the study (9), the performance of the software has been validated by flow cytometry analysis. In particular, there was a significant correlation (p</=0.02) between Cibersort analysis and flow cytometry for RM T cells and monocytes as well as other immune cells, as discovered using blood samples from 27 adult subjects (9). Therefore, the cell proportions imputed using Cibersort in our study should be reliable, which is further supported by the extreme low p values (p = 0.000) for the cell proportions inferred for each sample (see Appendix 1). Moreover, our main purpose is to correlate the imputed cell proportions with survival outcomes, not to use the inferred cell proportions as a measure for disease diagnosis. Since the imputed cell proportions have a significant correlation with flow cytometry results, it is reasonable to project that our findings may not be significantly different from those achieved using “real” flow cytometry data.

Nevertheless, as a key limitation, our findings were based on in silico genomics data mainly of cross-sectional designs. It is necessary to further validate our findings using a longitudinal study design (ideally with flow cytometry data) or with mechanistic experiments on model animals.

Methods:

All of the analyses were performed using basic R functions and specific R packages for statistical analysis.

Specifically, we downloaded the four microarray datasets (i.e., the Series Matrix Files) from GEO (Gene Expression Omnibus) website (https://www.ncbi.nlm.nih.gov/geo/). The GEO accession numbers of the datasets are provided in the Results section. The gene expression count matrix was annotated for gene symbols using the platform files provided in GEO. For example, for the GSE28042 dataset, the platform file used for annotation is GPL6480.

LM22 is the signature genes file we used for Cibersort analyses (9). The file contains expression counts for 547 signature genes (547 rows) for 22 distinct human immune cells (22 columns). The annotated gene expression count matrix file (as the gene expression mixture file), together with the LM22 signature file, were the input files for Cibersort analyses (9).

The technical details of Cibersort analyses are provided in its website https://cibersort.stanford.edu/index.php as well as the published paper (9). In principle, the method involves a deconvolution model m=f x B, where m is the gene expression mixture file, B a gene expression profile (GEP) signature matrix and f a vector consisting of the unknown fractions of each cell type in the mixture. The process is essentially solving f through this system of linear equations based on information in m and B. The signature matrix, B, is built by expression data of purified or enriched cell populations. Specifically, for the LM22 signature matrix, GEP data was obtained from the public domain for 22 leukocyte subsets profiled on the HGU133A platform (9).

The resultant file from Cibersort analysis (9) is a matrix with rows corresponding to input samples and columns corresponding to the 22 specific cell types. Please see Appendix 1 for details. Inside the matrix are the proportions of a specific cell type in a specific sample. The resultant file also contains a p value for each sample, which is statistical significance of the deconvolution result across all cells for a sample.

We exported each resultant file from Cibersort analysis (9) as a .csv file and then merged the file by the sample name with the design matrix of a study. The design matrix was extracted from the top section of a Series Matrix File downloaded from GEO. The merged file has rows corresponding to specific samples from human patients and columns corresponding to the variables, including patients’ characteristics, e.g., age, gender, race, IPF status, DLco, etc., as well as cell proportions for the 22 specific immune cells. This merged file was used in survival and other downstream analyses.

For survival analyses, we used R package, Survival (17, 18), and implemented the “survreg” function, assuming exponential distribution for the survival time. The Cox proportional hazard regression analysis was also performed using the Survival package with the “coxph” function. The survival time was modeled as a function that is determined by a cell’s proportion as well as other key covariates, such as age and gender.

For the Kaplan-Meier curve analyses (Figures 1–5) to compare subjects with different quartiles (or ranks) of a cell’s proportion, the Survfit function in the Survival package (17, 18) was used.

We also used “pairwise_survdiff” function of the “survminer” package to compare subjects in terms of survival whose cell proportions belong to the 1^st or the 4^th quartile (or rank) (Figures 1–5) in order to infer the Log-RANK p values.

To compare a cell’s proportion between IPF and normal control subjects while controlling for other key covariates, such as age, gender and race, we used R’s glm function that modeled IPF status as the dependent variable (using logit link function in binomial family) and a cell’s proportion and the covariates as independent variables. We also used R’s lm function to model relationship between DLco and a cell’s proportion, controlling for covariates, such as age, gender and ethnicity.

Box plots in Figures 6 and ​and77 were drawn using R package ggpubr. Scatterplots in Figures 8 and ​and99 were drawn using R package ggplots (19).

Supplementary Material

Acknowledgement

Footnotes

Literature Cited