Identification and validation of a 44-gene expression signature for the classification of renal cell carcinomas

Background Renal cancers account for more than 3% of all adult malignancies and cause more than 23,400 deaths per year in China alone. The four most common types of kidney tumours include clear cell, papillary, chromophobe and benign oncocytoma. These histological subtypes vary in their clinical course and prognosis, and different clinical strategies have been developed for their management. Some kidney tumours can be very difficult to distinguish based on the pathological assessment of morphology and immunohistochemistry. Methods Six renal cell carcinoma microarray data sets, including 106 clear cell, 66 papillary, 42 chromophobe, 46 oncocytoma and 35 adjacent normal tissue samples, were subjected to integrative analysis. These data were combined and used as a training set for candidate gene expression signature identification. In addition, two independent cohorts of 1020 RNA-Seq samples from The Cancer Genome Atlas database and 129 qRT-PCR samples from Fudan University Shanghai Cancer Center (FUSCC) were analysed to validate the selected gene expression signature. Results A 44-gene expression signature derived from microarray analysis was strongly associated with the histological differentiation of renal tumours and could be used for tumour subtype classification. The signature performance was further validated in 1020 RNA-Seq samples and 129 qRT-PCR samples with overall accuracies of 93.4 and 93.0%, respectively. Conclusions A 44-gene expression signature that could accurately discriminate renal tumour subtypes was identified in this study. Our results may prompt further development of this gene expression signature into a molecular assay amenable to routine clinical practice.


Background
According to the newest Globocan 2012, renal cancers are the 17th most common malignancy, accounting for more than 3% of adult malignancies and causing approximately 23,400 deaths per year in China alone [1,2]. In 2011, the overall incidence of renal cancers in China rose to 3.35 cases per 10 5 people, and the estimated mortality rate was 1.12 deaths per 10 5 people [3]. According to the 2016 World Health Organization (WHO) classification, there are 16 subtypes of renal cell carcinoma (RCC), a family of carcinomas that arise from renal tubule epithelia [4]. Currently, the four most common types of kidney tumours include clear cell RCC (ccRCC), papillary RCC (pRCC), chromophobe RCC (chRCC) and benign oncocytoma [4]. These histological subtypes vary in their clinical course and outcomes, and different clinical management strategies have been developed for their treatment. Among patients with the four most common types, patients with ccRCC have the worst prognosis, and there are differences between the prognosis of patients with pRCC and chRCC [5]. Different genetic alterations induce the development of renal tubules into RCCs of varying histological subtypes that exhibit different gene expression patterns or mutations, thus providing specific molecular candidates for targeted therapy (e.g., mTOR, VEGF, KIT, and checkpoint inhibitors) [6]. Improving the molecular understanding of the mechanisms underlying RCC subtypes has facilitated the development of targeted therapies and biomarkers in response to treatment [6]. Distinguishing between some types of kidney tumours based on morphology and immunohistochemistry can be very difficult for pathologists, while the correct identification of these subtypes is important for making precise decisions regarding therapeutic regimens.
Recent studies focused on microarray profiling of different RCC subtypes to develop accurate diagnostic RCC biomarkers. Using microarray analysis of renal tumours, claudin-7 mRNA, a distal nephron marker, was overexpressed in chRCC compared with that in oncocytoma, ccRCC, and pRCC [7]. Further immunohistochemical analysis of two independent cohorts showed that claudin-7 expression was detected in 67 and 100% of chRCCs, 0 and 7% of ccRCCs, 28 and 90% of pRCCs, and 26 and 45% of oncocytomas [8,9]. These studies revealed the potential of claudin-7 as a biomarker for distinguishing chRCC from the remaining three RCC subtypes and indicated the accuracy of microarray technology for detecting diagnostic biomarkers. Compared with classifying diseases using a single gene marker, simultaneously quantifying the expression of numerous genes may potentially capture the complex physiopathology underlying tumourigenesis and the development of specific RCC subtypes. Several studies have used microarray technology to identify gene expression signatures for the classification of RCCs. Chen and coworkers published a four-gene panel that could classify RCC subtypes with an estimated prediction accuracy of 96% [10]. Youssef and colleagues also reported a classification system using miRNA signatures with a maximum of four steps that had sensitivities of 97% for distinguishing normal cells from RCC, 100% for the ccRCC subtype, 97% for the pRCC subtype, and 100% accuracy in distinguishing the oncocytoma subtype from the chRCC subtype [11].
In this study, to identify novel gene biomarkers for the classification of RCC subtypes, we performed an integrative analysis of six microarray data sets (n = 295). The selected genes in the training set were validated in 1020 RNA-sequencing samples from The Cancer Genome Atlas (TCGA) database and then tested in 129 independent specimens by qRT-PCR. A 44-gene signature was identified and validated as being highly sensitive and specific for the classification of RCCs.

Gene expression database curation
Gene expression data sets of 1315 renal tumours with histologically confirmed subtypes and adjacent normal tissues were collected from public data repositories (e.g., ArrayExpress, Gene Expression Omnibus (GEO), and TCGA data portal) and curated to form a comprehensive RCC transcriptome database. Array-based gene expression profiling of 295 tissue samples obtained from six GEO data sets (GSE12090, GSE15641, GSE19949, GSE8271, GSE7023 and GSE19982) was mainly conducted on two different Affymetrix oligonucleotide microarray platforms, GeneChip Human Genome U133A Array and U133Plus 2.0 Array. Detailed descriptions of the specimen characteristics and clinical features are provided in the original studies [12][13][14][15]. The sequence-based gene expression profiles of 1020 tissue samples (including 534 ccRCC, 291 pRCC, 66 chRCC and 129 normal kidney samples) were generated on an Illumina HiSeq 2000 RNA sequencing platform and retrieved from the cBioPortal for Cancer Genomics [16]. The gene expression profiles consisted of transcriptomic data for 20,500 unique genes, and clinical information for the selected samples was retrieved from the "Clinical Biotab" section of the data matrix based on the Biospecimen Core Resource IDs of the patients.

Microarray data processing and normalization
Gene expression data analysis was performed using R software and packages from the Bioconductor project [17][18][19]. We used the Single Channel Array Normalization (SCAN) approach from the SCAN-UPC package to process Affymetrix microarray data [20,21]. Upon normalising each raw CEL file, SCAN outputs probe-level expression values. We further used the custom mapping files from the BrainArray resource to summarise probe-level intensities directly to gene-level expression values [22]. Thus, probes mapping to multiple genes and other problems associated with older generations of Affymetrix probe designs were avoided. After normalization, we applied the ComBat approach to adjust for batch effects [23].

Gene signature identification and performance assessment
To identify a gene expression signature, we used the support vector machine-recursive feature elimination (SVM-RFE) algorithm for feature selection and classification modelling [24]. For multi-class classification, a one-versusall approach was used by which multiple binary classifiers were first derived for each subtype. The results are reported as the subtype classifying the test sample with the highest confidence. For each specimen, the predicted subtype was compared with the reference diagnosis, and a true positive result was indicated when the predicted subtype matched the reference diagnosis. When the predicted subtype and reference diagnosis did not match, the specimen was considered a false positive. For each subtype on the panel, sensitivity was defined as the ratio of true positive results to the total positive samples analysed, while specificity was defined as the ratio (1 -false positive)/(total tested -total positive).

Biological network and functional enrichment analysis
Enrichment analysis of Gene Ontology and molecular pathways was performed using the Lynx Systems Biology Tool [25]. All significance tests were two-sided, and a false discovery rate less than 0.05 was considered significant. Biological network analysis was performed with Networ-kAnalyst software [26,27]. Protein-protein interaction information was retrieved from the IMEx Interactome Database [28]. A dense network was connected by retaining only the seed proteins as well as minimum essential non-seed proteins to study the key interactions.

qRT-PCR analysis
We included 121 renal tumour samples and 8 nontumour kidney tissues for qRT-PCR analyses. Written informed consent was obtained from all participants. The study was approved by the Ethics Committee of Fudan University Shanghai Cancer Center (FUSCC), China. Of the 121 tumours, 26 were ccRCC, 40 were chRCC, 28 were pRCC, and 27 were oncocytoma. Total RNA was isolated from formalin-fixed paraffin-embedded (FFPE) tissue sections using a FFPE Total RNA Isolation Kit (Canhelp Genomics, Hangzhou, China). Briefly, the paraffin sections were placed in sterile 1.5-ml microcentrifuge tubes, deparaffinized with 100% xylene, and washed twice with 100% ethanol. The deparaffinized tissue was digested with proteinase K at 56°C for 15 min and then incubated at 80°C for another 15 min to partially reverse nucleic acid crosslinking. The samples were treated with DNase and eluted in 40 μl RNase-free water. The concentration of total RNA was spectrophotometrically determined using total absorbance at 260 nm, and the purity was quantified using the A260/A280 ratio. RNA samples with A260/A280 ratios of 1.9 ± 0.2 were included in this study.
For each sample, cDNA was generated from isolated total RNA using a High-Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (Applied Biosystems, Foster City, CA, Unites States). Primers and MGB probes for the tested gene candidates and control gene were designed using Primer Express software (Applied Biosystems). Subsequently, the expression level of gene candidates was analysed on an Applied Biosystems 7500 Real-Time PCR system using TaqMan Gene Expression Assays (Applied Biosystems). The PCR program was initiated at 95°C for 10 min, followed by 40 thermal cycles, each at 95°C for 15 s and at 60°C for 1 min.

Establishment of the RCC Transcriptome database
To create a RCC transcriptome database for subtype classification, we performed a systematic search of major biological data repositories (e.g., ArrayExpress, GEO, and TCGA) to collect gene expression data sets from ccRCC, pRCC, chRCC, oncocytoma and adjacent normal tissue samples. Overall, we accumulated the gene expression profiles of 1315 tissue samples to form a comprehensive RCC transcriptome database. To identify a reliable gene expression signature, we adopted a training-testing-validation approach in this study. First, the microarray-based gene expression profiles of 295 specimens were retrieved from the database and curated to form a training set. Second, two independent sets were used to test and validate the classification performance of the gene expression-based signature; one was composed of the sequence-based gene expression profiles of 1020 specimens (Test Set 1), and the other was composed of the gene expression profiles of 129 specimens that were analysed with qRT-PCR (Test Set 2). Figure 1 depicts the three distinct phases of our study design, and Table 1 summarises the clinical characteristics of the samples in the study.

Identification of a 44-gene signature in the training set
The training set consisted of 106 ccRCC, 66 pRCC, 42 chRCC, 46 oncocytoma and 35 adjacent normal tissue samples. After the data normalization and annotation steps, a matrix of 12,263 unique genes in 295 samples (≈ 3.5 million data points) was prepared for downstream bioinformatics analyses. Extracting a subset of informative genes from high-dimension genomic data is a critical step for gene expression signature identification. Although many algorithms have been developed, the SVM-RFE approach is considered one of the best gene selection algorithms. For each subtype, we used the SVM-RFE approach to (1) evaluate and rank the contributions of each gene to the optimal separation of a specific subtype from other subtypes; (2) select the top 10 ranked genes as the most differentially expressed for that subtype; (3) repeat the process for each subtype, and obtain 5 lists of the top 10 gene set. After removing redundant features, 44 unique genes (listed in Table 2) were obtained and used to cluster the 295 training set samples. The average linkage hierarchical clustering method was performed where the metric of similarity was Pearson's correlation between the 44-gene expression profiles of the samples. As shown in Fig. 2a, the samples were clustered into five groups that closely followed the histological subtypes. Among the four tumour subtypes, the oncocytoma and chRCC samples clustered together, whereas the ccRCC samples were more similar to pRCC samples.

Functional enrichment and biological network analysis
We further investigated whether the 44 candidate genes exhibited biological features relevant to renal carcinogenesis. As shown in Table 3, the most significantly enriched gene categories are involved in insulin-like growth factor binding, transmembrane transport of small molecules, cocaine, amphetamine addiction, etc. Interestingly, seven of the 44 candidate genes (ASS1, DEFB1, IGFBP6, LCN2, SERPINA5, UMOD and VCAN) were indeed overrepresented in the "Renal-cell cancer" gene set (p < 1.4 E-5). More specifically, AQP6, CLDN8 and KRT7 were overrepresented in the "Renal oncocytoma" gene set (p < 6.1 E-6). We also explored the underlying biological networks of these 44 candidate genes. We used the 44 genes as seeds to generate a minimum protein-protein interaction network. As shown in Fig. 3, the network includes 33 genes of the 44-gene set and is centred on essential nodes such as APP, ASS1, ATF2, CRYAB, HNF1A, S100A2 and UBC. Enrichment analysis revealed that the most significant molecular networks were the TGF beta signalling pathway, Androgen receptor signalling pathway, Transcriptional misregulation in cancer, etc. (Table 4).

Performance assessment with 5-fold cross-validation
As an initial step, we assessed the performance of the classifier using 5-fold cross-validation within the training set. In 5-fold cross-validation, we created the training and testing sets by splitting the data into five equally sized subsets. We treated a single subsample as the testing set and the remaining data as the training set. We then ran and tested models on all five datasets and averaged the estimates. Given the limited sample size of the training set, we repeated the 5-fold cross-validation process 1000 times and estimated the average classification accuracy and corresponding 95% confidence interval (95% CI). The 44-gene expression signature showed an overall accuracy of 95.7% (95% CI: 0.912 to 1.00) with notable variation between different subtypes. Sensitivities ranged from 88.0% (chRCC) to 98.1% (ccRCC). Using this internal validation of the training set, these data provided a preliminary estimate of classification performance.

Independent validation in renal Tumours profiled with next-generation sequencing
The final classification model of the 44-gene expression signature was established using the entire training set and then applied to an independent validation set comprising 534 ccRCC, 291 pRCC, 66 chRCC and 129 adjacent normal tissue specimens profiled with nextgeneration sequencing (Test Set 1). The hierarchical clustering of 44 genes and 1020 samples revealed distinct patterns between ccRCC, pRCC, chRCC and adjacent normal samples (Fig. 2b). With the 44-gene  The detailed sensitivities and specificities are listed in Table 5.

Clinical validation of the 44-gene signature by qRT-PCR analysis
Microarray and RNA-sequencing data provide a global assessment of transcriptomic variations, but their resolution and accuracy are limited in individual gene analyses, and they remain difficult to use in clinical practice. qRT-PCR is generally considered the "standard procedure" assay for measuring individual gene expression and often used to confirm the findings of microarray and RNA-sequencing analyses. Hence, we further evaluated the expression levels of 44 genes by qRT-PCR in an independent cohort of 121 RCC tumours (comprising 26 ccRCC, 28 pRCC, 40 chRCC, and 27 oncocytoma) and 8 normal kidney tissues (Test Set 2). Figure 2c shows the hierarchical clustering of the 44 genes and 129 samples based on the qRT-PCR data. As seen in the figure, distinct patterns were observed between four tumour subtypes and adjacent normal samples. With the 44-gene expression signature, 29 samples were classified as ccRCC, 25 as pRCC, 39 as chRCC, 26 as oncocytoma and 10 as normal kidney tissues. Overall, the gene expressionbased assignments reached 93.0% overall agreement with the reference diagnoses (120 of 129; 95% CI: 0.868 to 0.966). Sensitivities ranged from 89.3% (pRCC) to 100% (normal tissue), while specificities ranged from 96.1% (ccRCC) to 100% (chRCC and pRCC). The detailed sensitivities and specificities are listed in Table 5.

Discussion
Due to the comprehensive development of highthroughput microarray and next-generation sequencing technologies, as well as the comprehensive efforts of systematic cancer genomics projects, numerous genomic data sets were utilised in our research. In this study, we identified a 44-gene expression signature for the accurate and robust classification of RCC subtypes (ccRCC, pRCC, chRCC, and oncocytoma). The 44-gene expression signature demonstrated an overall accuracy of 95.7% for 4 RCC subtypes by cross-validation of the training set profiled with the high-throughput microarray and 93.4% in an independent test set of 1020 RCC and normal kidney samples profiled with next-generation sequencing. Furthermore, we tested the signature on an independent cohort by qRT-PCR. An overall accuracy of 93.0% was achieved with the 129 RCC samples with 4 subtypes and normal specimens. This signature may serve as a reliable diagnostic tool to aid pathologists with the growing unmet need for RCC classification. Kidney tumour subtypes are characterised by different genetic mutations and chromosomal variations and thus present different gene expression profiles. Numerous molecules have been reported as capable of distinguishing kidney tumour subtypes. For example, vascular cell adhesion molecule 1 (VCAM1) was reportedly significantly upregulated in ccRCC and pRCC, whereas it was downregulated in chRCC and oncocytoma [29]. Furthermore, positive immunoreactivity of the metastasis suppressor protein KAI1 was often detected in chRCC specimens and rarely in ccRCC and oncocytoma specimens [30], and GST-alpha mRNA expression was higher in most ccRCCs than in other kidney tumours [31]. However, in addition to being unable to consistently distinguish RCC subtypes based on regular microscopic morphology, single molecules seldom exhibit extensive power for classifying all 4 major renal tumour subtypes. Therefore, comprehensive analysis of multiple gene expression panels is necessary for the classification of renal tumour types.
Based on the expression patterns of 44 genes, we classified the 4 most common renal tumour subtypes, ccRCC, pRCC, chRCC, and oncocytoma, with sensitivities ranging from 88% (chRCC) to 98% (ccRCC) in the training set, 90.9% (chRCC) to 94.6% (normal tissue) in Test Set 1, and 89.3% (pRCC) to 100% (normal tissue) in Test Set 2. In addition, the diagnostic histological classification accuracy was higher than that obtained with any of the genes used alone. The chRCC and oncocytoma samples displayed almost identical gene expression profiles for MAL, TMEM255A, RHCG, ATP6V0A4, STAP1, and DEFB1, as demonstrated by both RNA microarray and RNA sequencing, which is in agreement with the known fact that chRCC and oncocytoma are related neoplasms [32]. However, because chRCC is potentially malignant, and oncocytoma appears to be a benign mimic of RCC [4,33], the potential subtle difference in gene expression is expected, and the distinction between both subtypes has important clinical significance. Thus, we proposed that biomarkers identified by gene expression profiles accumulated from large cohorts indeed help to discriminate important and difficult differential diagnoses.
Several studies have reported the promise of gene or protein expression-based signatures in the classification of RCC subtypes. Unlike many studies in which samples were often collected from single central or ethnic cohorts, our approach exploited tumour samples from two large databases; samples extracted from the GEO database were used for construction of the classification panel, and samples from the TCGA database were extracted for testing our 44-gene expression signature. In addition, we further validated our 44-gene expression signature in an independent Chinese cohort using qRT-PCR. In a clinical scenario, the application of multi-centre, multi-ethnic data would greatly increase the reliability and universal applicability of our 44-gene expression signature. In this study, we showed that the 44-gene expression signature could reliably identify the tumour subtypes in 95.7% of the 295 samples tested. This accuracy is comparable to that of other signatures established by mRNA or miRNA biomarkers (ranging from 90 to 96%) [10,11,34]. The performance of this mRNA signature analysis by qRT-PCR also compares favourably with protein signature analysis by immunohistochemistry, the current clinical practice standard, which has shown 78-87% accuracy in identifying RCC samples using AMACR, CK7, and CD10 [35]. Moreover, analysis of the expression patterns of 44 genes by qRT-PCR classified the 4 most common renal tumour subtypes with 100% sensitivity in distinguishing normal from RCC, 96.2% for the ccRCC subtype, 92.5% for the chRCC subtype, 89.3% for the pRCC subtype, and 92.6% for the oncocytoma subtype; this signature is also comparable to other signatures (97% in distinguishing normal from RCC, 98%-100% for the ccRCC subtype, 93% for the chRCC subtype, 97-98% for the pRCC subtype, and 86% for the oncocytoma subtype) [10,11,34]. In routine clinical settings, the most commonly used diagnostic materials are FFPE samples; thus, further research is needed to successfully translate the 44-gene signature from gene expression microarrays and qRT-PCR to immunohistochemistry, thus allowing widespread access and applications in clinical diagnoses.

Conclusion
In conclusion, in the present study, we developed and validated a 44-gene expression-based signature for the classification of RCC subtypes. Our results may prompt further development of this gene expression signature into a molecular assay amenable to routine clinical practice. We foresee its application in cases wherein morphology and immunohistochemistry fail to distinguish between renal tumour subtypes. Further studies are   needed to determine the role of our gene expressionbased signature in personalised therapy choices and the prognosis of therapeutic outcomes for RCC patients with different subtypes.