2.2. Clinical and Biological Endpoints

Progression-free survival (PFS) represents the duration from the start of treatment until either the time when disease progression is observed, the date of death due to any reason, or the last follow-up date (censored). The overall survival (OS) rate is computed from the date of the patient’s diagnosis to either the date of death from any cause or the date of censoring, whichever occurs first [28]. Imaging techniques like CT scans, PET-CT scans, and MRIs are commonly employed to confirm disease advancement [11,12].

CD8 (cluster of differentiation 8) is a transmembrane glycoprotein that functions as a co-receptor for the T-cell receptor (TCR). The CD8 co-receptor serves a role in T cell signaling and facilitates interactions between cytotoxic T cells and antigens. CD8 T cells were evaluated to determine the status of immune infiltration in tumor samples. The CaloPix image analysis software (Tribun Health) was used to score IHC slides, with the results expressed as % of positive cells’ surface per square micrometer of tumor (μm²).

2.3. Radiomics Features Analysis

To extract features from the radiological scans, it is crucial to normalize the slice thicknesses of all CT images, which was done by interpolating them to a voxel size of 1 × 1 × 1 $\mathrm{m}{\mathrm{m}}^3$ through the B-spline interpolation method available in the SimpleITK (Simple Insightful Toolkit) library. Normalizing the slice thicknesses was essential to ensure consistency and standardized representation, achieve voxel alignment, and improve the accuracy of feature calculations. These measures contribute to reliable and comparable results, enabling meaningful comparisons across different patients and imaging studies. We used the Pyradiomics version 3.0.1 computational platform in Python [29] to extract the features from the pre-treatment CT scans. Pyradiomics is an open-source software library that allows for the extraction of radiomics features from medical imaging data in any format readable by ITK compliant with the image biomarker standardization initiative (IBSI https://arxiv.org/abs/1612.07003 accessed on 20 October 2020) [30]. This platform also provides several built-in features for pre-processing, such as normalization and resampling.

Radiomics features are generally categorized into three main groups based on their properties and calculation methods: (i) morphological or shape-based features—these features pertain to the shape and structure of the image, calculated in both 2D and 3D views. This category includes a subset of attributes describing the image’s shape, like volume, surface area, compactness, and sphericity; (ii) first-order or intensity-based features—these features quantify the intensity characteristics of the tumor based on the statistical properties of tumor intensity distribution within the region of interest (ROI). They include metrics like mean, median, and mode, along with indicators of distribution symmetry and diversity such as percentiles, skewness, kurtosis, and entropy; and (iii) second-order or texture-based features—these features consider the statistical dependencies and spatial relationships between adjacent voxels. Several methods are used to figure them out, including the gray-level co-occurrence matrix (GLCM), the gray-level run-length matrix (GLRLM), the gray-level size zone matrix (GLSZM), the gray-level distance zone matrix (GLDZM), and the neighborhood gray-tone difference matrix (NGTDM). For instance, GLCM quantifies the occurrence of pairs of pixel intensities at different spatial distances and orientations within the ROI, presenting them in matrix form. In addition, filter-based (e.g., wavelet, Laplacian of Gaussian, logarithm, exponential) radiomics features can be extracted from derived images using filters that can be applied to both texture- and intensity-based features.

2.4. Whole Slide Image (WSI) Pre-Processing

Pre-processing of WSIs is a crucial preliminary procedure before patch extraction, and various anomalies may occur during the process of digitizing WSIs. These anomalies observed in WSIs encompass various artifacts, such as noise and background interference. Prior to patch extraction, a thorough pre-processing of the WSIs is essential to methodically identify regions of interest. The ‘histolab’ module in Python was employed to carry out the WSI pre-processing methods. In order to clearly define the tissue region while eliminating artifacts and background interference, a series of procedures were implemented. These procedures consist of the following steps: conversion of the image to grayscale, application of Otsu thresholding, binary dilation, elimination of small holes, removal of small objects, etc. By using the aforementioned procedure, we obtained WSIs (without noise and artifacts) that exhibit a minimum of 80% tissue content. This procedure ensured that only the relevant part of the WSIs, known as the region of interest, was preserved for the downstream analyses. As the goal of this study was to investigate the relationship between radiomics features, pathomics features, and clinical endpoints, it was critical to consider only those areas of the WSI that have a significant number of nuclei.

This focus on regions with a substantial number of nuclei is driven by the need to ensure that the extracted patches are representative of the tissue’s cellular structure and pathology. Therefore, the pre-processing steps were designed to enhance the visibility and contrast of the nuclei within the tissue, facilitating their accurate identification and analysis. Moreover, this meticulous pre-processing approach aids in minimizing the potential biases that could arise from non-tissue elements or artifacts present in the WSIs. By concentrating on areas rich in nuclei, the study aims to establish a more robust and direct correlation between the extracted features and the clinical endpoints, thereby providing more reliable and clinically relevant insights. This careful selection of regions for analysis is not only crucial for the accuracy of the study but also enhances the reliability of the conclusions drawn regarding the interplay of radiomics and pathomics features in the context of clinical outcomes.

2.5. WSI Segmentation and Patching

The extracted patches from the WSI measure 512 × 512 pixels with a 20× magnification. All patches of this size containing a nucleus percentage of 15% or higher were extracted [31]. The extraction of these patches was accomplished using the ScoreTiler class from histolab (available online: https://histolab.readthedocs.io/en/latest/readme.html, accessed on 20 October 2020). The ScoreTiler class computes a score for each patch extracted from the WSI based on the detected nuclei count. It incorporates multiple techniques like threshold-based techniques, color space conversion [32,33], watershed transformation [34,35], and morphological operations [36,37] to facilitate nucleus detection. Patches were subsequently divided into bins of different scales (16 × 16 px, 32 × 32 px, and 64 × 64 px). The estimation of cell density, which measures the spatial arrangement of cells within a tissue, is conducted for individual patches at a specific scale, as previously described by Alvarez-Jimenez et al. [4]. The cell density was determined by the segmentation of nuclei and the allocation of a certain gray level to each patch, which was determined by the estimated number of nuclei contained within that region [4,5]. The resulting output consists of a map that displays the density of cells at different scales, providing a visual representation of the distribution of cell density across different areas within the patch. The approach utilized for nucleus segmentation in this study was the Watershed algorithm [35,38]. Prior to using the algorithm, several preprocessing stages were conducted on the image patches. We removed small white noises present in the image, as well as small holes in the object (nuclei), to ensure that the region near the center of the object was considered as the foreground region, while those farther away from the object were regarded as the background. Subsequently, the area that represents the nuclei was identified by using the distance transform with an appropriate threshold. We then isolated the area that was not the nuclei by applying a dilation function, which shifts the object’s boundary towards the background [39]. Following the isolation process, the Watershed algorithm was employed. This algorithm functions by discerning the high and low points in the grayscale image, thus effectively outlining each nucleus, even within regions where they are densely packed. Such segmentation is pivotal for accurately estimating cell density, as it facilitates the precise identification and enumeration of individual nuclei in each patch. Post-segmentation, the nuclei were categorized according to their size and shape attributes, enhancing the analysis of cell density. This categorization is crucial for differentiating among various cell types and conditions, which is essential for a deeper comprehension of the biological activities in the tissue. Employing this combination of techniques allows for an exhaustive and nuanced examination of the tissue microenvironment, yielding critical insights beneficial for both diagnostic and research objectives.

2.6. Extraction of Pathomics Features

Haralick features constitute a collection of statistical descriptors employed to characterize the texture or patterns discernible within an image. These features were derived from the gray-level co-occurrence matrix (GLCM), which is a matrix recording the frequency of occurrences of distinct pairs of pixel intensity values at predefined relative positions within an image. The utilization of Haralick features in the analysis of cell density maps enables a quantitative depiction of the spatial distribution and patterns of cell density. These features were obtained from each patch for each scale. From each matrix of cell density values, four co-occurrence matrices were computed for each scale (horizontal, vertical, minor diagonal, and main diagonal).

2.7. Radiopathomics Analysis

To investigate a possible association between radiomics and pathomics features in immunotherapy-treated NSCLC patients, an integrative study design was developed as presented in Figure 1. The initial row delineates the steps carried out on the pre-treatment CT scans for radiomics analysis, while the subsequent row outlines the steps used for pathomics analysis. After performing radiomics and pathomics feature extraction, a comprehensive analysis was conducted to explore the relationship between radiomics and pathomics, radiomics and clinical outcomes (overall survival, progression-free survival, and CD8), and pathomics and clinical outcomes (overall survival, progression-free survival, and CD8). The association between imaging-derived features was computed using the Spearman correlation. In order to evaluate the variations in correlation distributions between radiomics and pathomics features and clinical outcomes, the Kolmogorov-Smirnov (KS) test was utilized. This statistical measure, obtained from the scipy.stats package in Python, is employed to quantify the similarity between two distributions.

Open in a separate window

Figure 1

Process flow of the radiopathomics evaluation employed in this research: The upper row illustrates the sequential radiomics analysis process, from initial radiological imaging to feature extraction and statistical analysis. The lower row details the corresponding pathomics analysis, starting from histological slide preparation to intricate feature extraction, with both pathways converging on the assessment of clinical endpoints like OS (Overall Survival), PFS (Progression-Free Survival), and CD8 cell count distributions.

2.8. Feature Reduction Analysis

Initially, data normalization was done by subtracting the mean and dividing by the standard deviation. This type of data normalization, referred to as z-score normalization, is vital as it ensures that all data sets operate on the same scale, which is crucial for correlation studies. A correlation filter was then constructed to reduce feature redundancy. This correlation filter is based on Spearman’s correlation coefficient, which is known for its ability to detect monotonic relationships and its robustness against outliers [4]. A threshold of |r| = 0.9 was established, and all pairs of correlations between pathomics and radiomics features with an absolute correlation value greater than 0.9 were examined individually [5,40]. In summary, when two features exhibit the absolute value of the correlation coefficient (r) exceeding 0.9, the algorithm assesses the mean absolute correlation of each variable, and the variable with the highest mean absolute correlation is subsequently eliminated. Spearman’s correlation was also used to compute the correlation between radiomics features and clinical endpoints (OS, PFS, and CD8), as well as between pathomics features and clinical endpoints (OS, PFS, and CD8). To control for the false positives, we performed multiple hypothesis testing using the False Discovery Rate (FDR) method.

3.2. Feature Extraction and Reduction

Using the Pyradiomics pipeline, a total of 851 radiomics features were extracted from the annotated ROI of the pre-treatment CT scans. For the pathomics part, a total of 13 haralick features were computed for each co-occurrence matrix, resulting in a cumulative count of 156 features per patch. In order to obtain the features at the WSI level, five statistical measures (mean, median, variance, kurtosis, and skewness) were computed over all patches inside each WSI. Consequently, WSI is characterized by a total of 260 distinct attributes for each scale.

The removal of highly correlated variables will enhance the precision, stability, and interpretability of the correlation analysis between radiomics and pathomics features, which was crucial for obtaining reliable and informative results. Following the application of the correlation filter to remove redundant features, the number of pathomics features decreased from 260 to a total of 88, distributed as 17 features at the 16 × 16 px scale, 30 features at the 32 × 32 px scale, and 41 features at the 64 × 64 px scale. For the radiomics features, after applying the correlation filter, only 156 out of 851 extracted features were retained, with the other 695 being excluded due to their strong mutual correlation (|r| > 0.9). These sets of 88 and 156 pathomics and radiomics features were utilized for downstream analyses.

3.3. Association between Pathomics and Radiomics Features

Figure 2A displays the correlation heatmap between pathomics and radiomics features, while Figure 2B illustrates the comparison of the percentage positive correlations and negative correlations. For the scale of 16 × 16 px, the range of correlation coefficient values was found to be between −0.60 and 0.58. Positive correlations account for 49% (i.e., 1314) and negative ones for 51% (i.e., 1355) of the total correlations. The most strongly correlated radiomics features were the wavelet features. On the pathomics side, features such as “sum variance, sum average, and Information Measure of Correlation 2” were predominant. For the scale of 32 × 32 px, correlation values ranged from −0.57 to 0.57. Positive correlations account for 47% (i.e., 2195) of the total, whereas negative correlations account for 53% (i.e., 2514). The wavelet features from the radiomics and the Information Measure of Correlation 1 and 2, sum average, and difference variance from the pathomics are the features that were most correlated to each other. For the 64 × 64 px scale, correlations between the features extracted from the cellular density map and the radiomics features range between −0.64 and 0.58. Negative correlations account for 54% of the total (3481), while positive correlations account for 2950. Wavelet features among the radiomics features were the most strongly associated with pathomics. To summarize our findings, we found significant cross-scale associations between radiomics- and pathomics-derived handcrafted features.

Open in a separate window

Figure 2

(A) Detailed correlation matrix heatmap. (A) Displaying Spearman’s Rank Correlation Coefficients between pathomics and radiomics features across multiple scales: The color gradient from red to blue represents the spectrum of positive to negative correlations, respectively, segmented by feature category (Shape, Intensity, Texture, and Filter) and spatial resolution (Scale 16 × 16, Scale 32 × 32, Scale 64 × 64). The intensity of the color indicates the strength of the correlation, providing insights into the multidimensional relationship between imaging characteristics and histopathological patterns. (B) Bar Chart illustrating the proportion of positive and negative correlations between radiomics and pathomics features at varied analytical scales: This chart quantitatively contrasts the prevalence of positive (blue) and negative (green) correlation coefficients, providing a clear visual representation of the directional relationships across distinct spatial resolutions in the dataset.

The significant radiopathomics association given by the highest correlation coefficient (r > 0.5) is presented in Table 2.

We also investigated the association of pathomics features with different groups of radiomics-based features. Figure 3 illustrates the heatmap of correlations between various groups of radiomics features and the common correlated pathomics features across all scales. In the variable naming, “S1” refers to the 16 × 16 scale, “S2” to the 32 × 32 scale, and “S3” also refers to the 64 × 64 scale. Filter-based features, as well as pathomics features like correlation, information measure of correlation 1, and sum average, were found to be significant in this association study. These findings suggest that there are meaningful relationships between textural and structural characteristics extracted from radiological and pathological data. The presence of correlations between these different types of features highlights the potential for synergistic information when combining radiomics and pathomics analyses, which can lead to a more comprehensive understanding of the underlying biological and clinical aspects of the studied NSCLC cases.

Open in a separate window

Figure 3

Detailed heatmap visualization of correlation patterns between radiomics categories and common pathomics features across all analytical scales: this comprehensive figure delineates the extent of correlation, with shades of red and blue representing the degree of positive and negative correlation, respectively. Each row corresponds to a specific feature, categorized by shape, intensity, texture, or filter, and each column represents the correlation of that feature with common pathomics features at different scales, offering a nuanced view of the interplay between radiological imaging attributes and pathological characteristics.

3.4. Association between Imaging Features and Survival Endpoints

The association between imaging features, namely, radiological and pathomics with the two survival endpoints, PFS and OS, was examined. The range of the absolute value of the correlation coefficient between pathomics features and OS was found to be between −0.34 and 0.34, while for PFS, the range was between −0.32 and 0.33. Conversely, the correlation coefficients between radiomics features and OS ranged from −0.26 to 0.37, and for PFS, they ranged from −0.23 to 0.27 (Figure 4). Within the pathomics feature space, we found 54 and 34 features to be positively and negatively correlated with PFS, respectively, while 47 and 41 features were found to be positively and negatively correlated with OS, respectively. Similarly, within the radiomics latent space, we obtained 82 and 74 to be positively and negatively correlated with PFS, respectively, and 83 and 73 features to be positively and negatively correlated with the OS, respectively.

Open in a separate window

Figure 4

Correlation distribution between pathomics and radiomics features and survival endpoints (left—OS and right—PFS).

Furthermore, the KS test was employed to examine the similarity between the distribution of correlation values produced by radiomics and pathomics features and the survival endpoints. The KS test yielded a statistic of 0.202 with a p-value of 0.017 for comparing the distribution of radiomics and pathomics with PFS. The analysis revealed a significant difference in the distribution of correlation coefficients between radiomics and pathomics features concerning PFS. In contrast to PFS, the analysis of correlations with OS showed no strong evidence of a significant difference in the distributions of correlation coefficients. The KS statistic for OS was 0.137, with a p-value of 0.217. These results indicated that the distribution of the correlations between radiomics and pathomics features and OS is similar.

3.5. Association between Imaging Features and CD8 Counts

A study on the association between imaging features and CD8 endpoints was conducted. The analysis revealed that the correlation between pathomics features and CD8 ranges from −0.35 to 0.33, and in the case of radiomics features, it spans from −0.34 to 0.41, as illustrated in Figure 5. The number of positive and negative correlations between pathomics features and CD8 was the same, with 44 for each. In the case of radiomics features, 106 positive correlations were observed against 50 negative ones. To compare the two distributions, we performed the KS test. We obtained the KS statistic of 0.196 with a p-value of 0.022, suggesting that there is a statistically significant difference between the two distributions and that the two sets of imaging features exhibit a distinct association with the CD8 immune phenotype. The main variables influencing the correlations between pathomics features and clinical endpoints were identified as contrast, correlation, sum variance, sum average, and Information Measure of Correlation 1. However, the radiomics features most influencing correlations with clinical endpoints predominantly belonged to the filter-based features.

Open in a separate window

Figure 5

Correlation distribution between pathomics and radiomics features and CD8.

We would like to thank the Imagia-Canexia-Health company for providing the annotated images, which was done through their Evidens platform. We would like to thank the Quebec Heart & Lung Institute Research Center biobank for providing the specimens. All clinical data collected as well as wet laboratory experiments were made possible through the Oncotech grant supported by MEDTEQ+, TransMedtech Institute, Oncopole, and the Fonds de Recherche Québec-Santé.