Integrative Analysis of Histopathological Images and Genomic Data Predicts Clear Cell Renal Cell Carcinoma Prognosis

Histopathologic images confer important information for diagnosis, staging, and prognosis for cancers and are being used extensively by pathologists in clinical practice. With the recent availability of digital whole-slide images (1), automated computational histopathologic image analysis systems have shown great promise in diagnosis and the discovery of new biomarkers for cancers such as breast (2–4), lung (5, 6), brain (7), and colon cancers (8). In comparison with human inspection, computerized image analysis has great potential to improve efficiency, accuracy, and consistency. Besides histopathologic images, molecular characteristics, such as genetic alterations and gene expression signatures, are also widely adopted for predicting clinical outcomes for cancers (9, 10). Therefore, an interesting scientific question is the relationship between morphologic and genomic features while an important translational question is if the integration of these two types of features can lead to more accurate prediction of patient outcome. This has been previously explored in various cancers including breast, ovarian, and glioblastoma, and led to new insights into the relationship between cancer tissue morphology and genetic changes such as PTEN mutations (3, 11–13).

To study these issues, matched histopathologic images and genomic datasets for cancers are needed. Fortunately, The Cancer Genome Atlas (TCGA) project not only provides an extensive collection of genomics and clinical outcome data for large cohorts of patients of more than 30 types of cancers, but also hosts a large collection of matched histopathologic images for solid tumor samples. Currently, more than 24,000 histopathologic images are available at the TCGA data portal and can be visualized at the Cancer Digital Slide Archive (CDSA, http://cancer.digitalslidearchive.net/; ref. 14).

Quantitative analysis of these images and integration with genomics data require innovation in integrative genomics and call for techniques from bioimage informatics, genomics, and bioinformatics. We previously developed a computational framework for quantifying morphological features from large histopathologic images (4, 15) and genomics visualization tools for integrating imaging, clinical, and genomic features to predict patient outcomes (4, 16, 17). Therefore, to further promote this emerging integrative genomics field straddling bioimage informatics and genomics and ensure wide utilization of valuable large datasets, we demonstrate an integrative genomics workflow on the less well-studied renal cancers. The analysis tools are publicly available and can be adopted as building blocks for other integrative genomics workflows (see Materials and Methods section).

Renal cell carcinoma (RCC) is the most common type of malignant neoplasm arising from kidney in adults, responsible for approximately 90% to 95% of all cases (18). It can be categorized into the following histologic subtypes: clear cell, papillary, chromophobe, collecting duct, and unclassified RCC based on the Heidelberg classification system (19). In this study, we focus on clear cell RCC (ccRCC), which is the most prevalent subtype, accounting for 80% to 90% of all RCCs (20). In clinical practice, tissue sections are examined under a microscope by pathologists to make a diagnosis and predict prognosis. The clinical behavior of ccRCC is quite diverse, ranging from slow-growing localized tumors to aggressive metastatic disease (9). Therefore, prognostic markers play a crucial role in stratification of patients for personalized cancer management, which could avoid either overtreatment or undertreatment (21). For instance, patients classified into high-risk group may benefit from closer follow-up, more aggressive therapies, and advance care planning (5, 22). Currently, prognostic markers for ccRCC in routine clinical use consist mainly of tumor stage, nuclear grade, and presence of necrosis (23–25). However, cancer is a highly heterogeneous disease. The prediction accuracy of traditional clinical factors remains limited for individual patients, especially for early-stage patients, and also relies on the experience of pathologists. Therefore, there is a need for more effective markers for predicting prognosis of ccRCC.

Using the large cohort of ccRCC patients from TCGA, hundreds of cellular morphologic features can be extracted from hematoxylin and eosin (H&E)-stained whole-slide images, characterizing nucleus size, shape, texture, and the spatial relationship between nuclei. In this article, we demonstrate how image features correlate with coexpressed gene signatures and developed an automated prognostic model that could predict patient's survival risk for patient stratification, using a combination of quantitative image features and eigengenes. To the best of our knowledge, this is the first study to couple histopathologic images and genomic data to predict ccRCC clinical outcome and our results indicate that the integration of imaging and genomic features can lead to improved prognosis prediction for early-stage (stages I and II) ccRCC patients than existing clinical markers.

Processed data (extracted quantitative imaging features, combined gene expression data, etc.) and code with annotations, comments and instructions are available at https://github.com/chengjun583/image-mRNA-prognostic-model.

ccRCC patient samples used in our study included matched H&E-stained whole-slide images, transcriptome, somatic mutation, and clinical information, which were acquired from TCGA data portal at NCI Genomic Data Commons (26). Patients with missing or too short (i.e., less than 30 days) follow-up were excluded. Microscopic images (×20 and ×40 magnification) were obtained from TCGA. The demographic and clinical characteristics for the selected 410 patients are summarized in Table 1.

One challenge for this study was the lack of other large cohorts of ccRCC with matched histologic image and genomic data. Thus, instead of using a second dataset for validation, we applied cross-validation in every step of downstream of the machine learning analysis as described below.

Figure 1 outlines our data analysis workflow for both imaging and genomic data for both univariate and multivariate analyses with details of each major step being described in the following sections.

Our image feature extraction pipeline consists of three steps: nucleus segmentation, cell-level feature extraction, and aggregation of cell-level features into patient-level features (Fig. 1A). Rich pathologic information is present in stained cell nuclei that require segmentation to facilitate subsequent analyses. For this task, a recently proposed approach by Phoulady and colleagues (27) was used, which is an unsupervised segmentation method requiring no parameter learning or training data because the parameters are set adaptively. Next, 10 types of cell-level features were extracted for each segmented nucleus, characterizing nucleus size, shape, texture, and distance to neighbors. These cell-level features are nuclear area (denoted as area), lengths of the major and minor axes of cell nucleus, and the ratio of major axis length to minor axis length (major, minor, and ratio), mean pixel values of nucleus in RGB three channels respectively (rMean, gMean, and bMean), and mean, maximum, and minimum distances (distMean, distMax, and distMin) to neighboring nuclei in the Delaunay triangulation graph (28). The Delaunay triangulation graph was constructed on the basis of the locations of segmented nuclei. In this graph, each nucleus was a node and connected to neighboring nuclei. Finally, for each type of cell-level features, a 10-bin histogram and five distribution statistics (i.e., mean, SD, skewness, kurtosis, and entropy) were adopted to aggregate the numerous cell-level features extracted from a patient into patient-level features; 150 patient-level features were generated in total. Taking the cell-level feature, area, as an example, corresponding 15 patient-level features were denoted as area_bin1 to area_bin10 for the 10 histogram features, and area_mean, area_std, area_skewness, area_kurtosis, and area_entropy for the five distribution statistics. For other cell-level features, corresponding patient-level features were named in the same way. Area_bin1 represents the percentage of very small nuclei over the entire slide for a patient while area_bin10 indicates the percentage of very large nuclei in the patient sample. Skewness is a measure of the asymmetry of the data distribution around the sample mean, kurtosis is a measure of how outlier-prone a distribution is, and entropy is a statistical measure of randomness.

Additional description about aggregation of cell-level features into patient-level features is provided in the Supplementary Material. A qualitative example of nucleus segmentation results is shown in Supplementary Fig. S1.

mRNA expression profiles for the ccRCC tumors in TCGA were transformed from Illumina HiSeq 2000 RNA-seq readcounts to normalized reads per kilobase per million (RPKM). Although our first goal was to establish the relationships between gene expression data and the imaging features, the large number of genes posed a challenge to obtaining sufficient statistical power. Therefore, instead of focusing on individual genes, we first carried out gene coexpression network analysis (GCNA) to cluster genes into coexpressed modules and summarized each module as an “eigengene” using the protocol described in ref. 29 (Fig. 1B). Modules are clusters of highly interconnected/correlated genes. The eigengene of a module is defined as the first principle component, which can be considered a representative of the gene expression profiles in a module. This approach not only substantially improves statistical power (30), but also allows us to focus on important biological processes or genetic variations associated with the coexpressed gene modules, making the results more interpretable than individual genes as the coexpressed modules are often strongly associated with a specific gene group participating in the same biological process or located on the same chromosomal band.

Although there are many algorithms for performing GCNA including the well-known WGCNA package (31), we applied our recently developed weighted network mining algorithm called local maximum quasi-clique merging (lmQCM; ref. 32). Unlike WGCNA, which uses hierarchical clustering and does not allow overlap between modules, our algorithm is a greedy approach allowing genes to be shared among multiple modules, consistent with the fact the genes often participate in multiple biological processes. In addition, we have shown that lmQCM can find smaller coexpressed gene modules that are often associated with structural mutations such as copy number variation in cancers (32). The lmQCM algorithm has four parameters γ, α, t, and β. Among these parameters, γ is the most influential, as it determines if a new module can be initiated by setting the weight threshold for the first edge of the module as a subnetwork. In the lmQCM algorithm, we transformed the absolute values of the Spearman correlation coefficients between expression profiles of genes into weights using a normalization procedure adopted from spectral clustering, for which we have shown to be effective in previous studies (33). In practice, we found with γ = 0.30, t = 1, α = 1, and β = 0.4 the algorithm yielded 15 coexpressed gene modules (Supplementary Table S1) with balanced sizes and clear biological interpretation based on enrichment analysis (Supplementary Table S2).

We built a lasso-regularized Cox proportional hazards (lasso-Cox) model (R package “glmnet”) to calculate the risk index of each patient (34), based on the cellular morphologic features and eigengenes (Fig. 1C). Lasso penalty (i.e., L1 penalty) can induce sparsity and thus select an informative subset of features. To validate our method, we used a two-level cross validation (CV) strategy. After each patient was used as a test sample and classified into a low-risk or high-risk group, we used Kaplan–Meier estimator and log-rank test to test if these two groups had distinct survival. Additional description of the training and prediction process is provided in Supplementary Material.

To screen survival-associated features, for each patient-level feature we divided patients into two groups (low and high groups) where the median of each feature was used as a cut-off point. Kaplan–Meier estimator was used for patient stratification, and P value was calculated with the log-rank test, where P < 0.05 was considered significant. For the initial survival analysis, because our initial goal was screening, we did not apply multiple test compensation such as FDR control to obtain more candidate features. The lasso-Cox model was learned on the selected survival-associated features. Cox proportional hazards regression model was fitted, and 95% confidence intervals were computed to determine the prognostic values of our lasso-Cox risk indices and other known prognostic factors. Correlation was computed using Spearman rank correlation coefficients. Enrichment analysis of coexpressed gene modules was carried out using Toppgene (35). All the survival analyses were performed using R package “survival.”

To investigate which specific image features and eigengenes are associated with patient survival, we tested for each feature the statistical significance of difference in overall survival between low- and high-risk groups that were stratified by the median of feature values. Log-rank test results revealed that 33 image features and 6 eigengenes were significantly related to prognosis (P < 0.05). The log-rank test results of all survival-related variables are listed in Table 2 and the Kaplan–Meier survival curves for some variables are shown in Fig. 2A–E.

After examining these survival-associated variables, we found many of them were connected to stroma tissue. Stromal cells such as fibroblasts are typically spindle-shaped with elongated nuclei and therefore characterized by long major axes and/or large ratio between major and minor axes. As shown in Table 2, Fig. 2A and B image features such as major_bin8, major_bin9, ratio_bin8, ratio_bin9, ratio_std, and major_std were negatively related to prognosis, that is, patients with large values of these variables had worse prognosis than other patients. Large values of these variables imply a high percentage of stromal cell nuclei in whole-slide images (in terms of major_std, and ratio_std, large values of these variables mean that the major axis length and the ratio of major axis length to minor axis length are spread out in a wide range, indicating a high percentage of stromal cell nuclei). In other words, patients with high percentage of stromal tissue are related to poor prognosis for ccRCC in our study.

In addition to histopathologic images, gene expression data also corroborated that stroma played an important role in tumor prognosis. Enrichment analysis showed that gene module 2 was enriched with extracellular matrix genes (Supplementary Table S2), which is consistent with our knowledge that the tumor microenvironment plays critical roles in tumor development (2, 3). Kaplan–Meier survival curves demonstrated distinctly different outcomes for low- and high-expression groups (log-rank test P = 0.024), where high expression of eigengene 2 was associated with poor prognosis.

In the previous sections, we showed that many individual features derived from histopathologic images and genomic data stratified patients with distinct prognosis. We next investigated whether the integration of all identified survival-associated features would provide better prognostic prediction. We built a lasso-regularized Cox proportional hazards model to select the most informative features and calculate a risk index for each patient. On the basis of the risk indices, patients were divided into a low- or high-risk group by the median. The lasso–Cox model provided significantly better patient stratification than that using individual features (Fig. 2D–F, log-rank test P values = 2.23e−5, 7.46e−6, and 8.79e−10 for the most significant image feature, rMean_bin10, the most significant eigengene expression, eigengene3, and lasso–Cox model, respectively). Among the 33 survival-associated image features and six survival-associated eigengenes, eight image features and five eigengenes were selected: rMean_bin6, major_bin9, area_bin5, gMean_bin10, ratio_bin7, ratio_bin8, ratio_bin9, major_bin1, eigengene1, eigengene3, eigengene9, eigengene11, and eigengene13 (enrichment analyses of survival-related gene modules are listed in Supplementary Table S2). Both image features and eigengenes appeared in the final selected feature set, and most of the pairwise mutual information values between them are smaller than the ones between significantly correlated image features and eigengenes (Supplementary Fig. S2), suggesting that histopathologic images and genomic data complement each other in predicting survival outcome.

Genotype is one of the three factors that determine phenotype, the other two being inherited epigenetic factors and noninherited environmental factors. Therefore, tumor characteristics or morphology is very likely to have some relationships with gene expression data. To find out these relationships, we calculated Spearman rank correlation coefficients between each pair of 33 survival-associated image features and all eigengenes for the 15 modules. The heatmap of the correlation matrix is shown in Fig. 3.

As can be seen from the heatmap, eigengenes 2, 3, 9, and 11 significantly correlated with many image features (statistically significant after Bonferroni correction). The gene module 2 was enriched with extracellular matrix genes, which explained why it positively correlated with image features such as ratio_bin8, ratio_bin9, ratio_std, major_bin9, and major_std that describe the percentage of stromal cells. Gene module 3 was enriched with acid metabolic process and transmembrane transporter activity. Genes in this module play a central role in renal functions such as organic anion transport (36). Patients with low expression of this eigengene were related to poor prognostic outcome (log-rank test P = 7.46e−6; Fig. 2E), implying impaired renal function. This eigengene also negatively correlated with images features representing the amount of stromal cells such as ratio_bin9, major_bin9, ratio_std, and major_std. Gene module 9 was highly enriched with cell cycle and mitosis genes. In fact, genes in this module are frequently observed to coexpress in multiple types of cancers (37). High expression of this eigengene indicates that the tumor is more aggressive, and it was negatively related to patient prognosis (log-rank test P = 1.19e−4). Cells become bigger when they come into mitotic phase, which was in line with our observation that the gene module 9 was significantly and positively correlated with image features such as area_bin5, area_bin6, and area_std. The top molecular functions of gene module 11 by were frizzled binding and G-protein–coupled receptor (GPCR) binding. GPCRs represent the largest family of cell-surface molecules involved in signal transduction. Experimental and clinical data indicate that GPCRs have a crucial role in cancer progression and metastasis (38). Patients with high expression of gene module 11 had significantly worse outcome than other patients (log-rank P = 1.33e−3). Similar to gene module 2, module 11 also significantly correlated with many image features that describe stroma cells, such as ratio_bin8, ratio_bin9, and ratio_std. Survival analysis results and enrichment analysis results for all survival-associated eigengenes are summarized in Table 2 and Supplementary Table S2, respectively.

Using univariate and multivariate Cox proportional hazards analysis, we performed a comprehensive comparison between the lasso–Cox risk index and other known prognostic biomarkers, including two clinical variables, grade (G1+G2 vs. G3+G4), stage (I+II vs. III+IV), six gene expression signatures (39, 40), CSNK2A1, SPP1, DEFB1, PECAM1, EDNRB, TSPAN7, and five somatic mutation genes (26, 41–46), VHL, PBRM1, BAP1, SETD2, TP53. Patient subtyping for gene expression signatures was carried out by using the median as cut-off point. In terms of genes with somatic mutation, patients were classified as mutant or wild type. Of these factors, only grade, stage, lasso–Cox risk index, DEFB1, EDNRB, and TSPAN7 were associated with survival by univariate Cox proportional hazards analysis (Table 3). DEFB1 encodes β-defensin, which belongs to a family of antimicrobial peptides produced by white blood cells and epithelial cells. Rabjerg and colleagues (40) suggested that DEFB1 might be a tumor suppressor gene, but our results revealed that high expression of this gene predicted a worse prognosis with very weak significance (P = 4.99e−2; HR = 1.41, and 95% confidence interval = 1.00–1.98). EDNRB is a member of the endothelin axis, and TSPAN7 is a member of the transmembrane 4 superfamily. Wuttig and colleagues (39) showed that EDNRB and TSPAN7 might be suppressors of tumor progression and metastatic tumor growth, which is in agreement with our results that high expression of these two genes predicted a better prognosis. Subsequently, multivariate Cox proportional hazards analysis demonstrated that lasso–Cox risk index was an independent prognostic factor (P = 2.31e−4, HR = 2.26, 95% confidence interval, 1.46–3.49), as well as stage and TP53 (Table 3).

As shown in Table 3, tumor stage is the most effective prognostic factor, but its capability of stratifying early-stage (i.e., stage I and II) ccRCC patients is very limited (Fig. 4A and B). The Kaplan–Mejer curves of stages I and II are intertwined (log-rank test P = 0.962), which may be attributed to the less significant morphologic differences between stages I and II tumors and/or large subtyping variations among pathologists.

However, the image features and eigengenes can successfully stratify early-stage patients with distinct survival outcomes. Log-rank testing of each of the 165 variables (150 image features and 15 eigengenes) revealed that 13 image features and 2 eigengenes were associated with survival (Supplementary Table S3). Survival curves of three variables are shown in Fig. 4C–E. In addition, we also trained a lasso–Cox prognostic model using the above 15 variables related to survival. Figure 4F shows the survival curves stratified by the lasso–Cox risk index (log-rank test P = 0.014). Compared with individual variables, integrating image features and eigengenes did not improve the accuracy of prognostic prediction for early-stage patients while there indeed was a very significant improvement when using all patients. This is because the death rate in early-stage patients is much lower than that in all patients (18.5% vs. 32.9%), and high death rate is key to ensuring prediction accuracy of lasso–Cox model. If all patients were used in the lasso–Cox model to predict early-stage patient prognosis, the performance was improved (log-rank test P = 8.65e−3). The two eigengenes associated with the prognosis of the early-stage patients corresponded to coexpressed gene modules 3 and 13. The gene module 3 was highly enriched with genes related to kidney functions such as organic acid metabolic process (⁠|P = 5.702\ \times {10^{ - 18}}$|⁠), ammonium ion metabolic process (P = 6.612 × 10⁻⁹), and anion transport P = 5.994 × 10⁻⁸). This observation suggests that the physiologic functions for kidney can be potential prognostic markers for early-stage patients. Besides gene module 3, gene module 13 contains 10 genes. Interestingly, all the 10 genes locate on the same chromosome, straddling chromosome 14q11 to 14q32, implying potential copy number variation on 14q may be related to the prognosis of kidney patients.

Because our analysis relies on parameters for the machine learning algorithms and choices of cross validation (CV) methods, we also examined the choice of various parameters, especially the choice of number of clusters K for the cellular features. Supplementary Fig. S3 shows the log-rank test P value as a function of the number of clusters in K-means algorithm. Supplementary Fig. S3 suggests that lasso–Cox model can achieve very low P values when K ranges from 8 to 14. We also compared leave-one-out CV with k-fold CV. Supplementary Fig. S4 shows that as k increases, the P value tends to continuously decline. This is because in k-fold CV a large k means we have more training samples, and thus the learned model is likely to perform better especially when the whole data set is not very large. As a result, we chose K = 10 in the K-means algorithm, and we used leave-one-out CV in our experiments.

To our knowledge, this is the first study to predict the survival outcomes of ccRCC patients using a combination of quantitative morphologic features extracted from whole-slide tissue images and gene expression signatures. In this study, we developed an automatic image analysis pipeline to extract hundreds of cellular morphologic features, and found cellular morphology was highly linked to coexpressed gene signatures. For example, image features characterizing the amount of stromal cells positively correlated with extracellular matrix genes. SD of nuclear area correlated with genes that regulate cell cycle and mitosis. In addition, a powerful prognostic model was built to predict the survival outcomes of ccRCC patients using these two types of data. The performance of the integrated prognostic model significantly outperformed that of individual image or genomic features, which indicates that image data are complementary to genomic data for predicting patient prognosis. Using multivariate Cox regression, we verified that the risk index generated by our model was a prognostic factor independent of tumor grade, stage, and other known molecular markers. For early-stage patients, besides the imaging data, the genomic data suggest that the kidney functions and status of 14q may be predictors of the survival time for these patients.

Recent studies have underscored the important contribution of stromal gene expression and morphologic phenotypes to cancer growth and progression for breast cancer (2, 3, 47). The implication of tumor stroma to prognosis could be different for different cancer types. For instance, high percentage of tumor stroma is associated with poor prognosis in triple-negative disease but good prognosis in estrogen receptor–positive disease (48, 49). Here, we found in this study that for ccRCC both image features and gene expression signatures revealed that a large percentage of tumor stroma predicted poor prognosis.

The high resolution of whole-slide tissue images poses a great computational challenge to researchers. For this reason, many previous studies only focused on selected views in tissue microarrays or a few representative image tiles in whole-slide images (2, 5). Because tumor is a highly heterogeneous disease, image features extracted from a much larger area of the tumor would be more likely to ensure the robustness of the derived prognostic model. Our prognostic model was established on the fully automated quantitative image features that were extracted from whole-slide histopathologic images, which could avoid biases or discrepancies arising from only using a small portion of the tumor.

Our study is limited to only one large ccRCC patient cohort as it is difficult to find other cohorts that have matched histopathologic images, gene expression profiles, and survival information. However, the performance of our prognostic model was strictly assessed by cross validation. The model selection was performed by 10-fold cross validation on the training set, and then the selected model was applied to the held-out test samples to predict risk indices. Another technical contribution of this work lies in the fact that we used only the cryohistologic images from TCGA. Usually for each TCGA solid tumor sample, two histopathologic images are generated—the H&E stained diagnostic image and the cryohistologic image from a slice of tissue immediately adjacent to the tissue used for generating the omics data. Thus, due to spatial proximity, the cryohistologic image is a more accurate reflection of the molecular profiles of the tissue for the omics data. However, due to processing artifact, many of these images appear damaged and cannot be processed for tissue features using previous methods, preventing accurate characterization on the tumor morphology. Here we showed that the cell nucleic features suggestive of stromal cells indeed correlated well with the gene expression profiles of extracellular matrix and stromal genes, suggesting in such images, albeit for the artifacts affecting texture analysis, the cell nucleic features can still be used.

Finally, although our study focused on predicting survival for ccRCC patients, we believe that the workflow of integrative analysis of histopathologic images with genomic data could be easily applied to other cancer types or to predict response of specific treatments, which would allow for better patient management and cancer care.

No potential conflicts of interest were disclosed.

Conception and design: J. Cheng, X. Ye, K. Huang

Development of methodology: J. Zhang, Y. Han, X. Wang, Y. Meng, Q. Feng, K. Huang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): X. Wang, K. Huang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Cheng, J. Zhang, X. Ye, A. Parwani, Z. Han, K. Huang

Writing, review, and/or revision of the manuscript: J. Cheng, J. Zhang, Y. Han, X. Wang, X. Ye, A. Parwani, Z. Han, K. Huang

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Y. Meng, Z. Han, K. Huang

Study supervision: Z. Han, Q. Feng, K. Huang

This work was partially supported by an NCI ITCR grant 5 U01CA188547 to K. Huang, Leidos grant 15 × 014 to K. Huang, the Science and Technology Project of Guangdong Province, China (No. 2015B010131011 to Q. Feng), and Shenzhen Peacock Plan (No. KQTD2016053112051497).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.