Quantitative CT and COPD: cluster analysis reveals five distinct subtypes with varying exacerbation risks

Abstract

Background

The heterogeneity of chronic obstructive pulmonary disease (COPD) is increasingly recognized. To characterize the heterogeneity of COPD, we aimed to identify subtypes related to quantitative CT by using principal component analysis (PCA) and cluster analysis.

Methods

The data of 1879 participants in the SPIROMICS study were obtained from the NHLBI Biologic Specimen and Data Repository Information Coordinating Center. A combination of PCA and k-means clustering was used to analyze the data from these participants in the SPIROMICS study. We randomly split the samples into training and validation sets. Clusters were evaluated for their relationship with acute exacerbation risk throughout the entire follow-up period. The results of the training set were confirmed in the validation set. To avoid sampling errors, we conducted 10 random sampling cycles. Normalized mutual information (NMI) was applied in every cycle to evaluate the stability of clustering.

Results

We identified five clusters related to quantitative CT characterized as follows: (1) male-dominated low disease impact cluster, (2) obesity with relatively high symptom burden cluster, (3) airway wall lesion cluster, (4) lung upper region zone-predominant emphysema cluster, (5) severe emphysema cluster. There are significant differences in acute exacerbation risk among these five clusters.

Conclusions

Cluster analysis identified 5 clusters related to quantitative CT of all participants in the SPIROMICS cohort with significant differences in baseline characteristics and acute exacerbation risk. The stability of clustering results was validated through NMI in 10 sampling cycles. In addition, dimensionality reduction results showed high reproducibility in different studies.

Background

Chronic obstructive pulmonary disease (COPD) is now one of the top three causes of death worldwide and most of these deaths occur in low-income and middle-income countries [1, 2]. At present, the diagnosis of COPD, which is most widely used, is confirmed based on FEV1/FVC, and the severity of the disease is classified based on the FEV1% predicted [3]. However, GOLD 2024 describes COPD as a heterogeneous lung condition characterized by chronic respiratory symptoms resulting from abnormalities of the airways and/or alveoli [4]. This definition emphasizes the heterogeneity of COPD and suggests that abnormal lung changes in early COPD are a continuous process that may also exist in subjects who have not yet experienced standard airflow obstruction, such as the Pre COPD or PRISM population [5,6,7,8]. In fact, this method of diagnosis and classification based on the degree of airflow obstruction has led to improved diagnosis and treatment of the disease [9, 10], but it can also lead to significant overlap between the different disease features and the proposed subtypes, which may not reflect phenotypic heterogeneity. With the increasing availability of CT imaging, CT-based measurements are becoming a reliable and objective method for assessing the risk of acute exacerbations in COPD [11]. Therefore, while maintaining the diagnostic approach for COPD, researchers have continued to verify multiple subtypes by integrating CT images with clinical variables [12,13,14], which have important consequences for clinical management, such as asthma-COPD overlap [15] and upper lobe-predominant emphysema [16].

With the increasing number of measurements related to COPD, researchers have applied unsupervised algorithms to discover COPD-related subtypes [17]. The advantage of this approach is that it can uncover complex relationships among diverse variables, thereby identifying potential subtypes. The two main types of unsupervised machine learning algorithms are the clustering and dimensionality reduction algorithms [12, 18]. Disease axes, which are generated by specific algorithms such as dimensionality reduction, represent continuous trajectories of disease phenotypic features [17, 19]. Unlike discrete clusters that group individuals into distinct subtypes, disease axes provide continuous measures that are composed of many contributing variables, making them more suitable for situations where the variable set follows a continuous distribution [18, 20]. Principal component analysis (PCA) is a widely used technique for reducing dimensionality by generating linear combinations of features that maximally explain the observed population variance and thus can be applied to define disease axes. However, because unsupervised algorithms are purely data-driven and do not rely on predefined labels, the reproducibility and stability of the obtained clusters or disease axes in different studies have not been extensively recognized.

Therefore, to identify subtypes related to quantitative CT and explore the reproducibility of disease axes and clusters, we employed principal component analysis (PCA) and clustering on clinical data obtained from the SPIROMICS cohort. Additionally, we assessed the reproducibility of the obtained disease axes and clusters.

Methods

Research data

The data analysis in this study utilized the open-access SPIROMICS (Subpopulations and Intermediate Markers in COPD Study) (Clinical Trial Registry: ClinicalTrials.gov, Identifier: NCT01969344, registered on October 25, 2013) dataset from the NHLBI Biologic Specimen and Data Repository Information Coordinating Center. The SPIROMICS study has previously been described in detail [21]. Briefly, between November 2011 and January 2015, 2982 participants, aged 40–80 years, were enrolled in a multicenter prospective observational study funded by the NIH [21] and distributed across four enrollment strata: never-smokers, smokers without COPD, smokers with mild or moderate COPD, and smokers with severe COPD. We excluded individuals with missing data (Supplementary Data: Fig. S1). The median follow-up period for the 1879 subjects in this study was 1611 days.

Sample splitting

To evaluate the effectiveness of clustering solutions, the SPIROMICS data were randomly split into equally sized training and validation sets. All subsequent model building was performed on the training data. The validation set was used to validate clustering characteristics.

Dimensionality reduction

To understand the phenotypic spectrum of COPD, we applied PCA to standardized individual feature, spirometric data and CT quantitative data from the SPIROMICS training set. (Supplementary Data: Table S1). We aimed to include enough principal components in the clustering analysis to account for a sufficient proportion of the variance. However, since the explained variance ratios for principal component 9 and those thereafter were similar and relatively low, we chose to retain the first eight principal components. (Supplementary Data: Fig. S2)

To better interpret the principal components, we applied varimax rotation to the first 8 principal components using the method of maximizing variance. Varimax rotation does not alter the relative positions of data points, and therefore does not affect the results of k-means clustering. The principal component scores were calculated by multiplying the rotated principal component weights with the standardized original variables and summing.

Identifying clusters

We chose the k-means clustering method for this research, because it allows to verify the clustering results in the validation set by using the centers of the identified clusters. Then, we selected k = 5 as the optimal number of clusters, based on the greatest differentiation in clinical characteristics and acute exacerbation risk among the clusters. We used logistic regression to calculate the odds ratios (OR) of exacerbation for each subtype. To validate the significance of the subtyping, we adjusted the GOLD stage and calculated the OR of exacerbation for each cluster. Cox regression analysis was performed to estimate risk of exacerbation and adjusted for race.

Clusters validation

We computed the principal component scores for the validation set using the principal component score formulas derived from the training set. Based on the centers learned by the k-means algorithm in the training sample, we assigned the validation set to the nearest cluster center and allocated clusters in the validation sample accordingly. We utilized Wilcoxon rank sum tests to examine the differences in clustering features between the training and validation samples. K-means clustering was performed using the k-means function. The PCA was conducted using the psych package in R language. All statistical methods were executed in R language 4.3.2. Overview of methods is shown in Fig. 1.

Fig. 1

Overview of clustering, principal component analysis, and validation

Cross-validation

To avoid sampling errors, we conducted 10 random sampling cycles. In each cycle, cross validation was applied to each cycle (Details in Supplementary Data). We used normalized mutual information (NMI) to compare the two clustering solutions for each independent dataset.

Results

The characteristics of the training and validation samples are shown in Supplementary Data: Table S2, demonstrating their comparability.

The PCA with varimax rotation identified 8 principal components, explaining 73% of the variance in all variables. Detailed PCA loadings are shown in Fig. 2. PCA loadings represent the correlation between the original variables and the principal components, indicating how much each variable contributes to each component. Higher loadings reflect stronger contributions. We defined disease axes using these principal component scores.

Fig. 2

Principal component analysis description

Rotated principal component 1 (RC1), high representing by quantitative CT measurements of emphysema and air trapping, was interpreted as representing a multidimensional air trapping disease axis. Rotated principal component 2 (RC2), high representing by lung function measurement representation, was interpreted as a negative lung ventilation function disease axis. Rotated principal component 3 (RC3), mostly correlated with airway wall area percentages (loading scores 0.44 to 0.70), was interpreted as an airway wall lesion disease axis. Rotated principal component 4 (RC4), mostly correlated with high Pi10 index (loading scores 0.62 to 0.96), was interpreted as an airway wall thickness disease axis. Pi10 was calculated by regressing the square-root wall area on internal perimeter of included airways to predict the square-root wall area of a single hypothetical airway with internal perimeter of 10 mm. Notably, airway wall area percentages and Pi10 appear as independent components in the principal component analysis, representing principal component 3 and principal component 4, respectively, with low correlation between the two variables. This indicates that changes in Pi10 and airway wall area percentages are not linearly related in the progression of COPD. Rotated principal component 5 (RC5) to rotated principal component 8 (RC8) were not explained due to their low variance contribution.

Through PCA and cluster analysis, we identified five clusters and ranked them according to the acute exacerbation rate, with the lowest exacerbation rate cluster labeled as cluster 1.

Cluster 1: male-dominated low disease impact

Cluster 1 represents 26% of the SPIROMICS training sample and is characterised by the highest predicted FVC% and FEV1%, the lowest CAT score and mainly composed of males (Male individuals account for 80% of this cluster) (Table 1). The majority of individuals in cluster 1 is primarily composed of control subjects and individuals classified as GOLD stage 0 and GOLD stage 1. (Table 2)

Table 1 Cluster characteristics in training set, k = 5

Table 2 The composition of COPD stages and exacerbation risk for each cluster

In terms of quantitative CT data, this cluster exhibits the highest segmental airway wall area, the highest lumen area and total bronchial area compared to other clusters, while simultaneously having the lowest airway wall area percentage.

Cluster 2: obesity with relatively high symptom burden

Cluster 2 represents 14% of the SPIROMICS training sample and characterised by the highest BMI among the 5 clusters and a relatively high burden of COPD symptoms (CAT score greater than 10) (Table 1). This cluster mainly includes members of GOLD stage 0, GOLD stage 2 and the control group. (Table 2)

Compared to cluster 1, this cluster has a decrease in FEV1% predicted by 0.09. However, due to a decrease in predicted FVC% predicted of 0.14, the FEV1/FVC ratio of cluster 2 increases by 0.05. Additionally, this cluster has the lowest level of emphysema among the 5 clusters, with the skewness and kurtosis of the lung density histogram being the lowest among the 5 clusters, which may suggest a possible association with pulmonary fibrosis [22, 23]. Female individuals account for 67% of this cluster.

Cluster 3: airway wall lesion

Cluster 3 represents 37% of the SPIROMICS training sample and characterised by the lowest airway wall area among the five clusters. However, the percentage of airway wall area in this cluster is relatively high, suggesting that the increase rate of airway wall area towards the airway lumen exceeds the decrease rate of airway wall area, resulting in a net increase in wall area percent (Table 1). This cluster mainly includes members of GOLD stage 0, GOLD stage 1, GOLD stage 2 and GOLD stage 3 (Table 2). Females represent 59% of this cluster.

Cluster 4: lung upper region zone-predominant emphysema

Cluster 4 represents 12% of the SPIROMICS training sample and characterised by emphysema with marked upper zone-predominance. Air trapping in the upper regions are more severe compared to the lower regions, with the area of air trapping in the upper regions being more than twice that of the lower regions. Furthermore, this cluster exhibits the highest standard deviation in CT values, which may be related to the predominance of air trapping in the upper lung area of the cluster (Table 1). This cluster mainly includes members of GOLD stage 2 and GOLD stage 3. Females represent 50% of this cluster.

Cluster 5: severe emphysema

This cluster represents 12% of the training sample and is characterized by severe emphysema, airflow obstruction and the highest CAT scores. This cluster has the lowest BMI among the clusters and predominantly consists of males (Table 1). This cluster mainly includes members of GOLD stage 2, GOLD stage 3 and GOLD stage 4. Female individuals account for 32% of this cluster.

Validation of the clusters

All displayed baseline characteristics of the clusters have undergone Wilcoxon rank sum tests between the training and validation sets. Values marked with * indicate Wilcoxon rank sum tests p-values less than 0.05. The similarity of clustering features between the training and validation samples suggests that clustering can be reliably reproduced in separate data samples. The clustering characteristics of the validation set can be seen in Supplementary Data: Table S3. In cross validation, the results showed that the median NMI value was 0.66 for both the training and validation sets, which confirms the reproducibility of clustering results. (Supplementary Data: Fig. S3)

Acute exacerbation risk of each cluster

The exacerbation rate across the follow-up period in each cluster is shown in Table 2. Unless stated otherwise, the exacerbations mentioned in this study include mild, moderate, and severe exacerbation. Using cluster 1 as a reference, the hazards ratio and OR of acute exacerbation for other clusters are shown in Table 2.

To determine whether the association observed with these clusters and acute exacerbation risk was driven by severity of airflow obstruction, we repeated the cluster association tests adjusting for GOLD stage. The associations with exacerbation remained significant (p < 0.01). This suggests that the discovered clusters provide information on the risk of acute exacerbation independent from COPD severity as defined by GOLD.

Discussion

By employing PCA and cluster analysis on clinical data obtained from the SPIROMICS cohort, we distinguished 5 clusters as follows: (1) male-dominated low disease impact cluster, (2) obesity with relative high symptom burden cluster, (3) airway wall lesion cluster, (4) lung upper region zone-predominant emphysema, (5) severe emphysema cluster. The clinical characteristics and acute exacerbation risk of these clusters were validated through validation set.

This analysis reveals novel insights into COPD subtypes, particularly describing two clusters: the male-dominated low disease impact cluster and the airway wall lesion cluster, which have not been extensively described in prior research. Several factors contributed to the identification of these subtypes: (1) this study includes the most comprehensive set of clinical variables compared to previous COPD clustering studies, encompassing lung function tests, emphysema, airway wall thickness, and individual characteristic variables. (2) our study population also includes subjects with no airflow obstruction. Our work also addresses the issue of reproducibility in clustering analysis across independent samples derived from the same cohort. Furthermore, we evaluated whether the correlation of clinical characteristics within the clusters matched those observed in previous studies, thereby enhancing the credibility of the clustering.

This study offers a novel insight into airway wall progression in COPD, suggesting that airway wall thickness alone may not accurately reflect airway pathology. In the male-dominated low disease impact cluster, which has the highest airway wall area, the exacerbation rate during follow-up was the lowest. Conversely, in the airway wall lesion cluster, despite having the lowest airway wall area, it also has the smallest total bronchial wall and lumen areas and a higher exacerbation rate than clusters 1 and 2. This indicates that COPD progression may involve not only thickening of the airway wall towards the lumen but also destruction of the airway wall, manifested as a decrease in total bronchial wall area.

The marked gender differences between clusters are of particular interest because gender was not included in the dimension reduction process in machine learning as an original variable. However, the variability caused by gender differences was captured by other selected variables, such as FEV1, indicators related to airway wall thickness and so on [24, 25].

This study further validates previous discoveries concerning subtypes. Firstly, Peter J. Castaldi et al. identified a subtype characterized by emphysema in the upper lung region using k-means clustering in the COPDGene cohort, which closely aligns with our cluster 4. However, our cluster 4 in the training set exhibits higher airflow obstruction levels compared to theirs [13]. Secondly, we identified a cluster characterized by obesity with a high symptom burden, similar to the PRISM cohort. This cluster shows a decrease in FEV1% predicted compared to cluster 1, but due to a synchronous decrease in FVC% predicted, it has the highest FEV1/FVC ratio among the 5 clusters. Our results partially confirm PRISM’s features, such as higher BMI and higher proportion of females [8, 26]. Moreover, this cluster is primarily composed of individuals without COPD, with 64% of the patients being classified as GOLD stage 0. Thirdly, subtypes related to severe emphysema have been proposed in previous studies, and our study proposes a severe emphysema subtype with high airflow obstruction and high airway wall thickness [13, 14].

Additionally, this study confirms that the disease axis exhibits relatively high reproducibility in different studies. Compared with the factor loading done by KINNEY et al., where their factor 1 predominantly represented quantitative CT measures of air trapping and emphysema, and factor 2 predominantly represented lung function test measures, our study exhibits the similar results [19]. In addition, rotated principal component 3 (RC3) and rotated principal component 4 (RC4) in our study differ from factors 3 and 4 identified by Kinney et al. This discrepancy is attributed to the inclusion of a larger number of Pi10 and wall area percent indicators in our study.

This work has several limitations. Firstly, the study was based on a single cohort, and the external validity of the identified COPD subtypes needs to be confirmed in other independent datasets. Validation in different populations, regions, and clinical settings is crucial to assess whether the identified subtypes and their associated exacerbation risks are generalizable. Secondly, due to missing data, we excluded subjects with incomplete information. Although our sample size remains large and potential bias is minimal, this exclusion may still introduce some bias compared to the overall cohort.

Conclusion

In summary, through the utilization of PCA and k-means cluster analysis, we identified five clusters and their stability was validated in both the validation set and cross-validation. This will contribute to further describing and studying subtypes related to quantitative CT. In addition, this method shows high repeatability in independent data samples derived from the same cohort, providing a feasible approach for future research to address the issue of subtype repeatability obtained through unsupervised learning.