Network analysis based on TCGA reveals hub genes in colon cancer.

Colorectal cancer (CRC) is the third most widespread cancer in the world. Although many advances have been made in molecular biology, novel approaches are still required to reveal molecular mechanisms for the diagnosis and therapy of colon cancer. In this study, we aimed to determine and analyse the hub genes of CRC. First, we explored the mRNA and microRNA (miRNA) expression profiles of colon carcinoma, then we screened target genes of differentially expressed miRNAs and obtained the intersection between differently expressed genes and target genes. Gene Ontology (GO) classification and KEGG pathway analysis of differently expressed genes were performed, and gene-miRNA and TF-gene-miRNA networks were constructed to identify hub genes, miRNAs, and TFs. In total, 3436 significant differentially expressed genes (1709 upregulated and 1727 downregulated) and 216 differentially expressed miRNAs (99 upregulated and 117 downregulated) were identified in colon cancer. These differentially expressed genes were significantly enriched in GO terms and KEGG pathways, such as cell proliferation, cell adhesion, and cytokine-cytokine receptor interaction signalling pathways. GCNT4, EDN2, and so on were located in the central hub of the co-expression network. MYC, WT1, mir-34a, and LEF1 were located in the central hub of the network of TF-gene-miRNA. These findings increase our understanding of the molecular mechanisms of colon cancer and will aid in identifying potential targets for diagnostic and therapeutic usage.

Introduction

Colorectal cancer (CRC) is the third most widespread cancer in the world; over 1 million newly diagnosed cases and over half a million CRC-related deaths are estimated to occur each year [1, 2]. It is generally accepted that CRC mainly develops through two different genetic pathways: the chromosomal instability (CIN) pathway, which is responsible for approximately 70–80% of all CRC cases [3], and the microsatellite instability (MSI) pathway, which is related to inactivation of the mismatch repair (MMR) gene system, leading to inactivation of mutated tumour suppressor genes. Recently, a third pathway has been identified; this pathway is called the serrated pathway due to the serrated appearance of tumours upon histological analysis. This pathway involves hypermethylation of specific DNA regions near the promoter genes [4]. Although much progress has been made in molecular biology, novel approaches are still required to reveal the molecular mechanisms of CRC.

MicroRNAs (miRNAs) are short, 18–27-nucleotide, non-coding RNA sequences that function in post-transcriptional regulation by binding mRNA. Many studies have shown that miRNAs are involved in cell differentiation, proliferation, and apoptosis and act as tumour suppressor genes or proto-oncogenes [5]. Transcription factors (TFs) are modular proteins that regulate gene transcription by binding to the promoter region of target genes via their DNA-binding domain [6]. Both TFs and miRNAs are involved in the regulation of various biological processes, including cell proliferation, differentiation, and apoptosis.

With the rapid development of high-throughput technologies, such as RNA sequencing, gene expression analysis based on high-throughput data has become more important for exploring the biological mechanisms of disease. The Cancer Genome Atlas (TCGA) provides an additional high-throughput sequence dataset from thousands of tumours and patient-matched normal tissues, which may contribute to improvements in diagnostic methods, treatment standards, and, ultimately, cancer prognosis [7]. In this study, the expression profiles of genes and miRNAs were downloaded from TCGA and analysed. First, the differentially expressed genes and miRNAs were identified, and then Gene Ontology (GO) classification and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis were performed. Last, the gene-miRNA co-expression network and TF-miRNA-gene network were established to study the molecular mechanisms and find new strategies for the prevention, diagnosis, and treatment of colon cancer (Fig. 1).

Fig. 1

Flowchart of study

Material and methods

Data resources and preprocessing

The level 3 miRNASeq and RNAseqV2 data of colorectal adenocarcinoma (COAD) were downloaded from the TCGA database (http://cancergenome.nih.gov/) by the TCGA Assembler package in R (version 3.1.1) [8, 9]. The miRNA and mRNA expression profiles were generated using Illumina HiSeq 2000 miRNASeq and Illumina HiSeq RNASeq sequencing platforms, respectively, and the RNA-Seq by the Expectation-Maximisation (RSEM) software package [10]. Raw count data and read count data were used to represent the mRNA and miRNA expression levels, respectively. There were 280 samples in total, including 272 colon adenocarcinomas and eight normal tissue samples; these samples were classified into two cohorts: a tumour group and a normal group. Each sample included the corresponding miRNAseq and RNAseq data. Because the data were obtained from TCGA, approval by an Ethics Committee was not needed. This study meets the TCGA publication guidelines.

Identification of differentially expressed genes and miRNAs

Differentially expressed genes and miRNAs were selected based on their fold change and adjusted p-values, which were generated by the DESeq package [11]. The inclusion criteria were set as follows: 1) FDRPadj< 0.05 and |log2(fold change)| > 1; and 2) gene and miRNA expression level > 0.

Screening of target genes of differentially expressed miRNAs

To understand the regulatory relationships between the differentially expressed miRNAs and genes, the target genes were predicted with the miRWalk2.0 database [12] (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/index.html) and those that had been validated by experiments were chosen. The intersections of genes from target genes and differentially expressed mRNAs were then obtained to construct the miRNA-gene co-expression network.

GO classification and KEGG pathway analysis

We used the DAVID tool (Database for Annotation, Visualisation, and Integrated Discovery) [13] for GO classification and KEGG pathway analysis of the differentially expressed genes. The DAVID tool is a web-accessible program that provides a comprehensive set of functional annotation tools for investigators to understand the biological meaning behind large lists of genes (http://david.abcc.ncifcrf.gov/). We selected the top 10 GO terms by their p value for analysis. We used R packages “ggplot2” and “clusterProfiler” to annotate modules with GO terms and to perform KEGG pathway enrichment analysis, respectively.

Construction of the miRNA-gene interaction network

First, we constructed the co-expression network from the identified differentially expressed genes and miRNAs based on their topological properties, such as a distance measure, using the WCGNA (Weighted Gene Co-expression Network Analysis) [14] package in R. Because we applied WGCNA to the gene and miRNA expression data in our analysis, the soft power threshold of β = 16 was selected, which is the smallest value that reaches level 0.9 on the independence scale (Fig. 2). Finally, a network graph was constructed and visualised using Cytoscape v3.3.0 [15] (http://cytoscape.org/). The core region was identified with MCODE (Molecular Complex Detection) [16], a Cytoscape plugin. MCODE can detect densely connected regions in large protein-protein interaction networks.

Fig. 2

Analysis of network topology for various soft-thresholding powers

TF-miRNA-gene network

TFmiR was used (http://service.bioinformatik.uni-saarland.de/tfmir) to identify the TFs that might regulate these differentially expressed genes and miRNAs and to construct a TF-miRNA-gene network. TFmi [17] is a web server used for constructing and analysing disease-specific TF and miRNA co-regulatory networks.

Results

Differentially expressed genes and miRNAs

In total, 3436 differentially expressed genes (1709 upregulated and 1727 downregulated) and 216 differentially expressed miRNAs (99 upregulated and 117 downregulated) were identified in colon cancer in comparison with normal tissues, with cut-offs of padj < 0.05 and |log2(fold change)| > 1. Volcano plots and heat maps were generated to select differentially expressed genes and miRNAs, as shown in Figs. 3–6. The top ten most differentially expressed genes and miRNAs are also shown in Table 1.

Fig. 3

Volcano plots for expression of differentially expressed genes (Red represents the up-regulated genes, while green represents the downregulated genes)

Fig. 6

Heat map showing hierarchical gene clustering analysis of the differentially expressed microRNAs. Values used are normalised with log10. (The right longitudinal axis: the names of miRNAs, the left longitudinal axis: the clustering information of miRNAs, the upper abscissa axis: the clustering information of samples, the under abscissa axis: the sample code of cancer and normal. Red represents the up-regulated miRNAs while green represents the down-regulated miRNAs).

Table 1

The top ten most differential genes and miRNAs

Gene symbol	Log2(fold change)	Padj	Change
CA7	–6.0	4.37E-41	down
SCARA5	–4.7	4.37E-41	down
CLEC3B	–4.3	1.53E-35	down
KRT80	7.1	1.63E-35	up
SLC4A4	–4.3	2.60E-35	down
TRIB3	3.5	5.79E-33	up
SPIB	–4.4	3.10E-29	down
CADM3	–4.9	4.26E-29	down
ADAMDEC1	–3.9	4.26E-29	down
C2orf88	–3.6	5.76E-29	down
hsa-let-7b	–3.2	2.44E-43	down
hsa-mir-99b	–2.9	1.65E-33	down
hsa-mir-197	–4.6	4.89E-33	down
hsa-mir-328	–5.2	1.33E-32	down
hsa-mir-125a	–3.6	3.69E-32	down
hsa-let-7d	–3.4	3.69E-32	down
hsa-mir-139	–4.9	3.25E-27	down
hsa-mir-374a	7.8	2.49E-26	up
hsa-mir-21	6.1	2.94E-24	up
hsa-mir-181a-1	–2.9	6.82E-23	down

Volcano plots for expression of differentially expressed microRNAs (Red represents the up-regulated miRNAs, while green represents the downregulated miRNAs)

Heat map showing hierarchical gene clustering analysis of the differentially expressed genes. Values used are normalised with log10. (The right longitudinal axis: the names of genes, the left longitudinal axis: the clustering information of genes, the upper abscissa axis: the clustering information of samples, the under abscissa axis: the sample code of cancer and normal. Red represents the up-regulated genes, while green represents the down-regulated genes)

Functional enrichment analysis of the differentially expressed genes

We performed GO classification and KEGG pathway enrichment analysis for differentially expressed genes. We found some significant biological processes, including regulation of cell proliferation, cell adhesion, biological adhesion, cell cycle process, response to wounding, and others. The 10 top-ranked GO categories are shown in Fig. 7. Our results revealed that the differentially expressed genes were mainly involved in the cytokine-cytokine receptor interaction signalling pathway, cell adhesion molecule (CAM) signalling pathway, calcium signalling pathway, neuroactive ligand-receptor interaction signalling pathway, complement and coagulation cascade signalling pathway, and others (Fig. 8). In addition, 84 genes were enriched in the cytokine-cytokine receptor interaction signalling pathway (Fig. 9).

Fig. 7

GO analysis of differently expressed genes. (Y axis: Statistically Enriched GO biological process terms, X axis: the number of differentially expressed genes in each GO Term)

Fig. 8

KEGG pathway enrichment analysis. (Y axis: Statistically Enriched KEGG pathway, X axis: the number of differentially expressed genes in each pathway)

Fig. 9

The signalling pathway of cytokine-cytokine receptor interaction signalling pathway. (Red represents the up-regulated genes, while green represents the downregulated genes)

Analysis of the miRNA-gene regulatory network

WGCNA is a well-established method designed for constructing co-expression networks. Many network analysis studies have reported that a weighted network retains more information and is more robust than an unweighted network [18, 19]. According to the node connectivity, genes can be further classified into hub genes. Hub genes are very important nodes and represent a small proportion of nodes with maximal information exchange with other nodes in the gene co-expression network. In total, there were 106 nodes and 121 pairs of interaction relationships in the network (Fig. 10). The network included 13 miRNAs and 93 genes. To further explore the hub genes in detail, we identified a highly-connected core region with MCODE. Through MCOMD algorithm analysis, we identified six molecular complexes whose correlation integral values were higher than 2. The complex with the highest score contained seven nodes and 13 edges. The hub genes includedGCNT4 (glucosaminyl transferase 4, core 2), EDN2 (endothelin 2, also called ET-2), miR-1295, miR-615, and others (Fig. 11).

Fig. 10

Co-expression networks between miRNA and genes. (Red represents the genes, while green represents microRNAs)

Fig. 11

Molecular complexes obtained by MCOMD algorithm analysis. (Red represents the genes, while green represents microRNAs)

We constructed the TF-miRNA-gene network of the differentially expressed miRNAs and genes using TFmiR, which is a web server, as shown in Fig. 12. LEF1 (lymphoid enhancer-binding factor 1) is the core TF that regulates the most differentially expressed miRNAs and genes. The most central mRNAs were MYC (also called c-myc) and WT1, and the most central miRNA was has-mir-34a.

Fig. 12

TF-miRNA-gene network in colon cancer. (Rectangles represent differentially expressed genes, hexagons for transcription factors, and ellipses for miRNAs. The lines indicate the interaction relationship. Hub genes are marked in red)

Discussion

In this study, we analysed the differences in the gene expression profiles of normal and colon cancer samples and investigated the functions and pathways with which differentially expressed genes were mainly associated. We also constructed a miRNA-gene regulatory network and TF-miRNA-gene network based on the closely related mRNA and miRNA expression data to better clarify the cellular mechanisms of colon cancer at the molecular level.

Here, we found that some hub genes, such asGCNT4, EDN2, and miR-1295, are massively altered in colon carcinoma. These genes and miRNAs can potentially be applied not only as novel biomarkers but also as therapeutic targets. For example, GCNT4 belongs to the glucosaminyl (N-acetyl) transferases family. Altered glycosylation is considered a universal cancer hallmark. Though little is currently known about GCNT4, the other genes (GCNT2 and GCNT3) of the glucosaminyl transferases family were identified to be related to colon cancer. Nakamura et al. showed that the expression of GCNT2 was aberrant, which is explained by promoter DNA hypermethylation. Hypomethylation of the GCNT2 variant 2 can reflect lymph node metastasis of CRC in the tumour and in normal tissues [20]. Additional, low GCNT3 expression was proven to be a promising prognostic biomarker for colon cancer, which could be used to identify early-stage colon cancer patients at high risk of relapse [21]. Therefore, we speculated that the GCNT4 may also have significant relationships with colon cancer. EDN2 encodes a member of the endothelin protein family of secretory vasoconstrictive peptides. The overexpression of ET-2 and ET-3 significantly suppressed the migration and invasion of human colon cancer cells [22]. According to the miRNA-gene network, GCNT4 and EDN2 were regulated by miR615 and miR-1295. MiR-615 was reported to be a tumour suppressor in many cancers [23-25]. Because of promoter hypermethylation, miR-615 was abnormally downregulated in pancreatic cancer cells, which suppressed pancreatic cancer cell proliferation, migration, and invasion by directly targeting IGF2 and other genes, such as the proto-oncogene JUNB, and by interfering with the insulin signalling pathway [24]. The functions of miR-1295 are unclear.

KEGG pathway analysis showed that these genes were mainly enriched in cytokine-cytokine receptor interaction. Cytokines are soluble extracellular proteins that are crucial intercellular regulators and mobilisers of cells engaged in innate and adaptive inflammatory host defences; cell growth, differentiation, and death; angiogenesis; and development and repair processes aimed at restoring homeostasis. The aberrant expression of some cytokines may be associated with colon cancer; for example, IL-6 is greatly upregulated in colon cancer, and its effect on epithelial cells is multifaceted [26, 27]. Interleukin 6 aggravated tumorigenesis through STAT3 signalling in colon epithelial cells, which enhanced epithelial cell survival and proliferation [28]. Interleukin 6 also downregulated the tumour suppressor p53 and upregulated the oncogene c-myc in colon epithelial cells and EMT [29].

Gene expression is usually regulated in one of two ways: regulation of the transcription from DNA to RNA by TFs, or regulation of RNA stability by miRNAs [30]. In our study, LEF1, MYC, WT1, and miR-34a were located in the central hub of our TF-miRNA-gene network. As a transcript factor, LEF1 is involved in the Wnt signalling pathway, which has been demonstrated to be involved in CRC. LEF-1 expression is involved in the presence of KRAS mutations and has prognostic value in colorectal carcinoma as a trend of worse overall survival [31]. LEF1 can target CDH1 and CTNNB1 genes, which closely interact in EMT activation [32]. LEF1 expression levels are associated with infiltration depth, lymph node, distant metastases, and advanced TNM stages of colon cancers as well as poor overall survival rate in patients with colon cancer [33]. As a negative regulator, miR34a was reported to inhibit colon cancer proliferation and metastasis by targeting platelet-derived growth factor receptor α (PDGFRA) [33]. MiR-34a inhibits cell migration and invasion by regulating the SIRT1/p53 pathway in humans [34, 35]. MYC mutations, overexpression, rearrangement, and translocation have been associated with a variety of haematopoietic tumours, leukaemias, and lymphomas, including Burkitt’s lymphoma. MYC plays an essential role in the regulation of many physiological processes, including cell cycle control, apoptosis, protein synthesis, and cell adhesion [36]. Among its related pathways are the PI3K-Akt signalling pathway and MAPK signalling pathway. Diseases associated with WT1 include Denys-Drash syndrome and Frasier syndrome. Among its related pathways are the integrated pancreatic cancer pathway and transcriptional misregulation in cancer. Although aberration expression of MYC and WT1 in colon cancer is frequently observed, future studies are required to explore the mechanism. To conclude, our study may provide useful information for understanding the TF-miRNA-gene network in CRC.

Overall, our findings will improve our understanding of the molecular mechanisms of colon cancer and aid in finding potential targets for diagnostic and therapeutic usage.