The RNA-Seq upper quartile normalized RSEM data was available for 539 ccRCC; all data is accessible via the NCI genome data commons and the Gene Expression Omnibus (https://gdc.cancer.gov/ and https://www.ncbi.nlm.nih.gov/geo/). GSE150404 has 60 samples with different stages. Clinicopathologic data for the corresponding patients, including gender, race, age, tumor location, histology classification, differentiation grade, tumor stage, and survival information, were also retrieved from the database. Only patients with both survival information and expression data available were included in this study. Analysis of the RNA data was split into miRNA analysis, mRNA signature analysis, and immune gene signature analysis. All subsequent analyses were performed in R open-source programming language. For differential expression analysis, RPKM values were calculated from RNaseq raw counts and upper quantile normalized. For hierarchical clustering and WGCNA, raw count data were processed and normalized using the variance stabilizing transformation algorithm implemented by the DESeq2 package.

Somatic exome sequencing data were downloaded for the 255 ccRCC. The tumors with sequencing data are accessible via the NCI genome data commons (https://gdc.cancer.gov/). The SMG that had been previously identified by the MutSigCV algorithm in the previous TCGA KIRC publications were used as the reference SMG.

All data is available from TCGA data. To minimize the influence of tumor purity, we dichotomized the methylation data with a beta value cut off 0.3 and used the Ward method to cluster the distance matrix computed with the Jaccard Index. Heatmaps were generated based on row and column orders calculated as above and colored by dichotomized beta values. All methylation genes we selected with a beta value of 0.3 or more were considered evidence for epigenetic silencing and associated with poor survival.

To find EMT-specific genes in ccRCC, gene expression data were analyzed from TCGA cohort. The EMT-related gene expression signature comprised 200 genes obtained from gene set “hallmark epithelial mesenchymal transition” in The Molecular Signatures Database (MSigDB, software. broadinstitute.org/gsea/msigdb). Genes were selected when the mRNA expression levels were either positively or negatively correlated with at least one of the well-known EMT markers: E-cadherin (CDH1), vimentin (VIM), N-cadherin (CDH2), and fibronectin 1 (FN1). Using a gene feature and its correlated genes, hierarchical clustering analysis was performed with the centered correlation coefficient as the measure of similarity ( Supplementary Table S1 ). The patient clustering result divided the patients into two subtypes (mesenchymal and epithelial). We also calculated EMT scores based on previous studies and obtained comparable results.

All the slides were tiled in non-overlapping 256×256-pixel windows at a magnification of 10× using the openslide library. The tiles with more than 50% background (for which all the values are below 220 in the RGB color space) were considered to have insufficient information and were dropped while tiling (11). For the remaining tiles, we yielded 1150952 tiles, with an average of 3326 tiles per slide (minimum 174 tiles, maximum 9500 tiles). We used the Macenko method for color normalization to convert all images to a reference color space for all these tiles (12).

For tumor detection in ccRCC, we selected 1000 tiles of necrosis, stroma and tumor from the tiles respectively, and randomly split them into training, validation and testing set at a ratio of 7:1.5:1.5. We trained a convolutional neural network with the Inception V3 model to classify necrosis, stroma, and tumor tiles by transfer learning. This network was then used to detect all the tumor tiles generated from the whole-slide images. This process resulted in 1031666 tiles, with an average of 2982 tiles per slide (minimum 170 tiles, maximum 7930 tiles).

For EMT subtype detection, we trained another Inception v3 model using just the tumor tiles detected before. As we already obtained the EMT subtype for patients’ slides using our gene signature, we labeled the tiles with Epi or Mes according to the slides’ EMT subtype. That is, the label of tiles was the EMT subtype of the slides the tiles generated from. Due to the small number of Epi slides, training and validation were performed slide-wise using six-fold cross-validation. That is, the tiles associated with the same slide were grouped into the same fold. Because there were different numbers of Epi and Mes slides (6 Epi and 52 Mes in each fold) and different numbers of tumor tiles of each slide, there was an imbalance between these two classes. To minimize such imbalance, while training, we extracted all the Epi slides and downsampled the Mes slides at a ratio of 1:2 (Epi : Mes), and then randomly selected 600 tiles per Epi slide and 300 tiles per Mes slide. However, for slide-wise validation, all the tumor tiles were used. After getting the classification of all tumor tiles of a given slide, the mean prediction value of all these tiles was regarded as the slide’s final prediction result.

Both the tumor detection and EMT subtype detection models were based on Inception V3 and trained using transfer learning. All the convolutional neural networks were pretrained on the ImageNet (www.image-net.org) database. Only the last softmax layer was changed and retrained by our dataset. We used Stochastic gradient descent (SGD) optimization with a learning rate of 0.01 for the classification model. All codes were implemented in Python 3.7 and ran on desktop workstations with Nvidia graphics-processing units (GPUs; NVIDIA 2080). Performance was scored using accuracy, AUC, specificity, and sensitivity.

Immune gene signatures were derived from previously published work (13). RSEM upper quartile normalized, log2 transformed, and mean centered RNA-seq data were matched to predefined immune gene signature clusters via Entrez IDs. Each gene signature was included in Supplementary Table S2 . Differential expression for each gene signature was analyzed between subtypes via one-way ANOVA. These P values were adjusted for multiple tests using the Benjamini-Hochberg procedure. For hazard ratio forest plots, a univariate Cox proportional hazards (CoxPH) model was used with signature/clinical variable as a continuous variable compared to patient overall survival.

CIBERSORT is an analytical tool that accurately quantifies the relative levels of distinct immune cell types within a complex gene expression mixture to characterize and quantify each immune cell subtype. Here, we applied the original CIBERSORT gene signature file, LM22, which defined 24 immune cell subtypes, and analyzed datasets from ccRCC. The data were normalized using the cubic spline algorithm. All samples were analyzed for immune cell profiles by CIBERSORT, the number of permutations being set to 100. Twenty-two immune cell types together with CIBERSORT metrics as Pearson correlation coefficient, CIBERSORT P value and root mean squared error (RMSE) were quantified for each sample. CIBERSORT P value reflects the statistical significance of the deconvolution results across all cell subsets and was useful for filtering out deconvolution with less significant fitting accuracy (https://cibersort.stanford.edu). Immune cell profile was calculated for each sample, and mean values were calculated.

The ESTIMATE algorithm was applied to the normalized expression matrix for estimating the stromal and immune scores for each clear cell renal cell carcinoma sample. Access to the deidentified linked dataset was obtained from TCGA in accordance with the database policy. For analyses of deidentified data from TCGA database, Institutional Review Board approval and informed consent were not required.

Student T test, Wilcoxon test, and Kruskal–Wallis test were utilized to compare continuous variables and ordered categorical variables, such as mutation load, neoantigen load, HRD score,

CTAs number, and ITH. Permutation test was conducted in the comparison of gene mutation frequencies among clusters. Correlation matrices were created with Pearson or Spearman correlation. Survival analysis was performed using the Kaplan–Meier method, and the survival of the clusters was compared using the log-rank test. For all analyses, significance was determined as a P value < 0.05 and corrected for multiple testing where specified. Univariate analysis was performed unless otherwise specified. Survival analyses were performed using GraphPad Prisma (GraphPad Software, Inc.) or by individually specified methodologies.

Data are available from the authors upon request but may require data transfer agreements. No personalized health information will be shared.

The code used during the current study are available from the corresponding author on reasonable request.

An EMT-related gene signature was developed and consisted of 63 genes ( Supplementary Table S1 ) and the flowchart of strategy to identify EMT-related biomarkers in ccRCC was shown in Figure 1A . Two distinct subtypes of ccRCC, mesenchymal (Mes) and epithelial (Epi), were identified using hierarchical clustering of gene expression data. The Kaplan-Meier plots showed that the patients in the Mes subtype showed worse overall survival than those in the Epi and Intermediate subtypes (P = 0.009; Figure 1B ).

The independent cohort (GSE150404) was used to validate the accuracy of the EMT gene signature (P < 0.0001; Supplementary Figure S1 ). Furthermore, higher grade (grade 3/grade 4) and stage (stage III/stage IV) samples were in the Mes group ( Figures 1C, D ). Univariate and multivariate analyses were performed and EMT gene signature was significantly associated with outcome in the multivariate analysis ( Table 1 ).

According to TCGA Research Network article, unsupervised clustering methods identified subsets in mRNA (m1-m4) and miRNA (mi1-mi4) expression datasets. In the mRNA cluster, m3 accounted for roughly 50% of samples in the Mes classification and was associated with the worst survival outcome ( Figure 1E ). Survival differences were also evident in miRNA-based subtypes ( Figure 1F ). DNA methylation (me1-me3) expression datasets were detected ( Figure 1G ), and almost one-half of the Mes group was me1 type and associated with significantly poorer survival.

Thorsson et al. identified six immune signature sets (C1-C6) associated with overall survival and progression-free interval (14). C3 had the best prognosis and was enriched in most ccRCC, the more mixed-signature subtypes, C4 and C6, had the least favorable outcome. In our model, compared with Epi, Mes had the most significant proportion of C4+C6 ( Figure 1H ).

Further, we quantified the EMT levels by calculating EMT scores described by Tan et al (15). Positive EMT scores corresponded to the mesenchymal phenotype, whereas negative scores reflected the epithelial phenotype. In general, All Mes subtype was characterized by positive EMT scores corresponding to their phenotype. In contrast, Epi phenotype had intermediate and low EMT scores ( Figure 1I ), reflecting their partial EMT and complete epithelial states. These results demonstrated the robustness of the prognostic value of EMT signature.

Recently, deep learning has overmatched humans in some medical data processing, especially the ability to predict some genetic information efficiently and economically using histopathology images (16, 17). To facilitate applying the EMT gene signature, we investigated whether deep learning could directly predict EMT status which according to our molecular subtype from H&E-stained slides.

In this study, we used 346 whole-slide images from TCGA and labeled them with the EMT types which according to our gene signature ( Figure 2 ). Before classifying the EMT status, 3000 tiles of tumor, necrosis, and stroma were selected and labeled by a urologist and a pathologist. We trained an Inception V3 model to recognize tumor from the other two classes by transfer learning, which yielded an average AUC of 0.99 at the tile level in five independent experiments. Then we used the tumor tiles detected by the model with best performance among the five experiments for EMT subtype classification. For the slide-wise six-fold validation of the EMT subtype, the AUC value of the ROC curve for all classes was 0.84 ( ± 0.07). The mean validation accuracy at slide level was 74.90%, mean specificity was 72.23%, and mean sensitivity was 75.32% ( Figures 3A – C ). The predicted EMT status and the classification probability, were aggregated to extract the heatmap for visualization by the tiles’ raw position of the slide ( Figure 3D ).

To detect which histopathological features were the most relevant in our algorithm to identify molecular subtypes, we used class activation mapping (CAM) (18). For the Epi subtype, these regions had a looser arrangement, big cell gap, nucleoli absent or inconspicuous which were mainly in pink, and granular eosinophilic cytoplasm areas ( Figure 4A ). Several studies suggested that granular eosinophilic cytoplasm reflects the presence of abundant mitochondria, which are required to supply energy for acid secretion (19, 20). For the Mes subtype, the regions were densely packed, surrounded by arborizing vasculature, the large multinucleate cells with empty cytoplasm (different in shape and size), and were often surrounded by abundant immune infiltration ( Figure 4B ).

We further investigated the molecular alterations that could be correlated with the histopathological features. Notably, the genes located in the Mes subtype were annotated to various immune associated pathways and biologic processes ( Supplementary Figure S2 ). In addition, we noticed that Mes group exhibited higher artificial intelligence (AI) score (> 0.5 Mes possibility; < 0.5 Epi possibility) and was positively correlated with immune infiltration (Spearman correlation = 0.275, P < 0.001, Figure 4B ). Moreover, gene set enrichment analysis (GSEA) validated that abnormality of the mitochondrion and varieties of metabolic process was enriched in the Epi cluster, which was in accordance with the detection above ( Figure 4A ).

Our tests showed promising results on tumor EMT classification from sections that could be predicted from H&E images. Furthermore, we analyzed the connection between DNA and RNA-level alternations and histopathological features to explore the clinical relevance of the classification.

Heatmapping demonstrated clustering of 89 highly variable DNA methylation markers that were hypermethylated in the tumor. A β-value R greater than 0.3 was considered hypermethylated, and the markers were associated with significantly poorer survival (13). The methylation Mes subtype had a pronounced increase in hypermethylation across the selected genome ( Figure 5A ).

Evaluation of metabolic differences was performed by 15 major metabolic processes. Expression of the Krebs cycle and the electron transport chain genes (complex I – complex IV) provided a clear distinction between the subtypes ( Figure 5B ).

AMP-activated protein kinase (AMPK) acts as an intracellular energy sensor and was significantly lower in Mes subtype compared to Epi subtype (P < 0.01; Figure 5B ), which negatively regulates fatty acid synthesis and positively regulates mitochondria production (21). The metabolic shift identified in the Mes group is reported to contribute to the Warburg metabolic phenotype, further enhanced malignancy, immune protection of cancer cells (22). Furthermore, the GSEA and AMPK complex genes expression validated some of these results ( Figures 5C, D ).

The pattern of somatic alterations was determined from the analysis of 104 samples (24 Epi and 80 Mes). Based on previous study methodologies, we identified the eight most significantly mutated genes (P < 0.00001) (23). As shown in Figures 6A, B , patients with PBRM1 mutation tended to be in the Mes subtype, while the Epi subtype had higher mTOR mutant rates (P < 0.05). PBRM1 encodes the bromodomain-containing protein BAF180, a subtype of the switch/sucrose non-fermentable (SWI/SNF) chromatin remodeling complex (24). Importantly, retrospective data have shown that patients with metastatic ccRCC harboring truncating mutations in PBRM1 experienced increased benefit from immune checkpoint therapy (25, 26). Moreover, the Epi subtype tended to have a better response to the mTOR inhibitor.

To further investigate the relationship between the EMT-associated subtype and the immune, we first established a microenvironment phenotype calculated by ssGSEA ( Supplementary Table S2 ). Mes subtype was characterized by high innate and adaptive immune cells infiltration ( Figure 6C ). The ESTIMATE algorithm showed that the Mes subtype had the highest immune and stromal scores for the analyzed cohort (P < 0.01; Figures 6D – G ).

To identify whether different clusters of EMT had distinct tumor extrinsic immune escape mechanisms (27), we focused on the adaptive immune cells in ccRCC ( Figure 6H ). Analysis revealed that the fraction of naïve B cells was higher in Mes subtype tissue than in Epi tissue. The three main T cell subpopulations in tissue (CD4+ memory resting T cells, CD8+ T cells, and follicular helper T cells) were increased in Mes compared to Epi. The results also showed that innate immune components, neutrophils, M2 macrophages and resting NK fractions increased. In general, Mes had abundant active innate and adaptive immune cells and immunosuppressive cells.

Consistent with these results ( Figure 6I ), Mes had a higher expression of chemokines, including CCL4, CXCL9, and CXCL10, which have been proven to attract dendritic cells and CD8 T cells (28). Overall, the increase of chemokines might contribute to the extrinsic immune escape of the Mes clusters.

We further investigated the two aspects of intrinsic immune escape: tumor immunogenicity and immune checkpoint molecule expression (29). We first compared the factors that were the primary source of tumor antigens: mutation load, neoantigen load, CTA (cancer-testis antigen) level, HRD (homologous recombination deficiency), ITH (intratumoral heterogeneity), and tumor antigen-presenting capability. In general, the difference in the tumor antigen burden among the clusters was significant between the two groups ( Figures 7A – E ). Considering that high tumor mutation burden (TMB) correlates with a greater probability of displaying tumor neoantigens, it is rational to hypothesize that the tumors with Mes subtype are more likely to respond to immune checkpoint inhibitors (ICIs) as this greater mutation load.

Therefore, we referred to a database of immune checkpoint molecules among clusters (30) ( Figure 7F ). In many cancers, the PD-1 pathway is involved in tumor evasion from immune activity (31). PD-1 (PDCD1) and CTLA-4 (CTLA4) levels were significantly elevated in the Mes subtype as compared with those in the Epi subtype (p < 0.01); Furthermore, we demonstrated that PD-1, CTLA4, LAG3, TIGIT, CD80, and CD86 and the expression of most checkpoint molecules were positively correlated ( Figure 7G ). In addition, analysis of overall survival of CTLA4 was positively correlated with poor patient survival ( Supplementary Figure S3 ).

To further test the clinical relevance, a heatmap of genes previously defined and representing angiogenesis and immune biology was investigated ( Figure 7H ) ( 32). Our data suggested the association of Mes subtype with improved response to antiangiogenic drugs. We further conducted confirmatory studies to evaluate molecular subtypes by deep learning method in 34 stage IV patients enrolled in the study from 2 January 2020 to 2 June 2021. Patients received more than one dose of pazopanib. In our cohort, only four patients in the pazopanib arm were identified as Epi term, the average score of the remaining patients was 0.72 (95% CI, 0.54–0.83) for Mes tiles and 0.28 (95% CI, 0.16–0.40) for Epi tiles ( Supplementary Table S3 ). In Mes-type patients, the objective response rates were 47% (10% complete response; 37% partial response); no Epi-type patients benefitted from the pazopanib. Of note, we observed a trend for improved efficacy (objective response rates) with Mes subtype, indicating that our approach can be beneficial to targeted cancer therapy.