Correlating Transcriptional Networks to Papillary Renal Cell Carcinoma Survival: A Large-Scale Coexpression Analysis and Clinical Validation.

We aimed to investigate the potential mechanisms of progression and identify novel prognosis-related biomarkers for papillary renal cell carcinoma (PRCC) patients. The related data were derived from The Cancer Genome Atlas (TCGA) and then analyzed by weighted gene coexpression network analysis (WGCNA). The correlation between each module and the clinical traits were analyzed by Pearson's correlation analysis. Pathway analysis was conducted to reveal potential mechanisms. Hub genes within each module were screened by intramodule analysis, and visualized by Cytoscape software. Furthermore, important hub genes were validated in an external dataset and clinical samples. A total of 5,839 differentially expressed genes were identified. By using WGCNA, we identified 21 coregulatory gene clusters based on 289 PRCC samples. We found many modules were significantly associated with clinicopathological characteristics. The gray, pink, light yellow, and salmon modules served as prognosis indicators for PRCC patients. Pathway enrichment analyses found that the hub genes were significantly enriched in the cancer-related pathways. With the external Gene Expression Omnibus (GEO) validation dataset, we found that PCDH12, GPR4, and KIF18A in the pink and yellow modules were continually associated with the survival status of PRCC, and their expressions were positively correlated with pathological grade. Notably, we randomly chose PCDH12 for validation, and the results suggested that the PRCC patients with higher pathological grades (II + III) mostly had higher PCDH12 protein expression levels compared with those patients in grade I. These validated hub genes play critical roles in the prognosis prediction of PRCC and serve as potential biomarkers for future personalized treatment.

INTRODUCTION

Renal cell carcinoma (RCC) is regarded as one of the most serious human diseases. In the US, there were an estimated 73,820 (∼4.18%) new cancer cases and 14,770 (∼2.43%) premature deaths in 20191. The disease encompasses several major subtypes, including clear cell RCC (ccRCC; 70%–80%), papillary RCC (PRCC; 10%–15%), chromophobe RCC (chRCC; 3%–5%), and renal oncocytoma (RO). Each subtype has different genetic characteristics and histological features2–6. PRCC is characterized by the presence of a fibrous vascular core of tumor cells arranged in a papillary pattern7, and there is no specific treatment option available for PCRR patients. Few studies investigate the underlying mechanisms of PRCC, delaying the development of new drugs. Besides, for metastatic PRCC patients, they commonly have an unfavorable prognosis, particularly compared to patients with ccRCC8. Therefore, it is necessary to identify new molecular markers that can be used to predict disease stage and clinical outcome in patients with PRCC. This can help understand its pathological mechanisms and provide personalized treatment.

To date, many biomarkers for PRCC have been identified, including several markers that can predict treatment outcomes and clinical prognoses, such as VHL, VEGF, CAIX, and HIF1α/2α mutations9. The mutation rates of these genes are low, and we still lack sufficient gene markers to predict prognosis for PRCC patients10. The rapid breakthroughs in genome-wide sequencing technology provide new ideas for the study of various clinical problems and pathological mechanisms related to cancer11. The Cancer Genome Atlas (TCGA) can help to improve diagnostic methods and ultimately improve the survival and prognosis of cancer patients. TCGA has generated an extensive collection of clinical information and genetic sequencing data that allows the systematic analysis of the underlying molecular mechanisms involved in the various clinical characteristics associated with cancer, such as tumor grade, pathological stage, histological type, diagnosis, and prognosis12. Weighted gene coexpression network analysis (WGCNA) explores the relationship between different genomes and clinical features by constructing a free-scale gene coexpression network13. The WGCNA algorithm is widely used to screen biological processes and therapeutic targets related to cancer as well as specific biomarkers associated with complex diseases14–16. Similarly, WGCNA is typically used to identify critical genes that are significantly associated with clinical indicators of tumor progression, including the stage, grade, and metastasis of different types of tumors17,18.

In this study, by generating a coexpression network analysis, we identified differences in expression between different subclassification and functional modules related to the prognosis of patients with PRCC. The flowchart of our work is illustrated in Figure 1. Additionally, the hub gene analysis may have important clinical significance and allow us to identify critical genes in each module that can be used as prognostic biomarkers or therapeutic targets for PRCC patients.

Figure 1

Flow diagram of strategy for data preparation, processing, analysis, and validation used in this study.

MATERIALS AND METHODS

Data Acquirements

All PRCC patients were recruited from the TCGA project (https://xenabrowser.net/datapages/?cohort=TCGA%20Kidney%20Papillary%20Cell%20Carcinoma%20(KIRP)&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443), which is a cohort study of PRCC19. In the current study, a cohort of 289 PRCC patients and 32 normal controls with matched mRNA expression profiles were used along with their clinical data.

Identification of Differentially Expressed Genes (DEGs) and Gene Set Enrichment Analysis (GSEA) Analysis

An in-depth analysis of DEGs between 289 tumor specimens and 32 normal controls was performed using the R package DESeq2 v1.20.020. Significant differentially expressed genes were selected with both p adj < 0.05 and |log2 (fold change)| ≥ 1. In addition, GSEA analysis was conducted based on gseapy v0.9.3, which could also display the gene distribution differences between the PRCC patients and normal controls.

WGCNA and Pathway Enrichment

Gene coexpression network analysis is a systems biology method and was performed by the WGCNA package21. In the coexpression network, the DEGs are represented by nodes, and the connectivity between the genes is defined by the correlation of the gene expression patterns22. Analysis of the network topology was used to determine the soft threshold power. The soft threshold power was used to successively calculate the gene coexpression similarity and the adjacency. The adjacency was then transformed into a topological overlap matrix (TOM). The eigengene is defined using the singular value decomposition (SVD) of the module expression matrix23,24. X(I) = (χ (I)) represents the gene expression matrix of the Ith module. The module gene corresponds to the indexes i = 1, 2,..., n I, and the ID of microarray samples corresponds to the indexes l = 1, 2, ..., m. The SVD of X(I) is expressed by the following formula:

The columns of orthogonal matrices U is called a left singular vectors, and the columns of orthogonal matrices V is called a right singular vectors. Σ value is only in the main diagonal, and we call it as singular values, while other elements are all equal to zero.

More precisely, U (I) represents a n (I) × m matrix with orthonormal columns; V(I) represents a m × m orthogonal matrix, and Σ(I) represents the m × m diagonal matrix of the singular values {|d (I) |}.

The matrices V(I) and Σ(I) are denoted by

Assuming that the singular values |d (I) | are arranged in a nonascending order, we call the first column of V(I) as the module eigengene:

Since the direction of each singular vector is undefined, we determine the direction of the characteristic gene by constraining the direction of each characteristic gene to make it positively correlated with the average gene expression of each module gene24. Finally, using the TOM for hierarchical clustering, the dynamic tree cut algorithm was applied to the module screening, and then the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed on the gene module, after which the modules related to PRCC were characterized.

Survival Analysis and Hub Gene Validation

We used the “survival”25 and “survminer”26 R packages to calculate the correlation between each module and overall survival (OS) and recurrence-free survival (RFS), with a value of p < 0.05 considered as statistically significant. In addition, the connections between the modules and the clinicopathological features were determined by Pearson’s correlation analysis, and a value of p < 0.05 was treated as statistically significant. The gene/eigengene expression differences in different subgroups were determined by t-test (two groups) or one-way analysis of variance (ANOVA) analysis (more than two groups). A hub gene is one of a series of genes that determine the characteristics of a module, and it has the highest degree of connectivity in a gene module. Based on their importance, we chose the top 10 hub genes within the significant modules and used an additional GEO dataset (GSE2748) to validate their usage for the prediction of clinical characteristics, such as prognosis and tumor grade, and confirmed their role in cellular function in vitro and in prognosis prediction by immunohistochemistry (IHC). In addition, the detail steps of the IHC were demonstrated in our previously published work27. The t-test (unpaired) was used to determine the protein expression differences between low- and high-grade PRCC patients. The PCDH12 antibody (catalog # PA5-20703) was purchased from Invitrogen (Carlsbad, CA, USA). And the IHC assay was performed on a tissue array containing 37 PRCC cases, which was bought from Xi’an AiDi Biotechnology Ltd. Co. (Xi’an, Shanxi Province, China). All patients gave written informed consent for the use of the tissue samples for research purposes. For the patients enrolled in the TCGA project, they have announced the ethics and policies on their own website (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga/history/policies). The research contents and research programs were reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Anhui Medical University (anyiyifuyuanlunshen-kuai-PJ-2019-09-11).

RESULTS

DEGs Screening and GSEA Analyses

The expression matrices for 289 PRCC samples and 32 normal controls were obtained from the TCGA database. Based on thresholds defined by |log2(fold change)| ≥ 1 and a p adj < 0.05, 5,839 DEGs were screened in the PRCC and normal control groups, of which 2,108 were upregulated and 3,371 were downregulated (Fig. 2A and B). Table S1 contains all the differential genes used in Figure 2B (see supplemental Table S1, available at https://pan.baidu.com/s/1xRGQBjHNJrQnwpWfrEQUPQ). We performed KEGG enrichment for these DEGs and found that these genes were enriched mostly within the cytokine–cytokine receptor interaction, neuroactive ligand–receptor interaction, and alcoholism pathways (p adj < 0.05) (Fig. 2C). These DEGs were then selected for subsequent analysis. In addition, we also conducted an online GSEA analysis, and the results suggested that genes in the apical surface, bile acid metabolism, fatty acid metabolism, KRAS signaling, and myogenesis pathways were significantly differentially expressed between PRCC cancer patients and normal controls (p adj < 0.05) (Fig. 2D–H).

Figure 2

Differentially expressed genes and pathway enrichment analyses. (A) Principal component analysis (PCA) for the papillary renal cell carcinoma (PRCC) tumor and normal tissues. (B) The differentially expressed genes between PRCC tumor and normal tissues (p < 0.05). (C) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of the differentially expressed genes (p < 0.05). (D–H) Gene Set Enrichment Analysis (GSEA) analyses between PRCC tumor and normal tissues (p < 0.05).

Coexpression Network Construction

WGCNA was performed on 5,839 DEGs from 289 samples (p adj < 0.05). The intergene connectivity in the gene network attained a scale-free network distribution (Fig. 3A). Therefore, we calculated the eigengene for entire modules and clustered them based on their correlation to further quantify their coexpression similarity. After merging two similar modules, a total of 21 modules were obtained (eigengene of the module, height > 0.2) (Fig. 3B).

Figure 3

Visualization of the gene network using a Heatmap plot. (A) The heatmap depicts the topological overlap matrix (TOM) showing all genes used in the analysis. The light red color represents a low overlap, and the progressively darker red color represents higher overlap. (B) Hierarchical eigengene diagram of all samples.

Finding Modules of Clinical Significance and Identifying the Hub Genes of Modules

The identification of the gene modules most significantly associated with different clinical characteristics has important biological implications, and the correlation between each module and the clinical features was observed in the module characteristic diagram (Fig. 4A), in which the largest correlation coefficient (r 2 = 0.5, p = 2 × 10−18) was observed between the pink module and the initial weights of the samples. In addition, we also observed that several modules were significantly associated with smoking (number of years) (negative: light cyan, purple, black, magenta, and light green modules; positive: green module) and new tumor events (negative: pink, salmon, midnight blue, and yellow modules) (p < 0.05). It can be inferred that these modules were of clinical significance. Interestingly, we found several modules significantly associated gender differences, including the tan, royal blue, blue, dark green, pink, gray, midnight blue, dark red, yellow, yellow-green, and light green modules (p < 0.05). Current work also analyzed the connection between sex and module eigengene expression, and the results are presented in Figure 4B. We then analyzed the correlations between these modules and clinicopathological features and found that the pink and yellow modules were significantly associated with the PRCC stage (p < 0.05) (Fig. 4C) and could serve as markers for predicting the tumor stage. In addition, we also found that the light cyan, light yellow, pink, gray, salmon, dark red, yellow, yellow-green, and light green modules were associated with new tumor events (p < 0.05) (Fig. 4D), which is partly consistent with the module characteristic diagram.

Figure 4

Correlation between module eigengene expression and clinicopathological features. (A) Module–trait relationships. (B) Module eigengene expression difference between males and females. (C) Module eigengene expression difference between different tumor stages. (D) Correlation between module eigengene expression and recurrent status. A value of p < 0.05 was regarded as statistically significant.

Prognostic Analyses and Network Construction

Pathological staging is an important factor that determines the prognosis of patients with PRCC. The higher the pathological staging, the worse the prognosis. The molecular biomarkers we screened for their correlation with pathological staging are also prognostic genes, in theory. In our work, we found that the gray, pink, light yellow, yellow, and salmon modules were significantly correlated with OS and RFS in PRCC (p < 0.05) (Fig. 5). When this was combined with the results shown in Figures 4 and 5, we found that the pink and yellow modules had the highest clinical significance. Therefore, we used Cytoscape software to visualize these two modules (Fig. 6). Consistent with the hub gene analysis results, we found that NOVA2, GPR4, MYCT1, and NOTCH4, which were the hub genes in the yellow module, were located at the center of the network and were significantly associated with other genes. Similar results were obtained for the pink module.

Figure 5

Survival analyses. (A–E) Correlation between the module eigengene expression and recurrence-free survival of PRCC patients. (F–J) Correlation between the module eigengene expression and overall survival of PRCC patients. A value of p < 0.05 was regarded as statistically significant.

Figure 6

Network analyses for the hub genes. (A) Network analysis for the hub genes in yellow module. (B) Network analysis for the hub genes in the pink module.

Pathway Enrichment and Independent Dataset Validation

By pathway enrichment analysis, we identified the significant KEGG pathways and GO terms that were enriched in the significant modules when the Benjamin-adjusted p value threshold was <0.05. As shown in Figure 7, for the gray 60 module, the module genes that were involved in the systemic lupus erythematosus, maturity onset diabetes of the young, alcoholism, chromosome maintenance, DNA methylation, signaling by nuclear receptors, and epigenetic regulation of gene expression event pathways were significantly enriched (p adj < 0.05).

Figure 7

KEGG analyses for the hub genes. The KEGG analyses for the hub genes in thw gray 60 module. (A) KEGG analysis. (B) GO_cellular component analysis. (C) Reactome analysis. A value of p < 0.05 was regarded as statistically significant.

In addition, based on the importance of the pink and yellow modules, we chose the top 10 genes (p < 0.05, excluding those genes that were not associated with OS or RFS) and used additional GEO datasets to further validate their importance. As shown in Figure 8, we found that PCDH12, GPR4, and KIF18A were continually associated with survival status in PRCC (p < 0.05) (Fig. 8A–C). Consistently, we also found that KIF18A, GPR4, and PCDH12 were significantly positively associated with histological classification type 2B (p < 0.05) (Fig. 8D), molecular classification class 2 (Fig. 8E), and grades 3 and 4 and were more highly expressed in tumor node metastasis stages 3 and 4 (p < 0.05) (Fig. 8G) and the M1 stage (p < 0.05, KIF18A and GPR4 only) (Fig. 8H). Based on the importance of these three genes, further IHC assay was performed on the PRCC clinical samples (Table 1). Moreover, we randomly chose PCDH12 for IHC validation, and the results suggested that the PRCC patients with higher pathological grades (II + III) mostly had high PCDH12 protein expression (p < 0.05), a result consistent with our above findings (Fig. 9).

Figure 8

External validation of the key hub genes within yellow and pink modules. (A–C) Survival analyses for the key hub genes in the yellow and pink modules. (D) Correlation between hub gene expression and histologic classification. (E) Correlation between hub gene expression and molecular classification. (F) Correlation between hub gene expression and pathological grade. (G) Correlation between hub gene expression and tumor–node–metastasis stage. (H) Correlation between hub gene expression and metastasis stage. A value of p < 0.05 was regarded as statistically significant.

Table 1

Clinical Characteristics of 37 Papillary Renal Cell Carcinoma (PRCC) Patients

Characteristic	n (%)
Gender
Male	24 (37.84)
Female	13 (35.14)
Age (years)
≤55	23 (62.16)
>55	14 (37.84)
Subtype
I	10 (27.03)
II	27 (72.97)
Tumor grade
I	30 (81.08)
II + III	7 (18.92)
TNM stage
I + II	32 (86.49)
III + IV	5 (13.51)
pT stage
T1 + T2	35 (94.59)
T3 + T4	2 (5.41)
pN stage
N0	35 (94.59)
N1	2 (5.41)
pM stage
M0	37 (100)
M1	0 (0)

Figure 9

Immunohistochemistry analysis for the PCDH12 expression in PRCC. A value of p < 0.05 was regarded as statistically significant. ***p < 0.001.

DISCUSSION

It is extremely important to investigate the molecular markers used for the diagnosis and treatment of PRCC patients. However, there is a lack of sufficient clinical use for the molecular markers currently associated with histological grade in PRCC. In this study, we used WGCNA to investigate PRCC-related coexpression modules and their clinical features. Significant modules, including the gray, pink, light yellow, yellow, and salmon modules, were shown to be correlated with OS and RFS in PRCC. The gray module was mainly related to pathways involved in alcoholism, necroptosis, and DNA complex packaging. Previous studies have reported that ethanol, which is the main component of alcoholic beverages, is metabolized into acetaldehyde, which is classified as carcinogenic to humans28. Alcohol consumption has now been identified as a risk factor for human cancer, but there is still a lack of direct evidence for PRCC. In addition, recent studies have shown that alcohol increases the number of cancer stem cells (CSCs) in cancer, which may be the basis of alcohol-induced tumors29–31. In addition, necroptosis involves the caspase-independent regulation of cell death that is executed by receptor-interacting protein 1 (RIP1), receptor-interacting protein 3 (RIP3), and mixed lineage kinase domain-like protein (MLKL)32,33. In recent years, tumor treatment based on tumor necrosis has been considered as a new antitumor treatment strategy34,35.

The pink module was mainly involved in the CDC6 association with ORC (origin complex), hepatitis C, and transcriptional activator activity, and it was observed to be negatively associated with age at initial diagnosis and new tumor events. Currently, sporadic evidence has shown a connection between the CDC6-related pathway and PRCC. CDC6 is a regulator of the initiation of DNA replication in human cells36,37. CDC6 may act as an ATPase switch that binds to Mcm2-7 in CDT1:CDC6:ORC. Human CDC6 protein levels decrease early in the G1 phase but remain unchanged throughout the cell cycle36,38. It has been reported that after cells enter the S phase, CDC6 is phosphorylated and excluded from the nucleus, which allows it to be easily ubiquitinated and degraded. Supplementation of CDC6 protein levels during G1 appears to be regulated by the E2F transcription factors38. In addition, the aberrant expression of genes involved in the regulation of transcriptional activator activity pathways could also influence tumorigenesis.

The yellow module was shown to involve genes mainly related to embryonic organ development, epithelial cell proliferation, renal system development, extracellular matrix, transcriptional activator activity, neuroactive ligand–receptor interaction, breast cancer, GPCR ligand binding, class A/1 signaling (rhodopsin-like receptors), and G alpha (q) signaling and was negatively associated with new tumor events and positively associated with the number of years of smoking. The yellow module-related pathways have been frequently reported to be involved in the tumor progression process. For example, the G protein coupled receptor (GPCR) superfamily is at the center of many different signaling pathways involved in various aspects of human physiology, and they are composed of approximately 800 different members39. Although they share a common seven-transmembrane (7TM) architecture40, they can recognize a variety of different ligands, including small molecules, peptides, lipids, and proteins, and may serve as potential therapeutic targets for cancer patients41–43.

Although the light yellow module could serve as a predictor for OS and RFS, the genes in this module may not be significantly enriched. Interestingly, the GO-Reactome results showed that genes in the light yellow module may be enriched in the inactivation, recovery, and regulation of the phototransduction cascade and the mine ligand-binding receptor pathways, and further functional studies are warranted to determine whether these pathways could influence PRCC initiation and progression. The salmon module was mainly enriched in the olfactory signaling and olfactory transduction pathways. We could only weakly connect the olfactory signaling or olfactory transduction pathways with cancers. Currently, emerging evidence suggests that the regulation of cell–cell recognition, exocytosis, migration, proliferation, apoptotic cycles, and pathfinding involve olfactory receptors and the related pathways44,45. In addition, there is growing evidence that olfactory receptors have great potential as diagnostic and therapeutic tools46. We also validated the stability of the hub genes using an external GEO dataset and found that most of the genes we identified were consistent with our findings. Furthermore, we employed the IHC assay and identified that the expression of PCDH12 was upregulated in PRCC patients with high pathological grade (II + III vs. I). However, until now, few studies have reported the function of PCDH12 in cancers. Future efforts are warranted to explore it.

In summary, using WGCNA and other methods combined with clinical data from PRCC patients to analyze the RNA-seq data, a set of biomarkers that can predict the prognosis of PRCC patients were identified. These results have important clinical implications that will contribute to personalized treatment.