Radiomic Profiling of Lung CT in a Cohort of Sarcoidosis Cases

Introduction

Sarcoidosis is a granulomatous interstitial lung disease which affects ~ 110 thousand individuals in the United States (1). Typical diagnoses is between 30-50 years of life, resulting in a decrease in quality of life and productivity (2). Pulmonary disease occurs in over 90% of those with sarcoidosis (3) with significant morbidity and mortality. Currently, visual assessment of chest radiography (CXR) is used to quantify lung abnormalities standardized via the Scadding staging system (4) which groups CXR findings into five groups from 0 to 4. Substantial variation exists in the radiographic patterns within each Scadding stage. Scadding stage has limited utility in predicting prognosis even in the extremes of the scale (5).

Chest high resolution computed tomography (HRCT) is increasingly used in clinical practice of sarcoidosis to monitor disease as it offers more detailed visualization of parenchymal abnormalities compared to CXR (4, 6–8). As with CXR, visual assessment of chest HRCT is used to evaluate abnormalities although there are limited standardized scoring metrics (4, 9, 10). This is due in part to the diverse and heterogeneous patterns present on chest HRCT in sarcoidosis, often with multiple patterns noted on one CT. These complexities result in high inter- and intra-rater variation (11). More automated systems that quantify sarcoidosis chest HRCT could decrease these sources of variation.

Radiomics is a field of study in which large numbers of quantitative features are extracted from medical images (12). A radiomics panel computes summary measures of the distribution of the Hounsfield units (HU) along with summary measures of the spatial relationships of neighbouring voxels (13). The result is a panel of quantitative measures characterizing the texture of the image.

Radiomic analysis has proved useful for quantifying HRCT in emphysema (7), idiopathic pulmonary fibrosis (14, 15), interstitial lung disease (16, 17), diffuse lung disease (18) and cancer (19). Ryan et al. (20) showed the potential utility of radiomics in sarcoidosis, comparing radiomic and other textural based measures between sarcoidosis patients and controls. It remains unclear whether radiomics also has the potential to differentiate varied phenotypes within sarcoidosis patients. Radiomic profiles within sarcoidosis patients that also correlate with pulmonary function or patient reported outcomes could indicate that radiomics may serve as a useful measure to track change in the lung parenchyma over time.

The goal of this research study is to develop a radiomic profile of sarcoidosis chest HRCT using a large, phenotypically-diverse population of sarcoidosis cases from the Genomic Research in Alpha-1 Antitrypsin Deficiency and Sarcoidosis (GRADS) study (21) using statistical clustering techniques and to investigate the clinical utility of the clusters by quantifying their association with pulmonary function testing (PFT) and patient reported outcomes (PROs). PFTs and PROs were chosen to capture both clinician and patient focused measures of disease (26,,27).

Methods

Radiomic Analysis

Details of the imaging acquisition and lung segmentation can be found in the online supplement. Radiomic features were calculated on the left and right lungs using RIA_lung from the lungct R package. These features include 44 first-order features and 239 GLCM features for each lung, for a total of 566 features. To calculate the GLCM features, the Hounsfield units (HU) from each HRCT scan were first discretized into 16 bins with equal relative frequencies; then, the features were calculated in each of 26 directions, assuming a voxel distance of one; finally, these features from all directions were summarized using the mean statistic. The exact computation for each radiomic measure can be found elsewhere (31).

The ComBat function (https://github.com/Jfortin1/ComBatHarmonization) was used to harmonize the radiomic measures across scanners (32–34), while preserving the biological variability in the following covariates: age, height, BMI, sex, race and clinical phenotype.

As radiomic features are high-dimensional and repetitive (Figure 1), we first used a decorrelation filter (see online supplement), to identify a subset of features for analysis. This reduced the number of variables from 566 to 97. We then used robust k-means (35) to cluster patients (R package, RSKC). For outlier robustness, the proportion of cases trimmed was set at 0.1; the optimal bound on feature weights was found to be 7.5, using the permutation approach observed in sparcl R package (36). We used a standardized Gap statistic to determine the optimal number of clusters using the cluster R package (37). A standardized Gap statistic is the typical gap statistic divided by its standard deviation. The variables with the top five weights from the RSKC clustering algorithm were then selected for further investigation (discriminative features).

• Download figure
• Open in new tab

Figure 1:

Heat map of the correlation between different radiomic measures for the entire population.

Results

Tables 1 and 2 show the characteristics of the study population, which included N=321 sarcoidosis cases, including 147 (45.8%) males and 233 (72.6%) whites, with an average age of 53 (SD=10) years, average height of 67.0 in (SD=4.2) and average BMI of 30.6 kg/m² (SD=6.5). Participants were spread across Scadding stages, with 43 (13.5%) in stage 0, 63 (19.8%) in stage 1, 92 (28.9%) in stage 2, 44 (13.8%) in stage 3 and 76 (23.9%) in stage 4. The average FVC was 2.62 L (SD=0.93), FEV1 was 3.57 L (SD=1.14), and DLCO 80.15 (SD=23.37). The population was predominantly of the non-obstructive type (72.7%).

View this table:

• View inline
• View popup

Table 1:

Patient demographics by radiomic grouping, ordered from least severe (0) to most severe (4) based on average FVC. Unless otherwise noted values are mean (SD).

View this table:

• View inline
• View popup

Table 2:

Summary measures of pulmonary function testing and self-reported outcomes by radiomic grouping

Radiomic Based Clustering of Sarcoidosis

Five groups (clusters) of patients were identified (Table 1). Ordered from highest FVC to lowest, 75 (23.4%) patients were in radiomic group 0, 57 (17.8%) in group 1, 51 (15.9%) in group 2, 58 (18.1%) in group 3 and 80 (24.9%) in group 4. The new radiomic-based grouping was associated with Scadding stage (p<0.001), but not a direct reflection of Scadding stage (Table 1).

PFT differed between radiomic groups (Figure 2; p<0.0001). For FVC, radiomic groups 3 and 4 had between 0.5 and 0.7 L lower average FVC compared to radiomic groups 0, 1, and 2. For FEV1, radiomic groups 2 and 3 and had an average FEV1 0.5 L lower than radiomic groups 0 and 1, which were similar. Radiomic group 4 had the lowest average FEV1 and was approximately 0.8 L lower than groups 0 and 1. DLCO has a similar pattern to FVC. Groups 2 and 4 demonstrated more obstructive disease compared to groups 0, 1, and 3 (57.4% and 38.0% vs. 13.5%-17.2%, respectively; p<0.001).

• Download figure
• Open in new tab

Figure 2:

Group differences in average PFT values from regression models with adjustment for demographics (circles) and then also Scadding stage (triangle). The significance of the radiomic group differences from group 0 are in green (p>0.05) and blue (p<0.05). The bottom left text in each panel is the overall p-value for the association between radiomic grouping and PFT outcome along with the R². The bottom right text if the overall p-value for the association between radiomic grouping and PFT after additional adjustment for Scadding stage along with the R².

Some PROs also differed between groups, but not consistently (Figure 3). Average shortness of breath (SOBQ) score differed between radiomic groups (p<0.0001). Radiomic group 4 had higher average SOBQ compared to radiomic groups 0 and 1 with radiomic group 4 being approximately 17.5 units higher than radiomic groups 0 and 1. In addition, average physical health (SF-12) was lowest for radiomic group 4 and linearly increased to radiomic group 0.

• Download figure
• Open in new tab

Figure 3:

Group differences in average PRO values from regression models with adjustment for demographics (circles) and then also Scadding stage (triangle). The significance of the radiomic group differences from group 0 are in green (p>0.05) and blue (p<0.05). The bottom left text in each panel is the overall p-value for the association between radiomic grouping and PFT outcome along with the R². The bottom right text if the overall p-value for the association between radiomic grouping and PFT after additional adjustment for Scadding stage along with the R².

The new radiomic groups remained significantly associated with PFT after adjusting for Scadding stage (Figure 2; p<0.0001). Scadding stage also remains significant for all PFT after controlling for radiomic group (p<0.0041; Figure 2). The radiomic group also remained associated with SOBQ score (p=0.032; Figure 3), while Scadding stage was no longer associated with SOBQ (p=0.072; Figure 3). None of the other PROs maintained a significant association with radiomic group after adjustment for Scadding stage (p>0.38).

The five most discriminatory radiomic features included kurtosis, which is a measure of shape of the distribution of the HU from an image, as well as four summary measures from the gray level co-occurrence matrix (GLCM). The GLCM measure spatial correlation and similarity of the HU in image voxels near each other. Image visualization of cases with different values of kurtosis and two GLCM measures are shown in Figure 4 with distribution of kurtosis noted in Figure 5 for maximum, median and minimum kurtosis levels in our population. The images visibly show the increase in observable parenchymal abnormalities (left panel Figure 4) with lower kurtosis.

• Download figure
• Open in new tab

Figure 4:

CT images in axial orientation for three patients with minimum, median, and maximum values for GLCM Gaussian (left column), the GLCM Inverse Gaussian, and the kurtosis.

• Download figure
• Open in new tab

Figure 5:

Distribution of HU from HRCT images with the maximum (32.2) (solid), median (9.6) (dashed) and minimum (0.4) (dotted) percentiles of the kurtosis distribution, respectively.

The discriminatory radiomic measures were jointly associated with FVC, FEV1, FEV1/FVC and DLCO (p<0.001; Table 3). The radiomic measures explained between 9 and 18% more variation in PFT than adjustment for demographics (age, race, sex, height and BMI) only. For comparison, Scadding stage explained between 0 and 11% more variation in PFT. Kurtosis was associated with FVC, FEV1/FVC and DLCO (p<0.0001). A 1-unit increase in kurtosis was associated with a 0.4 L (SE=0.07; p<0.0001) increase in FVC. Kurtosis was also positively associated with increasing DLCO (1.59, SE=0.68, p=0.020). Kurtosis was negatively associated with FEV1/FVC. A 1-unit increase in kurtosis was associated with a 7% (SE=0.01, p<0.0001) decrease in FEV1/FVC ratio.

View this table:

• View inline
• View popup
• Download powerpoint

Table 3:

Results of the regression analysis of the five discriminatory radiomic measures for PFT.^d

Validation of Clustering

Pairwise ARI values ranged from 0.3 to 1 and peaking around 0.5 (online supplement Figure 2 left). The maximum ARI of fit clusters with Scadding stage was 0.07. In the linear models, the maximum p-values for the significant of cluster in the demographic adjusted and demographic plus Scadding models were <0.0001 for both models. The corresponding proportions of p-values less than 0.01 were 100% for both models.

Bootstrap samples contained between 184 and 219 unique observations and contained 203 on average. Pairs of bootstrap samples contained between 100 and 159 unique overlapping observations and contained 129 on average or about 63% of the sample is used in each ARI calculation. Pairwise ARI values computed from unique overlapping observations ranged from 0.2 to 1 and had distribution peaking around 0.5 (online supplement Figure 2 right). The maximum ARI of fit clusters with Scadding stage was 0.17. For significance of cluster label in linear models fit to bootstrap samples, maximum p-values for the significant of cluster in the demographic adjusted and demographic plus Scadding models were <0.0001 for both models. The corresponding proportions of p-values less than 0.01 were 100%.

Taken together the validation analysis suggest fairly mild sensitivity of clustering to the analysis pipeline. There is more sensitivity to bootstrapped samples; however cluster groups remain a significant predictor of FVC across all sources of random perturbations in the pipeline.

Discussion

Radiographic manifestations in sarcoidosis are protean. As a result, image characterization is traditionally done visually and not standardized. The utility of visual assessment in sarcoidosis is limited by the intra-observer and inter-observer variability. Radiomics is a more reproducible and computationally efficient approach to characterize HRCT. We used radiomics to quantitatively characterize images from a large, phenotypically diverse cohort of sarcoidosis subjects and demonstrated that radiomics are associated with clinical and patient reported outcomes of disease. Using a common unsupervised learning approach (35), we identified five clusters of cases based on radiomic characterization of their chest CT. These five radiomic-based groups differed significantly by PFTs and several PROs. Notably, each radiomic group included a range of Scadding stages and radiomic group remained significantly associated with PFT after adjusting for Scadding stage. These data suggest that radiomics represents radiographic abnormalities that differ from Scadding stage.

Demographic characteristics of the individual explained a sizable amount of the variation observed in PFT (R²=50-60%). Radiomics also explained a significant amount of additional variation in PFT (8% to 18%). This additional amount of variation explained is consistent with other quantitative imaging approaches such as CALIPER (16) and those used for investigations of other lung conditions such as systemic sclerosis (38) and diffuse interstitial lung diseases (39).

Except in lung cancer research, much of the work on quantitative CT image analysis has focused only on the HU density (17, 38). Our radiomic panel included summarization of both the density and the spatial characterization. We note that more measures of the GLCM appear as discriminative in the cluster analysis than first-order densitometry measures. This implies that the texture (the spatial information) is perhaps more useful for differentiating various parenchymal abnormalities seen in sarcoidosis.

Although it remains speculative to the reasons for the association between decreased kurtosis and decreased FVC and DLCO, the following explanation is plausible. As kurtosis decreases, more severe parenchymal abnormalities are present that impact the function of the lung, as is shown in Figure 4. Severe parenchymal abnormalities in sarcoidosis are observed as higher HU values. In Figure 5 we observe the presence of a higher frequency of moderately high HU values dampening the peak in the HU distribution and leading to a lower kurtosis in this case. We also found a negative association between lung function and Gaussian weighted GLCM. Figure 4 highlights more abnormalities present with a higher Gaussian weighted GLCM value, as might be expected with this association.

The associations with PROs were less consistent, although important patterns emerged. Radiomic group 4 demonstrated decreased quality of life based on physical function in the SF-12 and worse shortness of breath on the SOBQ, findings that are consistent with the worse lung function in this group. Radiomic group explained substantially less variation in PROs (1 to 15%) compared to PFTs. While the tools we used to measure PROs are validated and used widely, they are not generally correlated with objective measures of lung function (40). As an example, fatigue is a known multi-factorial PRO in sarcoidosis and may impact a number of other PRO (22, 40). For example, fatigue can be associated with shortness of breath and the SOBQ but also with other organ involvement; these two measurements can be confounded by physical function and cardiac or even neurological involvement. Our finding that radiomic group is associated with decreased physical function and increased shortness of breath provides encouragement that more detailed characterization of lung abnormalities, like we developed here, could contribute to a better understanding of the current disconnect between objective measures of lung function and patient experiences.

This study had several strengths. To our knowledge, this was the largest group of sarcoidosis patients with research grade HRCT available allowing for a full 3-D based quantitative analysis. In addition to the high-quality and robust clinical data, the GRADS study also provides a comprehensive and consistently collected set of patient reported outcomes, which reflect important aspects of treatment decision-making processes (23). The results were also internally validated. The validation results taken together suggest our approach is not over fitting to these data.

This study is not without limitations. The GRADS study relied on enrolment from academic centres and could be skewed to a population that was referred for worse disease, was near one of the centres, which were primarily localized to the Eastern US. The demographic is predominantly white and of higher SES as a result and thus did not provide full representation of the spectrum of disease. In addition, many subjects had disease for a decade or longer and were treated. While the same protocol was used across study sites for obtaining the CT images we studied, the scanners themselves differed. We performed harmonization of the scans to mitigate the potential for systematic large differences in scans due to scanner type. In addition, the distribution of PFT measures is not dependent on scanner type such that we expect any differences in the radiomic measures due to scanner type are non-differential as they relate to PFTs. Finally, the number of radiomic clusters (groups) we identified may not directly translate to other populations of sarcoidosis patients or even be the optimal solution for this group of patients. The optimization algorithm we used is common, however, there is no single way to choose the optimal number of clusters.

Radiomic and other quantitative imaging approaches have a major strength in that they can be computed in a time efficient and reliable, reproducible automated procedure. The radiomics package in R took only 3 minutes to compute all the 3-D radiomics measures. This means large quantities of scans can have radiomic profiles computed for aiding in decisions on what scans the radiologist might prioritize for further consideration for visual or clinician evaluation.

In summary, this work provides evidence that in sarcoidosis radiomic quantification is useful for classifying abnormality of the lung along with pulmonary function and to a lesser degree patient reported outcomes. Future work should evaluate the potential of radiomics to capture small changes over time not easily detected by radiologic assessment to further evaluate the potential of radiomics to serve as a clinically useful quantitative measure of sarcoidosis presentation and progression.