Baodong Wang
Mei Li2- 1Department of Nephrology, Fifth Hospital of Shanxi Medical University (Shanxi Provincial People’s Hospital), Taiyuan, China
- 2Department of Laboratory Medicine, Shanxi Provincial Hospital of Integrated Traditional Chinese and Western Medicine, Taiyuan, China
Background: Identifying Kidney Renal Papillary Cell Carcinoma (KIRP) patients with high-risk, guiding individualized diagnosis and treatment of patients, and identifying effective prognostic targets are urgent problems to be solved in current research on KIRP.
Methods: In this study, data of multi omics for patients with KIRP were collected from TCGA database, including mRNAs, lncRNAs, miRNAs, data of methylation, and data of gene mutations. Data of multi-omics related to prognosis of patients with KIRP were selected for each omics level. Further, multi omics data related to prognosis were integrated into cluster analysis based on ten clustering algorithms using MOVICS package. The multi omics-based cancer subtype (MOCS) were compared on biological characteristics, immune microenvironmental cell abundance, immune checkpoint, genomic mutation, drug sensitivity using R packages, including GSVA, clusterProfiler, TIMER, CIBERSORT, CIBERSORT-ABS, quanTIseq, MCPcounter, xCell, EPIC, GISTIC, and pRRophetic algorithms.
Results: The top ten OS-related factors for KIRP patients were annotated. Patients with KIRP were divided into MOCS1, MOCS2, and MOCS3. Patients in the MOCS3 subtype were observed with shorter overall survival time than patients in the MOCS1 and MOCS2 subtypes. MOCS1 was negatively correlated with immune-related pathways, and we found global dysfunction of cancer-related pathways among the three MOCS subtypes. We evaluated the activity profiles of regulons among the three MOCSs. Most of the metabolism-related pathways were activated in MOCS2. Several immune microenvironmental cells were highly infiltrated in specific MOCS subtype. MOCS3 showed a significantly lower tumor mutation burden. The CNV occurrence frequency was higher in MOCS1. As for treatment, we found that these MOCSs were sensitive to different drugs and treatments. We also analyzed single-cell data for KIRP.
Conclusion: Based on a variety of algorithms, this study determined the risk classifier based on multi-omics data, which could guide the risk stratification and medication selection of patients with KIRP.
1 Introduction
Renal cell carcinoma (RCC) is a common malignancy of urinary system (1). Clear cell carcinoma (ccRCC) is the most common pathological type of RCC, accounting for about 70% of RCC, followed by Kidney Renal Papillary Cell Carcinoma (KIRP) and chromophobe renal cell carcinoma (1). Although surgical resection is a good treatment for renal cell carcinoma at early stage, 30% of patients are diagnosed in advanced stage due to delayed diagnosis, and 10% to 20% of patients develop metastatic kidney cancer due to postoperative recurrence, which eventually leads to death (2–4). At present, due to the lack of understanding of the pathogenesis of renal cancer, there is a lack of effective treatment for metastatic renal cancer. Therefore, to explore the pathogenesis of kidney cancer is of great significance for the treatment of kidney cancer. The search for biomarkers related to kidney cancer can help clinicians personalize patient treatment strategies and increase patient benefits.
Kidney cancer is a heterogeneous disease with multiple subtypes, multiple genes, different biochemical characteristics and multiple forms (5). KIRP, the second most common type of RCC, is a heterogeneous disease originating in the tubular epithelium of the kidney (6). The histological features of KIRP are the papillary arrangement of tumor cells, and the axis of the papilla is fibrous vascular tissue (7). In 1997, Delahunt and Eble divided KIRP into type 1 and type 2 according to histopathological characteristics and prognostic differences (8). In many studies, histological subtypes have been shown to be important prognostic predictors, with type 1 KIRP having a better prognosis than type 2 KIRP (9). Previous studies have reported that type 2 KIRP have higher nuclear grading, later staging, and poorer prognosis than type 1 KIRP (10). In contrast, Bigot et al. showed in a study of 486 patients with KIRP who underwent nephron-sparing surgery that the histological subtype of KIRP had no effect on postoperative tumor outcome (11). In conclusion, whether the histological subtype involved in type 1 or type 2 can be used as an independent prognostic factor is controversial, and correct histological phenotype and prognostic prediction are essential for the formulation of medical protocols.
Advances in sequencing technology and machine learning of all kinds have led to significant advances in the acquisition and analysis of omics data, which have deepened the understanding of tumors at the molecular level (12). Compared with a single type of data, omics data reflect the characteristics of biological individuals at multiple levels, which provides the possibility to delineate cellular molecular mechanisms in detail. Different levels of omics data reflect different relationships between genomic distribution, cancer occurrence, progression, and prognosis (13). At the same time, each omics data has its own advantages. For example, methylation chip data and lncRNA expression matrix have good tissue conservation, which can be used as efficient markers for the early diagnosis of specific tumor tissues (14). miRNA data are characterized by dissociation and can be used for non-invasive diagnosis and dynamic detection of disease (15). Common transcriptome, or mRNA sequencing, is the cheapest and most readily available, and is suitable for use in a wide range of cohort studies to explore general patterns in patient populations (16–25).
RCCs with different pathological types have different therapeutic methods and prognosis. In addition, existing targeted drugs are mainly used for ccRCCs, with unclear clinical efficacy in non- ccRCCs (26, 27). It is important to note that there is currently a lack of multi-omics prognostic molecular typing based on KIRP to guide the diagnosis and treatment of KIRP. In this study, the risk stratification of KIRP was studied by integrating multiple omics, and the differences of subgroups were analyzed in each single omics data to characterize the key events in the development of KIRP. The study provides a reference for precision medicine of KIRP.
2 Materials and methods
2.1 Extraction and preprocessing of multi-omics data for KIRP
The dataset for KIRP was downloaded from The Cancer Genome Atlas (TCGA) (28) and TCGA database had the multi-omics data for our analysis in this study. We acquired gene expression profile for transcriptomics (including mRNAs encoding protein, long noncoding RNAs as known as lncRNAs, microRNA known as miRNAs, data of methylation, and data of gene mutations). We applied TCGAbiolinks package of R application to acquire clinicopathologic information and multi omics-based data. We downloaded the gene expression profiles of 34 cases with KIRP from GSE2748 as the external validation cohort (29). The patients with KIRP in the GSE2748 dataset had the prognostic information (29). In addition, we searched and downloaded the single-cell RNA sequencing for KIRP from GSE152938 (30). There was a total of four KIRP samples and one normal kidney sample included in GSE152938 (30). The matrix for single-cell RNA sequencing was generated by R package Seurat (31).
2.2 Identification of multi omics-based cancer subtypes by integrative analysis
The MOVICS package aimed to show the multi-omics comprehensive clustering and visualization of cancer typing studies (32). There were ten algorithms included in the MOVICS package: CIMLR, iClusterBayes, MoCluster, COCA, ConsensusClustering, IntNMF, LRAcluster, NEMO, PINSPlus and SNF (32). For the multi omics-based data, we focused on the characteristic related to prognosis (OS). The OS-related features, including mRNAs, lncRNAs, miRNAs, methylation, and gene mutations were analyzed by Univariate Cox regression analysis, and we screened out features with the threshold P-value<0.05. Due to the small amount of mutation matrix and miRNA expression data, only the top 30 mutations and 200 miRNA data were extracted. We carried out analysis for Clustering Prediction Index (CPI) (33) and Gaps-statistics (33) to filtrate out the optimal number of cancer subtypes. We finally identified the multi omics-based cancer subtype (MOCS) based on consensus ensembles and high robustness, thus separating the patients with KIRP into different MOCSs.
2.3 Nearest template prediction validation
Nearest template prediction (NTP) algorithm could also be applied to cross-platform, cross-species and multi-class predictions without any optimization of analysis parameters (34). In this study, we also used NTP algorithm of CMScaller package to test the dependability and stability of MOCS subtypes via the external GSE2748 cohort.
2.4 Biological characteristics for MOCS subtypes
The gene sets (including immune-related pathways) were analyzed, and enrichment scores were calculated using gene set variation analysis (GSVA) from R package GSVA (35). The differentially expressed genes (DEGs) among the three MOCS subtypes were assessed using limma package (36). Pathway enrichment analysis was performed by clusterProfiler package with the employment of Biological Processes in Gene Ontology (GO) (37).
2.5 Calculation of immune microenvironmental cell abundance and immune checkpoint
Tumor Immune Estimation Resource (TIMER) is a website from which researchers can use RNA-Seq expression profile data to detect the infiltration of immune cells in tumor tissue (38). The TIMER provides the infiltrations of six kinds of immune cells (B cells, CD4+ T cells, CD8+ T cells, Neutrophil, Macrophages and Myeloid dendritic cells) (38). CIBERSORT (39) and CIBERSORT-ABS (40) algorithms were used to acquire the infiltrations of 22 kinds of immune cells. quanTIseq is a deconvolution tool developed specifically for RNA-seq data, enabling accurate quantification of unknown tumor content, as well as the immune cell component of the overall tissue (41). quanTIseq implemented a complete deconvolution process for analyzing RNA-seq data based on constrained least squares regression and a new eigenmatrix from 51 purified or enriched RNA-seq data sets, avoiding inconsistencies between mixtures and eigenmatrices (41). MCPcounter (42), xCell (43), and EPIC (44) algorithms (Estimate the Proportion of Immune and Cancer cells) were also used to assess the immune microenvironmental cell abundance. We estimated the infiltrating level of immune or stromal scores using ESTIMATE R package (45). In addition, DNA methylation of tumor-infiltrating lymphocyte (MeTIL) for TCGA- KIRP cohort was also calculated (46).
2.6 Evaluation of genomic mutation for MOCS subtypes
Mutation profiles of KIRP were acquired and we compared and visualized the difference of mutation among the MOCS subtypes utilizing Maftools package of R (47). We applied the Maftools function to analyze the oncogenic pathway and mutually exclusive or coexisting mutations (48). The loss and gain in genomic level was evaluated by GISTIC 2.0 algorithm (49).
2.7 Drug sensitivity profiles for MOCS subtypes
R package pRRophetic was employed to predict the drug sensitivity profiles for MOCS subtypes (50, 51). Subclass mapping was used to explore the immunotherapy of KIRP based on the literature published (52, 53).
2.8 Statistical analyses
R was used to conduct statistical analyses (v4.0.2). We also provide the codes of all methods used in this paper in Supplementary Code. P values or adjusted P values less than 0.05 were considered significant for all statistical comparisons.
3 Results
3.1 Three MOCSs were categorized for KIRP patients by MOVICS package
We discovered three MOCS subtypes for KIRP patients based on CPI analysis and Gaps-statistics, due to the optimal average statistic value with the number of MOCSs was found to be three (Figure 1A). Hence, patients with KIRP were divided into MOCS1, MOCS2, and MOCS3, indicating the robustness of the classification system (Figure 1B). The silhouette plot indicated that the silhouette score of MOCS1 was 0.70, the silhouette score of MOCS2 was 0.42, while the silhouette score of MOCS3 was 0.71, which substantiated that the MOCS subtypes were distinguishable and separated well from each other (Supplementary Figure S1A). From the Figure 1C, the consistency of the classification system for MOCSs were observed in consideration of four statistics (54) (Figure 1C), including Rand Index (RI), Adjusted Mutual Information (AMI), Jaccard Index (JI), and Fowlkes-Mallows (FM) (54). Furthermore, we observed that all the patients in the MOCS1 and MOCS1 were at AJCC Stage I and Pstage I (Figure 1C). Among the three MOCSs, we displayed the distribution of the multi-omics data for mRNA, lncRNA, miRNA, DNA methylation, and gene mutations as shown in the heatmap (Figure 1D). In the distribution diagram (Figure 1D), RBP4, MSLN, VSTM2L, FTCD, AC147651.5, RP11-23P13.6, RP11-326C3.2, RP11-124N19.3, CHL1-AS2, and RP11-807H17.1 were the top ten OS-related factors of transcriptome (mRNAs and lncRNAs). As for miRNA, hsa-mir-127, hsa-mir-1247, hsa-mir-1-1, hsa-mir-1-2, hsa-mir-1180, hsa-mir-1269a, hsa-mir-10b, hsa-mir-126, hsa-mir-105-1, and hsa-mir-105-2 were the top ten OS-related factors of miRNAs (Figure 1D). As for DNA methylation, cg16434331, cg06775420, cg25244238, cg06282596, cg02239902, cg22688012, cg23591302, cg03994717, cg06223834, and cg06234051 were the top ten OS-related factors (Figure 1D). SETD2, PBRM1, SYNE2, NF2, MET, LRP2, CUL3, PKHD1, TTN, and PCF11 were the top ten OS-related factors (Figure 1D). Further, we compared the outcome of clinical prognosis of patients with KIRP among MOCS1, MOCS2, and MOCS3. Patients in the MOCS3 subtype were observed with shorter overall survival time than patients in the MOCS1 and MOCS2 subtypes (Figure 1E), which was also observed for progression free survival time (Figure 1E). Using NTP algorithm, three MOCSs were also identified as predicted by the external GSE2748 cohort (Supplementary Figure S1B). Patients in the MOCS3 subtype were observed with shorter overall survival time (Supplementary Figure S1C).
Figure 1 Three MOCSs were categorized for KIRP patients by MOVICS package. (A) Determination of optimal cluster number through calculating CPI (blue line) and Gaps-statistics (red line) in TCGA- KIRP cohort. (B) Consensus heatmap based on outcomes from 10 multi-omics integrative clustering approaches with subtype number of three showing perfect diagonal rectangle. (C) Quantification of sample similarity using silhouette score based on the consensus ensembles result and alluvial diagram presenting the flow distribution among different multi omics-based cancer subtypes (MOCSs). (D) Comprehensive heatmap showing the detailed molecular landscape multi-omics data for mRNA, lncRNA, miRNA, DNA methylation, and gene mutations among the three MOCSs. (E) log-rank test for overall survival time and progression free survival time for patients with KIRP.
3.2 Biological characteristics for MOCS subtypes
Further, we depicted the molecular features characterization for MOCS subtypes. We computed the enrichment score of immune-related pathways (including Cell Functions, B Cell Functions, T Cell Functions, Leukocyte Functions, Pathogen Defense, Interleukins, TNF Superfamily, Chemokines, Cytokines, Regulation NK Cell Functions Complement, Antigen Processing, Cytotoxicity, Microglial Functions, TLR, Adhesion, Transporter Functions, Cell Cycle, Macrophage Functions and Senescence) based on GSVA analysis. We could find that MOCS1was negatively correlated with immune-related pathways (Figure 2A). As for other pathways, we found global dysfunction of cancer-related pathways among the three MOCS subtypes (Figure 2B). Generally, MOCS1 showed relatively lower enrichment level of Nature metabolism Hypoxia, Hu hypoxia signature, Exosomal secretion, Ferroptosis, MT exosome and exosome assembly (Figure 2B), suggesting the three MOCS subtypes were association with exosomes strongly. Biological processes of AXONEMAL DYNEIN COMPLEX ASSEMBLY, CILIUM MOVEMENT, AXONEME ASSEMBLY, INNER DYNEIN ARM ASSEMBLY, INTRACILIARY TRANSPORT, MICROTUBULE BUNDLE FORMATION, PROTEIN LOCALIZATION TO CILIUM, MRNA SPLICE SITE SELECTION, INTRACILIARY TRANSPORT INVOLVED IN CILIUM ASSEMBLY, and EXTRACELLULAR TRANSPORT were overactivated in MOCS1 (Figure 2C). Biological processes of RESPIRATORY ELECTRON TRANSPORT CHAIN, RESPIRATORY ELECTRON TRANSPORT CHAIN, ELECTRON TRANSPORT CHAIN, ORGANIC ACID CATABOLIC PROCESS, OXIDATIVE PHOSPHORYLATTON, GOTATP SYNTHESIS COUPLED ELECTRON TRANSPORT, COFACTOR METABOLIC PROCESS, COENZYME METABOLIC PROCESS, AEROBIC RESPIRATION, SMALL MOLECULE CATABOLIC PROCESS, and ALPHA AMINO ACID METABOLIC PROCESS (Figure 2C). Biological processes of CORNIFICATION, NEURON FATE SPECIFICATION, TONGUE DEVELOPMENT, AUTONOMIC NERVOUS SYSTEM DEVELOPMENT, INNERVATION, FORELIMB MORPHOGENESIS, APPENDAGE DEVELOPMENT, ENDOCARDIAL CUSHION MORPHOGENESIS, APPENDAGE MORPHOGENESIS, EYELID DEVELOPMENT IN CAMERA TYPE EYE (Figure 2C). In addition, we evaluated the activity profiles of regulons among the three MOCSs, thus highlighting the additional potential regulatory differences. The higher level of several regulon, such as ZNF683, IRF4, CEBPB, EPAS1, and TFE3 was observed in MOCS2 and MOCS3 (Figure 2D), indicating the important differentiators of epigenetically driven transcriptional networks among the three MOCS subtypes. GSVA analysis was carried out regarding metabolism-related pathways, we found that most of the metabolism-related pathways were activated in MOCS2 (Supplementary Figure S2A). Consistently, most of immune-associated signatures were enriched in MOCS2 (Supplementary Figure S2B).
Figure 2 Biological characteristics for MOCS subtypes. (A) Heatmap showing GSVA enrichment score of immune-related pathways among the three MOCS subtypes. (B) Heatmap showing GSVA enrichment score of cancer-related pathways among the three MOCS subtypes. (C) GO enrichment analysis showing the upregulated pathways and the downregulated pathways. (D) Heatmap showing the regulon distribution among the three MOCS subtypes.
3.3 Calculation of immune microenvironmental cell abundance and immune checkpoint
In consideration of the critical role of immunity in KIRP progression, we investigated the immune microenvironmental cell abundance and immune checkpoint among the three MOCS subtypes. Several immune microenvironmental cells were highly infiltrated in specific MOCS subtype. For instance, B cell in MOCS3, Macrophage M2 in MOCS2, NK cell in MOCS1 and so on (Figure 3A). As for the immune checkpoint genes, on the whole, MOCS3 was associated with higher levels of immune checkpoint genes (Figure 3B). MOCS3 was also associated with higher levels of MeTIL (Figure 3B). MOCS1 was found to be associated with lower levels of immune, stromal and ESTIMATE scores (Figure 3C). Additionally, we found that MOCS3 showed a significantly lower tumor mutation burden (TMB, Figure 3D). MOCS1 was found to be associated with lower signature score of CD8+ T effector, Immune checkpoint, APM, TME score A, Pan F TBRs, EMT2, EMT3, and TME score B (Supplementary Figure S3A). The level of RNAss, DMPss, ENHss, EREG.EXPss and HRD was found to be lower in MOCS1 (Supplementary Figure S3B).
Figure 3 Calculation of immune microenvironmental cell abundance and immune checkpoint. (A) Heatmap showing the immune microenvironmental cell profile for TCGA- KIRP cohort based on TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, EPIC algorithms. (B) Heatmap showing the profile for immune checkpoint genes, and DNA methylation of tumor-infiltrating lymphocytes (MeTILs). (C) Boxplot showing the distribution of immune, stromal and ESTIMATE scores. (D) Distribution of TMB and TiTv (transition to transversion) among the three MOCS subtypes. *** means P < 0.001 and **** means P < 0.0001.
3.4 Evaluation of genomic mutation for MOCS subtypes
The differences in copy number variations (CNV) among the three MOCS subtypes were compared, and the result revealed that the CNV occurrence frequency was higher in MOCS1 (Figures 4A–C). In detail, amplification in chr 2p, 2q, 3p, 3q, 4p, 7p, 7q, 12p, 12q, 16p, 16q, 17p, 17q, 18p, 18q, 20p, 20q, and 21q were higher in MOCS1 (Figure 4B). The above results were also proved by the total copy number alteration rate as shown in Figure 4C. MOCS1 displayed a higher rate in focal and arm-level mutation level gain (Figure 4D). Mutation patterns of the top 20 most frequently mutated genes among the three MOCSs were displayed in the waterfall plot (Supplementary Figure S4A), from which we could see that, TTN, MET, CUBN, SYNE1, HERC2, KIAA1109, MUC16, PKHD1, WDFY3, DNAH8, KMT2C, LRP2, MACF1, NEB, PCLO, SMARCA4, ANK3, COL18A1, DDX5, DYNC2H1 were the top 20 mutated genes for MOCS1; TTN, SETD2, MUC16, CUL3, KIAA1109, KMT2C, PBRM1, PCF11, BAP1, FAT1, KMT2D, PKHD1, KDM6A, LRBA, SRRM2, ARID1A, ASAP2, BIRC6, CENPE and CNOT1 were the top 20 mutated genes for MOCS2; NF2, TTN, TXNIP, BAP1, CAMK1D, CDH8, CMYA5, CREBBP, EBF2, HECTD4, ITGAL, KRAS, MAP1B, TAS1R2, TG, EIF4G3, FAT1, HELZ2, KDM6A, and SYNE1 were the top 20 mutated genes for MOCS3 (Supplementary Figure S4A). The synthetic lethal mutations in MOCS1, MOCS2, and MOCS3 were displayed in Supplementary Figure S4B. The potential druggable gene categories from the mutation dataset for MOCS1, MOCS2, and MOCS3 were shown in Supplementary Figure S4C, we found that ANK3, CUBN, LRP2, MET, PKHD1 and so on were the potential therapeutic targets for MOCS1; ARID1A, BAP1, CUL3, FAT1, KDM6A and so on were the potential therapeutic targets for MOCS2; BAP1, CAMK1D, CMYA5, CREBBP, FAT1 and so on were the potential therapeutic targets for MOCS3. The fraction of pathways and samples affected were the minimum among the three MOCSs (Supplementary Figure S4C).
Figure 4 Landscapes of copy number variations. (A) Comparison of overall copy number among all patients with KIRP, MOCS1, MOCS2, and MOCS3. (B) The amplification or deletion frequency in chromosome among the three MOCSs. (C) Bar-plot indicating the total alteration frequency among the three MOCSs. (D) Different burden of copy number gain at focal and arm-level among the three MOCSs. * means P < 0.05, ** means P < 0.01, *** means P < 0.001, and **** means P < 0.0001, ns, no significance.
3.5 Drug sensitivity profiles for MOCS subtypes
We collected drug response data reflected by the IC50 value via GDSC database. We observed that patients in MOCS3 were more sensitive to Crizotinib, Erlotinib, Pazopanib, Saracatinib, Sunitinib, and Temsirolimus (Figure 5A). We found that patients in MOCS1 were more sensitive to AS601245, Bosutinib, PAC.1, ABT.888, and Bleomycin (Figure 5B). Whereafter, we carried out subclass mapping and the results revealed that patients in MOCS2 were more likely to respond to anti-PD1 blockades (Figure 5C).
Figure 5 Drug Sensitivity Profiles for MOCS subtypes. (A) Estimated IC50 of Crizotinib, Erlotinib, Pazopanib, Saracatinib, Sunitinib, and Temsirolimus among the three MOCSs. (B) Estimated IC50 of CCT018159, AS601245, Bosutinib, PAC.1, DMOG, BMS.708163, A.443654, Sunitinib, ABT.888, and Bleomycin. (C) Subclass analysis manifested that MOCS2 were more likely to respond to anti-PD1 blockades. * means P < 0.05 and ** means P < 0.01.
3.6 Single-cell analysis
A total of 16 cell clusters were identified after gene filtering, normalization and principal component analysis, as shown in Figure 6A. There were nine specific cell types, including B cell, CD8+ T cell, Endothelial cell, Plasma cell, TAM cell, CAF cell, Dendritic cell, Fibroblast cell, pRCC cell (Figure 6B). In addition, a total of three cell clusters (C0, C1, and C2) were predicted by Scissor tool, as shown in Figure 6C. The bar graph displayed the fraction of specific cell types in each cell cluster predicted by Scissor tool (Figure 6D). C0 cluster was rich in TAM cell, CAF cell, Fibroblast cell, CD8+ T cell, Endothelial cell, and pRCC cell (Figure 6D). C1 cluster was rich in Dendirtic cell, Plasma cell, and B cell (Figure 6D). The correlation networks were generated to show the interactions among different cells (Figure 6E). The ligand–receptor pairs among cells were displayed in Figure 6F.
Figure 6 Cell–cell interactions for KIRP. (A) tSNE plot of the distribution of 16 samples; (B) tSNE plot of the distribution of nine cell clusters after clustering. (C) UMAP plot showing three subclusters (C0, C1, C2) of the KIRP. (D) The fraction of specific cell types in each cell cluster predicted by Scissor tool. (E) Circle plot showing the intercellular communication among major cell types in KIRP. (F) The ligand–receptor pairs among cells.
4 Discussion
Global cancer data show that RCC accounts for about 3%~5% of adult malignant tumors, and its incidence is higher in males than in females (4, 55). RCC is the 9th most common male cancer and the 14th most common female cancer worldwide (56). As for etiology, tobacco exposure of any kind is thought to be associated with the development of kidney cancer (57). In addition, diets high in fat, high in protein, low in fruits and vegetables, and increased intake of dairy products are associated with kidney cancer, but the relative risk is not very high, and many scholars have different views (57, 58). The current diagnosis and treatment problem of renal cancer is the heterogeneity of the tumor, which often leads to different prognosis of patients with the same stage and grade (59). In addition, the incidence of tumor resistance and metastasis is high in renal cancer, and the treatment options for these patients are extremely limited, resulting in a low 5-year survival rate (60). In view of the above problems, it is necessary to determine new diagnosis and treatment strategies to improve the survival rate of patients with kidney cancer.
TCGA is an oncology research initiative of The Cancer Genome Atlas and the National Human Genome Research Institute (28). The plan includes multifactorial data on common tumor tissues and prognostic information for patients. The data included pathological sections, cancer and para-cancer transcriptome, methylation chip data and genome data (28). The development of multi-omics has made it easy for researchers to deepen their understanding of cancer at the molecular level. At the same time, a large number of omics data also brings new challenges to analysts (61). It is particularly critical to reduce data noise and obtain key characteristics of tumor occurrence and development while preserving tumor characteristics (61). Few studies have attempted to establish a comprehensive model based on multiple omics data to predict prognosis and personalized drug selection in patients with KIRP. Therefore, it is particularly important to develop a comprehensive and robust prognostic and drug selection model for patients with KIRP to assist in prognostic prediction and guide personalized treatment. In this study, we conducted a comprehensive integrated analysis of multiple omics data, including mRNA, lncRNA, miRNA, DNA methylation profile and somatic mutation data, and constructed a classifier to evaluate the prognosis of patients with KIRP and assist drug selection. Omics data are complex, multi-layered, and high weaves, so a key goal of analyzing multi-omics data is to screen for valid predictors to predict phenotypic characteristics and thus elucidate the biological significance behind them. Another major difficulty in omics data processing is dimensionality reduction, omics noise elimination and overfitting avoidance. In this study, R package Survival was first used to screen the molecular features associated with patient prognosis in each omics for subsequent analysis. The classification of cancer patients into different molecular subgroups based on multi-omics data is an important problem in the context of precision medicine. MOVICS provides a unified interface to 10 state-of-the-art multiomics ensemble clustering algorithms and integrates the downstream analyses most commonly used in cancer typing studies, including characterization and comparison of identified subtypes from multiple perspectives and validation of subtypes in external corporations of multi-class predictions using two model-free methods. Patients with KIRP were divided into three multi omics-based cancer subtypes (MOCS1, MOCS2, and MOCS3). Patients in the MOCS3 subtype were observed with shorter overall survival time than patients in the MOCS1 and MOCS2 subtypes, therefore, the classification system can be used as an important prognostic tool. Similar prognostic outcomes were observed in independent external datasets. Therefore, the classification system established by us is reliable in prognostic assessment.
In recent years, molecular typing of kidney cancer has been emerging. Molecular typing of renal cancer from genomic changes, DNA methylation profiles, RNA and protein levels has revealed repeated mutations in the PI3K/AKT pathway, suggesting that this pathway is a potential therapeutic target (62). A large number of molecular typing studies of renal cancer have emerged based on single omics or specific gene sets, recently. Chen et al. integrated multi-omics data of all kidney cancer patients based on a single algorithm but did not include data of lncRNA data in the analysis (63). Ricketts et al. integrated the multi-omics data of kidney cancer for reclassification, while this study only conducted classification from the level of each omics, without realizing the real sense of integrated multi-omics data for classification (64). Although these studies provide new directions for the diagnosis and treatment of kidney cancer to some extent, they also have certain shortcomings. The classification methods used in most typing studies are relatively simple. These shortcomings make it difficult to apply these classification studies to clinical practice. In this study, ten robust clustering algorithms based on MOVICS package were used, combined with multiple omics information, to conduct multi-omics cross-validation for patients with KIRP. Further, intra-omics heterogeneity analysis was conducted at each omics level to crack the omics differences among patients with different prognostic characteristics. Specifically, patients in the MOCS3 subtype were observed with shorter overall survival time than patients in the MOCS1 and MOCS2 subtypes. Compared to the other two subtypes, MOCS1was negatively correlated with immune-related pathways. Global dysfunction of cancer-related pathways among the three MOCS subtypes were also observed. We also evaluated the immune microenvironmental cell abundance and immune checkpoint and compared the discrepancy among the MOCSs.
Our study unexpectedly found that these three MOCSs also have significant differences in sensitivity to molecularly targeted drugs. We observed that patients in MOCS3 were more sensitive to Crizotinib, Erlotinib, Pazopanib, Saracatinib, Sunitinib, and Temsirolimus; while patients in MOCS1 were more sensitive to AS601245, Bosutinib, PAC.1, ABT.888, and Bleomycin. Whereafter, the results of subclass mapping revealed that patients in MOCS2 were more likely to respond to anti-PD1 blockades. In recent years, the treatment of kidney cancer has evolved from non-specific immune approaches to targeted therapy of vascular endothelial growth factor (VEGF), and now to novel immunotherapies. Our study assessed therapeutic differences among different subtypes and therefore can be a potential therapeutic direction for patients with KIRP.
In summary, our study provides a new reference for molecular subtypes of KIRP risk. In this study, a robust prognostic and drug selection subtype system was constructed by integrating multiple omics data using multiple algorithms. However, there are still some limitations in our study. Firstly, multi-omics data used for molecular subtypes is difficult to be applied in clinical practice. Second, although we compared the enrichment pathway and drug sensitivity between subgroups, further experiments and external data sets are still needed for verification.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Author contributions
RL designed the study. BW and ML performed data analysis. BW drafted the manuscript. RL revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by the Provincial Special Supporting Fund Scientific Research Project of Shanxi Provincial People’s Hospital (No. SZ2019022 and No. XY2018011), the Shanxi Fundamental Research Program (No. 201801D221382), and the Project of Administration of Traditional Chinese Medicine of Shanxi Province (2023ZYYA019).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2023.1169395/full#supplementary-material
Supplementary Figure 1 | Identifying and verifying of the multi omics-based cancer subtype. (A) Quantification of sample similarity using silhouette score based on the consensus ensembles result. (B) Three MOCSs were also identified as predicted by the external GSE2748 cohort. (C) Comparison of the overall survival time for the three MOCSs.
Supplementary Figure 2 | Functional enrichment analysis of MOCS1, MOCS2 and MOCS3 subgroups. (A) Heatmap of metabolism-related enrichment scores among the three MOCSs. (B) Heatmap of immune-related enrichment scores among the three MOCSs.
Supplementary Figure 3 | The landscapes of specific immune scores among the three MOCSs. (A) Comparison of the signature score among the three MOCSs. (B) Comparison of the RNAss, DNAss, DMPss, ENHss, EREG.EXPss, and HRD among the three MOCSs.
Supplementary Figure 4 | Landscapes of somatic mutations and potential targets in the two subtypes. (A) Waterfall plot showing the mutation patterns of the top 20 most frequently mutated genes among the three MOCSs. (B) The synthetic lethal mutations in MOCS1, MOCS2, and MOCS3. (C) Potential druggable gene categories from the mutation dataset for MOCS1, MOCS2, and MOCS3. (D) The fraction of pathways or samples of oncogenic signaling pathways for MOCS1, MOCS2, and MOCS3.
References
1. Motzer RJ, Jonasch E, Michaelson MD, Nandagopal L, Gore JL, George S, et al. NCCN guidelines insights: Kidney cancer, version 2.2020. J Natl Compr Canc Netw (2019) 17:1278–85. doi: 10.6004/jnccn.2019.0054
2. Shinder BM, Rhee K, Farrell D, Farber NJ, Stein MN, Jang TL, et al. Surgical management of advanced and metastatic renal cell carcinoma: A multidisciplinary approach. Front Oncol (2017) 7:107. doi: 10.3389/fonc.2017.00107
3. McGuire S. World cancer report 2014. In: Adv nutr, vol. 7. . Geneva, Switzerland: World Health Organization, International Agency for Research on Cancer, WHO Press, 2015 (2016). p. 418–9. doi: 10.3945/an.116.012211
4. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin (2022) 72:7–33. doi: 10.3322/caac.21708
6. Thoenes W, Störkel S, Rumpelt HJ. Histopathology and classification of renal cell tumors (adenomas, oncocytomas and carcinomas). the basic cytological and histopathological elements and their use for diagnostics. Pathol Res Pract (1986) 181:125–43. doi: 10.1016/s0344-0338(86)80001-2
7. Troxell ML, Higgins JP. Renal cell carcinoma in kidney allografts: Histologic types, including biphasic papillary carcinoma. Hum Pathol (2016) 57:28–36. doi: 10.1016/j.humpath.2016.06.018
8. Delahunt B, Eble JN. Papillary renal cell carcinoma: A clinicopathologic and immunohistochemical study of 105 tumors. Mod Pathol (1997) 10:537–44.
9. Allory Y, Ouazana D, Boucher E, Thiounn N, Vieillefond A. Papillary renal cell carcinoma. prognostic value of morphological subtypes in a clinicopathologic study of 43 cases. Virchows Arch (2003) 442:336–42. doi: 10.1007/s00428-003-0787-1
10. Cancer Genome Atlas Research Network, Linehan WM, Spellman PT, Ricketts CJ, Creighton CJ, Fei SS, et al. Comprehensive molecular characterization of papillary renal-cell carcinoma. N Engl J Med (2016) 374:135–45. doi: 10.1056/NEJMoa1505917
11. Bigot P, Bernhard JC, Gill IS, Vuong NS, Verhoest G, Flamand V, et al. The subclassification of papillary renal cell carcinoma does not affect oncological outcomes after nephron sparing surgery. World J Urol (2016) 34:347–52. doi: 10.1007/s00345-015-1634-0
12. Sánchez-Gastaldo A, Kempf E, González Del Alba A, Duran I. Systemic treatment of renal cell cancer: A comprehensive review. Cancer Treat Rev (2017) 60:77–89. doi: 10.1016/j.ctrv.2017.08.010
13. Dimitrieva S, Schlapbach R, Rehrauer H. Prognostic value of cross-omics screening for kidney clear cell renal cancer survival. Biol Direct (2016) 11:68. doi: 10.1186/s13062-016-0170-1
14. Mariani MP, Chen JA, Zhang Z, Pike SC, Salas LA. MethylMasteR: A comparison and customization of methylation-based copy number variation calling software in cancers harboring Large scale chromosomal deletions. Front Bioinform (2022) 2:859828. doi: 10.3389/fbinf.2022.859828
15. Miranda-Poma J, Trilla-Fuertes L, López-Camacho E, Zapater-Moros A, López-Vacas R, Lumbreras-Herrera MI, et al. MiRNAs in renal cell carcinoma. Clin Transl Oncol (2022) 24:2055–63. doi: 10.1007/s12094-022-02866-z
16. Mishra SK, Thakran P. Intron specificity in pre-mRNA splicing. Curr Genet (2018) 64:777–84. doi: 10.1007/s00294-017-0802-8
17. Wahl MC, Will CL, Lührmann R. The spliceosome: Design principles of a dynamic RNP machine. Cell (2009) 136:701–18. doi: 10.1016/j.cell.2009.02.009
18. Sillekens PT, Beijer RP, Habets WJ, van Verooij WJ. Molecular cloning of the cDNA for the human U2 snRNA-specific a’ protein. Nucleic Acids Res (1989) 17:1893–906. doi: 10.1093/nar/17.5.1893
19. Kim YD, Lee J, Kim HS, Lee MO, Son MY, Yoo CH, et al. The unique spliceosome signature of human pluripotent stem cells is mediated by SNRPA1, SNRPD1, and PNN. Stem Cell Res (2017) 22:43–53. doi: 10.1016/j.scr.2017.05.010
20. Liu Y, Zhu X, Zhu J, Liao S, Tang Q, Liu K, et al. Identification of differential expression of genes in hepatocellular carcinoma by suppression subtractive hybridization combined cDNA microarray. Oncol Rep (2007) 18:943–51. doi: 10.3892/or.18.4.943
21. Dou N, Yang D, Yu S, Wu B, Gao Y, Li Y, et al. SNRPA enhances tumour cell growth in gastric cancer through modulating NGF expression. Cell Prolif (2018) 51:e12484. doi: 10.1111/cpr.12484
22. Fish L, Khoroshkin M, Navickas A, Garcia K, Culbertson B, Hänisch B, et al. A prometastatic splicing program regulated by SNRPA1 interactions with structured RNA elements. Science (2021) 372(6543):eabc7531. doi: 10.1126/science.abc7531
23. Bolduc F, Turcotte MA, Perreault JP. The small nuclear ribonucleoprotein polypeptide a (SNRPA) binds to the G-quadruplex of the BAG-1 5’UTR. Biochimie (2020) 176:122–7. doi: 10.1016/j.biochi.2020.06.013
24. Zeng Q, Lei F, Chang Y, Gao Z, Wang Y, Gao Q, et al. An oncogenic gene, SNRPA1, regulates PIK3R1, VEGFC, MKI67, CDK1 and other genes in colorectal cancer. BioMed Pharmacother (2019) 117:109076. doi: 10.1016/j.biopha.2019.109076
25. Yuan P, Ling L, Gao X, Sun T, Miao J, Yuan X, et al. Identification of RNA-binding protein SNRPA1 for prognosis in prostate cancer. Aging (Albany NY) (2021) 13:2895–911. doi: 10.18632/aging.202387
26. Négrier S, Gravis G, Pérol D, Chevreau C, Delva R, Bay JO, et al. Temsirolimus and bevacizumab, or sunitinib, or interferon alfa and bevacizumab for patients with advanced renal cell carcinoma (TORAVA): A randomised phase 2 trial. Lancet Oncol (2011) 12:673–80. doi: 10.1016/s1470-2045(11)70124-3
27. Ljungberg B, Bensalah K, Canfield S, Dabestani S, Hofmann F, Hora M, et al. EAU guidelines on renal cell carcinoma: 2014 update. Eur Urol (2015) 67:913–24. doi: 10.1016/j.eururo.2015.01.005
28. Wang Z, Jensen MA. & zenklusen, j. c. a practical guide to the cancer genome atlas (TCGA). Methods Mol Biol (2016) 1418:111–41. doi: 10.1007/978-1-4939-3578-9_6
29. Yang XJ, Tan MH, Kim HL, Ditlev JA, Betten MW, Png CE, et al. A molecular classification of papillary renal cell carcinoma. Cancer Res (2005) 65:5628–37. doi: 10.1158/0008-5472.Can-05-0533
30. Su C, Lv Y, Lu W, Yu Z, Ye Y, Guo B, et al. Single-cell RNA sequencing in multiple pathologic types of renal cell carcinoma revealed novel potential tumor-specific markers. Front Oncol (2021) 11:719564. doi: 10.3389/fonc.2021.719564
31. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol (2015) 33:495–502. doi: 10.1038/nbt.3192
32. Lu X, Meng J, Zhou Y, Jiang L, Yan F. MOVICS: An r package for multi-omics integration and visualization in cancer subtyping. Bioinformatics (2020) 36(22–23):5539–41. doi: 10.1093/bioinformatics/btaa1018
33. Hu X, Wang Z, Wang Q, Chen K, Han Q, Bai S, et al. Molecular classification reveals the diverse genetic and prognostic features of gastric cancer: A multi-omics consensus ensemble clustering. BioMed Pharmacother (2021) 144:112222. doi: 10.1016/j.biopha.2021.112222
34. Hoshida Y. Nearest template prediction: A single-sample-based flexible class prediction with confidence assessment. PloS One (2010) 5:e15543. doi: 10.1371/journal.pone.0015543
35. Hu J, Tzeng JY. Integrative gene set analysis of multi-platform data with sample heterogeneity. Bioinformatics (2014) 30:1501–7. doi: 10.1093/bioinformatics/btu060
36. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007
37. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
38. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res (2017) 77:e108–10. doi: 10.1158/0008-5472.Can-17-0307
39. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12
40. Zhu K, Xiaoqiang L, Deng W, Wang G, Fu B. Development and validation of a novel lipid metabolism-related gene prognostic signature and candidate drugs for patients with bladder cancer. Lipids Health Dis (2021) 20:146. doi: 10.1186/s12944-021-01554-1
41. Plattner C, Finotello F, Rieder D. Deconvoluting tumor-infiltrating immune cells from RNA-seq data using quanTIseq. Methods Enzymol (2020) 636:261–85. doi: 10.1016/bs.mie.2019.05.056
42. Sturm G, Finotello F, List M. Immunedeconv: An r package for unified access to computational methods for estimating immune cell fractions from bulk RNA-sequencing data. Methods Mol Biol (2020) 2120:223–32. doi: 10.1007/978-1-0716-0327-7_16
43. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18:220. doi: 10.1186/s13059-017-1349-1
44. Racle J, Gfeller D. EPIC: A tool to estimate the proportions of different cell types from bulk gene expression data. Methods Mol Biol (2020) 2120:233–48. doi: 10.1007/978-1-0716-0327-7_17
45. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612
46. Jeschke J, Bizet M, Desmedt C, Calonne E, Dedeurwaerder S, Garaud S, et al. DNA Methylation-based immune response signature improves patient diagnosis in multiple cancers. J Clin Invest (2017) 127:3090–102. doi: 10.1172/jci91095
47. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28:1747–56. doi: 10.1101/gr.239244.118
48. Jiang A, Song J, Fang X, Fang Y, Wang Z, Liu B, et al. A novel thinking: DDR axis refines the classification of ccRCC with distinctive prognosis, multi omics landscape and management strategy. Front Public Health (2022) 10:1029509. doi: 10.3389/fpubh.2022.1029509
49. Cancer Genome Atlas Research Network. Comprehensive molecular characterization of gastric adenocarcinoma. Nature (2014) 513:202–9. doi: 10.1038/nature13480
50. Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and. Vitro Drug sensitivity Cell lines. Genome Biol (2014) 15:R47. doi: 10.1186/gb-2014-15-3-r47
51. Geeleher P, Cox N, Huang RS. pRRophetic: An r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9:e107468. doi: 10.1371/journal.pone.0107468
52. McGranahan N, Furness AJ, Rosenthal R, Ramskov S, Lyngaa R, Saini SK, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science (2016) 351:1463–9. doi: 10.1126/science.aaf1490
53. Lu X, Jiang L, Zhang L, Zhu Y, Hu W, Wang J, et al. Immune signature-based subtypes of cervical squamous cell carcinoma tightly associated with human papillomavirus type 16 expression, molecular features, and clinical outcome. Neoplasia (2019) 21:591–601. doi: 10.1016/j.neo.2019.04.003
54. Zhang G, Liu Z, Deng J, Liu L, Li Y, Weng S, et al. Smooth muscle cell fate decisions decipher a high-resolution heterogeneity within atherosclerosis molecular subtypes. J Transl Med (2022) 20:568. doi: 10.1186/s12967-022-03795-9
55. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2021. CA Cancer J Clin (2021) 71:7–33. doi: 10.3322/caac.21654
56. Fernando A. Time to focus on the rare-encouraging progress in the management of non-clear cell renal cell carcinoma. Eur Urol (2017) 72:984–5. doi: 10.1016/j.eururo.2017.07.028
57. Scelo G, Larose TL. Epidemiology and risk factors for kidney cancer. J Clin Oncol (2018) 36:Jco2018791905. doi: 10.1200/jco.2018.79.1905
58. Gray RE, Harris GT. Renal cell carcinoma: Diagnosis and management. Am Fam Physician (2019) 99:179–84.
59. Ohashi R, Schraml P, Batavia A, Angori S, Simmler P, Rupp N, et al. Allele loss and reduced expression of CYCLOPS genes is a characteristic feature of chromophobe renal cell carcinoma. Transl Oncol (2019) 12:1131–7. doi: 10.1016/j.tranon.2019.05.005
60. Lalani AA, McGregor BA, Albiges L, Choueiri TK, Motzer R, Powles T, et al. Systemic treatment of metastatic clear cell renal cell carcinoma in 2018: Current paradigms, use of immunotherapy, and future directions. Eur Urol (2019) 75:100–10. doi: 10.1016/j.eururo.2018.10.010
61. Subramanian I, Verma S, Kumar S, Jere A, Anamika K. Multi-omics data integration, interpretation, and its application. Bioinform Biol Insights (2020) 14:1177932219899051. doi: 10.1177/1177932219899051
62. Chowdhury N, Drake CG. Kidney cancer: An overview of current therapeutic approaches. Urol Clin North Am (2020) 47:419–31. doi: 10.1016/j.ucl.2020.07.009
63. Nargund AM, Pham CG, Dong Y, Wang PI, Osmangeyoglu HU, Xie Y, et al. The SWI/SNF protein PBRM1 restrains VHL-Loss-Driven clear cell renal cell carcinoma. Cell Rep (2017) 18:2893–906. doi: 10.1016/j.celrep.2017.02.074
Keywords: kidney renal papillary cell carcinoma, prognosis, immune microenvironment, drug response, multi-omics
Citation: Wang B, Li M and Li R (2023) Identification and verification of prognostic cancer subtype based on multi-omics analysis for kidney renal papillary cell carcinoma. Front. Oncol. 13:1169395. doi: 10.3389/fonc.2023.1169395
Received: 19 February 2023; Accepted: 17 March 2023;
Published: 05 April 2023.
Edited by:
Nan Zhang, Harbin Medical University, ChinaReviewed by:
Aimin Jiang, Second Military Medical University, ChinaJinxiao Li, Huazhong University of Science and Technology, China
Copyright © 2023 Wang, Li and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Rongshan Li, rongshanli13@163.com