Identification of co‑expression modules and hub genes of retinoblastoma via co‑expression analysis and protein‑protein interaction networks.

Retinoblastoma is a common intraocular malignant tumor in children. However, the molecular and genetic mechanisms of retinoblastoma remain unclear. The gene expression dataset GSE110811 was retrieved from Gene Expression Omnibus. After preprocessing, coexpression modules were constructed by weighted gene coexpression network analysis (WGCNA), and modules associated with clinical traits were identified. In addition, functional enrichment analysis was performed for genes in the indicated modules, and protein‑protein interaction (PPI) networks and subnetworks were constructed based on these genes. Eight coexpression modules were constructed through WGCNA. Of these, the yellow module had the highest association with severity and age (r=0.82 and P=3e‑07; r=0.72 and P=3e‑05). The turquoise module had the highest association with months (r=‑0.63 and P=5e‑04). The genes in the two modules participate in multiple pathways of retinoblastoma, and by combining the PPI network and subnetworks; 10 hub genes were identified in the two modules. The present study identified coexpression modules and hub genes associated with clinical traits of retinoblastoma, providing novel insight into retinoblastoma progression.

Introduction

Retinoblastoma is a retinal embryo malignancy that occurs in childhood and is triggered by mutations in the retinoblastoma gene (RB1) in cells of the developing retina (1,2). Retinoblastoma accounts for 3% of all cancer types in children, and is the most common intraocular malignancy in children (3). Moreover, children with a family history of retinoblastoma have a significantly increased risk of developing retinoblastoma (4). When diagnosed at an early stage and treated with standard regimens of systemic chemotherapy and topical consolidation therapy it is possible to successfully treat retinoblastoma, with most patients maintaining normal vision in at least one eye (5). However, retinoblastoma can cause blindness or death if not treated at an early stage (6). In addition, the use of standard chemotherapy regimens may result in notable toxicities and the development of secondary malignancies, such as soft tissue sarcomas, brain tumors and osteosarcoma (7–9). Therefore, there is an urgent need to understand the molecular characteristics of RB1 expression in order to improve the quality of life of patients with retinoblastoma (4,10).

It has been confirmed that clinical risk assessment in retinoblastoma requires combining clinical traits. For example, tumor staging, pathological grade and tumor laterality show associations with the patient's age at diagnosis, overall survival and second malignancy in retinoblastoma (11). High-risk histopathologic features (HRPFs) are closely related to poor prognosis, including tumor invasion of the optic nerve, choroid or anterior chamber, suggesting the need for postenucleation adjuvant chemotherapy (12). Therefore, it is necessary to identify gene coexpression modules associated with these clinical traits.

Network-based analysis has been used to characterize clinically significant genes and constitutes a means of presenting a variety of biological data, such as protein-protein interactions, gene regulation, cellular pathways and signal transduction. The core elements of biological networks may be identified by measuring nodes based on their network features. In addition, weighted gene coexpression network analysis (WGCNA) has been widely applied to develop gene coexpression modules of various diseases. For example, Zhou et al (13) successfully applied coexpression analysis to identify gene coexpression modules associated with the prognosis of pancreatic carcinoma. However, there are few studies on the gene coexpression modules associated with the clinical traits of retinoblastoma using WGCNA.

In this study in order to comprehensively analyze the molecular characteristics of RB1 expression, the aim was to identify coexpression modules associated with the clinical traits of retinoblastoma based on WGCNA and pathways involving the genes in these modules, as well as identifying hub genes related to clinical traits by combining the findings with protein-protein interaction (PPI) network. This study will contribute to an improved understanding of the development of retinoblastoma.

Materials and methods

Raw data and pretreatment

The GSE110811 dataset was downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). The dataset contains the gene expression profiles of 28 retinoblastoma samples based on the GPL16686 platform (Human Gene 2.0 ST array; Affymetrix; Thermo Fisher Scientific, Inc.) (12). Annotation information of microarray data was used to match probes with corresponding gene information. Probes matching with more than one gene were excluded, and average expression values were calculated for genes matching with more than one probe. Corresponding clinical information was also obtained. The workflow chart is shown in Fig. 1.

Figure 1.

Workflow chart for the construction of the gene coexpression network.

Construction of the gene coexpression network by WGCNA

A gene coexpression network was built with the WGCNA algorithm (14). Before constructing the network, the number of genes with different expression threshold was calculated, and the appropriate soft threshold value was then determined. Cluster analysis of 28 samples was performed at the appropriate threshold value using the flashclust package in R. The function pickSoftThreshold was applied to the calculate scale-free topology model fitting index (R2) and mean connectivity values corresponding to different soft threshold values (β; ranging from 1–20). The connectivity suggests the degree of correlation to which a gene is related to other genes in the network. The appropriate power value was estimated when R2 was 0.85. Furthermore, the gene expression matrix was transformed into an adjacency matrix and a topological overlap matrix (TOM). The corresponding dissimilarity of TOM (1-TOM) was calculated. The topological overlap is a measure of gene biological similarity based on pairwise gene coexpression correlation. A cluster dendrogram was constructed based on 1-TOM, and the DynamicTreeCut package (http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/BranchCutting/) was used to assign branches and generate modules. The heatmap package was employed to visualize the TOM among 400 randomly selected genes. For module construction, the default for maxBlockSize was 5,000, and deepSplit was set to 2. After calculating module eigengenes (MEs) using the function moduleEigengenes, modules with similar expression patterns were merged by the function mergeCloseModules. During this process, MEs are the major components of the first principal component in a module with the same expression profile, which can reflect the entire features of module genes.

Module-trait relationship construction

Correlation between MEs and clinical traits was analyzed by Pearson correlation analysis. R>0.8 was considered to indicate a strong correlation. For each expression profile, gene significance (GS) was calculated as the absolute value of the correlation between the expression profile and each clinical trait. Module membership (MM) was defined as the correlation between the expression profile and each ME (15). Thus, genes with a high significance for clinical traits and MM were identified by two parameters.

Functional and pathway enrichment analyses

To explore the biological processes and pathways in which the genes indicated in the module participate, functional enrichment analysis was performed via the clusterProfiler package in R language, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (16). GO terms include ‘biological process’, ‘molecular function’ and ‘cellular component’. P<0.05 after correction was set as the cutoff criterion.

PPI network and core subnetwork construction

PPI network analysis is a powerful tool for understanding the biological responses occurring in retinoblastoma. In the PPI network, a protein is defined as a node, and the interaction between two nodes is defined as an edge. The size of a node represents a degree: The larger the node, the larger the degree (17). The thickness of an edge indicates a correlation: The thicker the edge, the higher the correlation. Genes in the indicated modules were entered into the STRING online database (http://string-db.org/), and the cutoff value was set to 0.17. Related pairs of genes were retrieved and visualized with Cytoscape (version 3.4.0) (18). As an open platform with numerous plugins, the function of Cytoscape is to expand visualization options and network analysis power. Cytoscape cytoHubba plugin was used to select the top ten hub nodes according to maximal clique centrality (MCC) (19), and the Cytoscape MCODE plugin was applied to identify core subnetworks with highly interconnected nodes (20).

Results

Gene expression data

The box plot shows the expression value of each gene in 28 retinoblastoma samples (Fig. 2A). The average expression value of all mRNAs in each sample was nearly consistent, suggesting that the GSE110811 dataset was well normalized. As shown in Fig. 2B, two outlier samples (GSM3017126 and GSM3017124) were removed.

Figure 2.

Gene expression data for 28 retinoblastoma samples. (A) Box plots of normalized gene expression data from the GSE110811 dataset. (B) Sample clustering was performed to detect outliers.

Coexpression module construction

When the soft threshold power value β was equal to eight, the connectivity between genes conformed to a scale-free network distribution (Fig. 3A-D). Different coexpression modules of retinoblastoma were identified by hierarchical clustering and dynamic branch cutting, which are shown as different colors in Fig. 3E. After merging similar modules, eight coexpression modules were generated. The number of genes in the eight-coexpression modules is listed in Table I. The gray module represents the set of genes that were not assigned to any module. Eigengenes were used as representative gene expression profiles, quantifying module similarity through eigengene correlation (Fig. 3F).

Figure 3.

Coexpression module construction. (A and B) Effect of different soft threshold power values on the scale independence and mean connectivity degree. (C and D) Histogram of connectivity distribution and the scale free topology when β=0.85. (E) Clustering dendrogram of genes with dissimilarity based on topological overlap, together with assigned module colors. (F) Network heatmap plot showing the topological overlap matrix among 400 randomly selected genes. Light color represents low overlap; dark color represents high overlap.

Table I.

Number of genes in the eight-coexpression modules.

Module colors	Frequency
Black	44
Blue	836
Brown	146
Green	91
Grey	697
Red	49
Turquoise	3,039
Yellow	98

To explore meaningful modules associated with clinical traits, correlations between MEs and the following clinical traits were analyzed: Age; HRPF; metastasis; mild; moderate; severe; unilateral; survival time (months); mutation germline (mutation G); mutation not known or test not done (mutation ND); mutation none; and mutation nongermline (mutation NG; Fig. 4). Among them, HRPF has been confirmed to be associated with a poor prognosis in retinoblastoma. Moreover, anaplastic grades (mild, moderate and severe) are able to predict risk stratification for retinoblastoma. As shown in Fig. 4, different colors represent different correlation coefficients; furthermore, the green suggests negative correlation and the red stands for positive correlation. It was found that the yellow module had the highest association with severe anaplasia and age (r=0.82 and P=3e-07; r=0.72 and P=3e-05) and that the turquoise module had the highest association with months (r=−0.63 and P=5e-04). The eigengene dendrogram and heatmap showed that the MEyellow module was highly related to severe anaplasia and age, and the MEturquoise module was highly related to months (Fig. 5A-J). Therefore, the two modules were selected as meaningful modules associated with clinical traits. The scatterplots of GS vs. MM in the MEyellow and MEturquoise modules are shown in Figs. 6A-F and 7A-G, respectively. The above results revealed that the yellow and turquoise modules were significantly associated with clinical traits.

Figure 4.

Module-trait relationship construction. Each row corresponds to an ME, and each column corresponds to a trait. Each cell contains the corresponding correlation and P-value. The table is color-coded by correlation based on the color legend. Anaplastic grades include mild, moderate and severe. Mutation_G, RB1 mutation germline; mutation_ND, not known or test not done; mutation_NG, RB1 mutation nongermline; mutation_None, no RB1 mutation identified; months, survival time in months; HRPF, high-risk histopathologic feature; ME, module eigengene.

Figure 5.

Eigengene dendrogram and heatmap of the identified eight modules. (A) Eigengene dendrogram and heatmap of eight modules and age. (B) Eigengene dendrogram and heatmap of eight modules and HRPF. (C) Eigengene dendrogram and heatmap of eight modules and metastasis. (D) Eigengene dendrogram and heatmap of eight modules and mild. Eigengene dendrogram and heatmap of the identified eight modules. (E) Eigengene dendrogram and heatmap of eight modules and moderate. (F) Eigengene dendrogram and heatmap of eight modules and months. (G) Eigengene dendrogram and heatmap of eight modules and mutation_ND. (H) Eigengene dendrogram and heatmap of eight modules and mutation_NG. Eigengene dendrogram and heatmap of the identified eight modules. (I) Eigengene dendrogram and heatmap of eight modules and severe. (J) Eigengene dendrogram and heatmap of eight modules and unilateral. The color changes from blue to red, and the correlation gradually increases. HRPF, high-risk histopathologic feature; mutation_ND, not known or test not done; mutation_NG, RB1 mutation nongermline.

Figure 6.

Scatterplot of module membership vs. gene significance in the yellow module (A) Correlation between module membership in yellow module and gene significance for age. (B) Correlation between module membership in yellow module and gene significance for severe. (C) Correlation between module membership in yellow module and gene significance for mild. (D) Correlation between module membership in yellow module and gene significance for moderate. (E) Correlation between module membership in yellow module and gene significance for months. (F) Correlation between module membership in yellow module and gene significance for mutation_NG. mutation_NG, RB1 mutation nongermline.

Figure 7.

Scatterplot of module membership vs. gene significance in the turquoise module. (A) Correlation between module membership in turquoise module and gene significance for moderate. (B) Correlation between module membership in turquoise module and gene significance for months. (C) Correlation between module membership in turquoise module and gene significance for mutation_ND. (D) Correlation between module membership in turquoise module and gene significance for mutation_NG. (E) Correlation between module membership in turquoise module and gene significance for age. (F) Correlation between module membership in turquoise module and gene significance for severe. (G) Correlation between module membership in turquoise module and gene significance for unilateral. mutation_ND, not known or test not done; mutation_NG, RB1 mutation nongermline.

Functional enrichment analysis of genes in the yellow and turquoise modules

KEGG and GO enrichment analyses were performed for the genes in the yellow and turquoise modules, and a significant difference in the biological processes enriched by the genes in the two-coexpression modules was observed. For GO analysis, genes in the yellow module are mainly enriched in negative regulation of neuron differentiation, axon guidance, neuron projection guidance, negative regulation of neurogenesis and negative regulation of nervous system development (Fig. 8A). However, no KEGG pathways enriched with genes in the yellow module were found. For GO analysis, genes in the turquoise module are mainly enriched in chromosome segregation, organelle fission, chromosomal region, centrosome, tubulin binding, ATPase activity and coupling (Fig. 8B). In KEGG enrichment analysis, genes in the turquoise module are enriched in protein processing in the endoplasmic reticulum, cell cycle, proteasome, DNA replication, mismatch repair and protein export (Fig. 8C).

Figure 8.

Functional enrichment analysis of genes in the yellow and turquoise modules. (A) GO analysis of genes in the yellow module; (B) GO analysis of genes in the turquoise module; (C) KEGG enrichment analysis of genes in the turquoise module. The color changes from red to green, and the P-value gradually increases. The sizes of the circles represent the count of enriched genes.

PPI network construction

Genes in the yellow modules were used to construct the PPI network (Fig. 9A); subnetworks were then constructed to explore more specific and detailed information in the PPI network, as identified using the Cytoscape MCODE plugin. Two subnetworks (scores=4.6 and 3) were obtained for the yellow module (Fig. 9B and C). Genes in the turquoise modules were used to construct the PPI network (Fig. 10A). In total, three subnetworks (scores=4, 3.714 and 3.111) were obtained for the turquoise module (Fig. 10B-D). The Cytoscape CytoHubba plugin was employed to find the top ten hub genes with a high degree of connectivity between the nodes, according to the MCC. The top 10 hub genes identified from the PPI network of genes in the yellow and turquoise modules are shown in Fig. 11A and B, respectively, and the top ten hub genes in the yellow and turquoise modules are listed in Tables II and III, respectively.

Figure 9.

PPI networks and subnetworks of genes in the yellow module. (A) PPI network of genes in the yellow module; (B and C) subnetworks identified from the PPI network of genes in the yellow module (scores=4.6 and 3, respectively).

Figure 10.

PPI networks and subnetworks of genes in the turquoise module. (A) PPI network of genes in the turquoise module; (B-D) subnetworks identified from the PPI network of genes in the turquoise module (scores=4, 3.714 and 3.111, respectively).

Figure 11.

Top ten hub genes identified from the PPI networks of genes in the yellow module and turquoise module. (A) Yellow and (B) turquoise modules.

Table II.

Top 10 hub genes identified from the PPI network of genes in the yellow module.

Rank	Gene name	Score
1	BCL11A	3
2	SSBP3-AS1	2
2	EBF3	2
4	KCNQ5-IT1	1
4	LOC101929633	1
4	CA2	1
4	ASAP2	1
4	ROBO1	1
4	INSM2	1
4	LGI2	1

BCL11A, BAF chromatin remodeling complex subunit BCL11A; SSBP3-AS1, SSBP3 antisense RNA 1; EBF3, EBF transcription factor 3; KCNQ5-IT1, KCNQ5 intronic transcript 1; LOC101929633, uncharacterized LOC101929633; CA2, carbonic anhydrase 2; ASAP2, ArfGAP with SH3 domain, ankyrin repeat and PH domain 2; ROBO1, roundabout guidance receptor 1; INSM2, INSM transcriptional repressor 2; LGI2, leucine rich repeat LGI family member 2.

Table III.

Top 10 hub genes identified from the PPI network of genes in the turquoise module.

Rank	Gene name	Score
1	NIPBL	35
2	LRPPRC	31
3	GPBP1	24
4	UBE3A	21
5	HERC1	8
5	RNF214	8
7	CCNH	7
8	DYRK1A	4
9	STAM	3
9	CCT2	3

NIPBL, NIPBL cohesion loading factor; LRPPRC, leucine rich pentatricopeptide repeat containing; GPBP1, GC-rich promoter binding protein 1; UBE3A, ubiquitin protein ligase E3A; HERC1, HECT and RLD domain containing E3 ubiquitin protein ligase family member 1; RNF214, ring finger protein 214; CCNH, cyclin H; DYRK1A, dual specificity tyrosine phosphorylation regulated kinase 1Al; STAM, signal transducing adaptor molecule; CCT2, chaperonin containing TCP1 subunit 2.

Discussion

The present study systematically analyzed gene expression for retinoblastoma using WGCNA and PPI networks. These findings identified gene coexpression modules and hub genes associated with clinical traits of retinoblastoma, providing novel insight into retinoblastoma progression.

In the present study, 8 coexpression modules were identified based on WGCNA. Among them, the yellow module exhibited the highest association with histopathologic grading and age. Histopathologic grading included mild, moderate and severe. Histopathologic grading was determined based on tumor size, growth pattern, level of differentiation, degree of apoptosis, grade of anaplasia, tumor seeding, extent of tissue invasion and the presence of retinocytoma. High parental age is associated with increased risk of sporadic hereditary retinoblastoma (21). Therefore, early detection is important in curing the disease using surgical treatment (22). It has been recorded that a delay of more than 6 months from the first clinical sign to diagnosis is associated with 70% mortality in developing countries (23). The turquoise module had the highest association with survival time (months), indicating that the module could become a potential prognostic factor. Therefore, the two modules were considered important in retinoblastoma. GO analysis revealed that the genes in the yellow module were mainly enriched in the negative regulation of neuron differentiation, axon guidance, neuron projection guidance, neurogenesis and development of the nervous system. Increasing evidence suggests that neuronal differentiation of neuroblastoma cell lines is induced by a number of factors, such as p73, Tropomyosin receptor kinase A and Ubiquitin C-Terminal Hydrolase L1 (24–26). According to GO enrichment analysis, the genes in the turquoise module are mainly enriched in chromosome segregation, organelle fission, chromosomal region, centrosome, tubulin binding, and ATPase activity and coupling. In KEGG enrichment analysis, it was shown that genes in the turquoise module were enriched in protein processing in the endoplasmic reticulum, cell cycle, proteasome, DNA replication, mismatch repair and protein export. Hub genes in the turquoise module were involved in these pathways. For example, among these pathways, hub gene CCNH was enriched in the cell cycle pathway. Overall, the genes in the two modules may play important roles in the development of retinoblastoma.

In this study, the top ten hub genes with a high degree of connectivity according to MCC in the yellow and turquoise modules were identified. The degree of a particular protein is related to the necessity of the gene, and proteins with higher numbers are more likely to be essential (27). As biological networks are heterogeneous, it is necessary to capture essential proteins using a variety of methods. The top 10 hub genes in the yellow module included: BAF chromatin remodeling complex subunit BCL11A (BCL11A); SSBP3 antisense RNA 1 (SSBP3-AS1); EBF transcription factor 3, KCNQ5 intronic transcript 1 (KCNQ5-IT1); uncharacterized LOC101929633 (LOC101929633); carbonic anhydrase 2 (CA2); ArfGAP with SH3 domain, ankyrin repeat and PH domain 2; roundabout guidance receptor 1 (ROBO1); INSM transcriptional repressor 2; and leucine rich repeat LGI family member 2 (LGI2). Those in the turquoise module included: NIPBL cohesion loading factor (NIPBL); leucine rich pentatricopeptide repeat containing (LRPPRC); GC-rich promoter binding protein 1; ubiquitin protein ligase E3A; HECT and RLD domain containing E3 ubiquitin protein ligase family member 1; ring finger protein 214; cyclin H (CCNH); dual specificity tyrosine phosphorylation regulated kinase 1Al (DYRK1A); signal transducing adaptor molecule; and chaperonin containing TCP1 subunit 2 (CCT2). However, the function of these genes in retinoblastoma remains unclear. A previous study found that BCL11A plays a critical role in several diseases as an oncogene (28). Consistent with the present study, it has been confirmed that in epigenetic complexes, the transcription factor BCL11A competes with histone H3 for binding to Retinoblastoma-binding proteins 4 and 7 (29). Therefore, BCL11A has a widespread role in the development of retinoblastoma (30). It has been reported that ROBO1 is localized to the cell membrane; in primary and metastatic prostate cancer its expression is lower (31). However, the role of ROBO1 in retinoblastoma has not yet been studied. The present results showed that ROBO1 in the yellow module may be associated with retinoblastoma, but this requires further validation. LGI2 is secreted by neurons and acts on members of the metalloproteinase-deficient ADAM neuronal receptor family, with roles in synaptic remodeling (32). Similarly, GO analysis in the present study indicated enrichment of the genes in the yellow module in synaptic guidance. Among the genes in the turquoise module, NIPBL plays a critical and regulatory role in cohesion loading on chromatin, as well as having roles in gene expression and transcriptional signaling (33,34). Zinc finger protein 609 may participate in the regulation of cortical neuron migration (35). Upregulated LRPPRC may promote tumourigenesis in various tumors, such as lung adenocarcinoma, esophageal squamous cell carcinoma, stomach, colon, mammary and endometrial adenocarcinoma, as well as lymphoma (36,37). Regardless, the present study is the first to report that LRPPRC may affect the development of retinoblastoma. Furthermore, CCNH, DYRK1A and CCT2 are involved in the regulation of the cell cycle, which is consistent with the present functional enrichment analysis results (38–42). Due to the unclear molecular mechanisms of retinoblastoma, the ability to treat the disease remains limited (43–45). Moreover, patients with hereditary retinoblastoma have an increased risk of subsequent malignant neoplasms after radiotherapy or chemotherapy (46,47). The present study revealed that these hub genes may be involved in the molecular mechanisms underlying the development and progression of retinoblastoma, though further validation is required.

Currently, WGCNA is widely used to analyze large-scale data sets and find modules for highly related genes. To our knowledge, the present study is the first to analyze the correlation between coexpression modules and clinical traits of retinoblastoma via WGCNA. In a previous study using WGCNA analysis four hub genes were identified as prognostic markers in another eye tumor, uveal melanoma (15) and another study identified four hub genes associated with the progression of uveal melanoma (48). However, identical genes associated with the two eye tumors, uveal melanoma and retinoblastoma, have not been found. The genes identified using WGCNA in uveal melanoma have been not used in a clinical setting; however, functional enrichment studies have not revealed the role of these genes. Since WGCNA has been widely used to construct coexpression modules, it is a useful method to study retinoblastoma.

However, the present study has several limitations that should be noted. First, this study was based on analysis of data extracted from a limited sample size; therefore, the genes identified cannot be generalized to the majority of patients with retinoblastoma. Second, the biological mechanisms of the hub genes, such as SSBP3-AS1, KCNQ5-IT1, LOC101929633 and CA2, are still unknown. Third, this study was limited by the absence of experimental evidence, so large-scale studies are needed for validation.

In the present study, two meaningful gene coexpression modules for retinoblastoma were identified based on WGCNA. Among the genes of the two-coexpression modules, hub genes were identified by PPI networks and subnetworks, and based on previous studies, these genes may play a critical role in retinoblastoma. Therefore, this research is helpful in understanding the molecular mechanisms of retinoblastoma.

In this study, 8 coexpression modules were constructed via WGCNA. The yellow module and turquoise module had highly significant associations with clinical traits. In addition, the genes in the two modules participate in multiple pathways in retinoblastoma. Ten hub genes in the two-coexpression modules were identified according to the PPI network. This study identified meaningful gene coexpression modules and hub genes associated with a number of clinical traits of retinoblastoma that may contribute to the development of the disease.