Co-expression Network Analysis Elucidated a Core Module in Association With Prognosis of Non-functioning Non-invasive Human Pituitary Adenoma

Non-functioning pituitary adenomas (NFPAs) are tumors with clinically challenging features since they have insidious progression. A complex network of gene interactions is thought to have roles in tumor formation and progression. Therefore, revealing the genetic network behind NFPA tumorigenesis is not only essential to attain further knowledge of tumor biology, but also plays a fundamental role in the development of efficacious treatment strategies. Differential co-expression network analysis is an outstanding approach for elucidation of groups of genes which show distinct co-expression patterns among phenotypes. In this study, we carried out a differential co-expression network analysis of NFPA-associated transcriptome dataset (n = 40) considering invasive (n = 22) and non-invasive (n = 18) phenotypes. Furthermore, we identified differentially co-expressed and co-regulated mRNA modules, which might be considered as potential systems biomarkers for NFPA prognosis and invasiveness. As a result, we have identified a novel 13-gene module, including CEACAM6, CYP4B1, EIF2S2, HID1, IFFO1, MYO18A, PDCD2, SGIP1, SWSAP1, and four unknown genes (A_24_P127621, A_24_P255786, A_24_P683553, and A_24_P916979), which was able to categorize the patients into two groups as invasive and non-invasive NFPA with distinct prognosis. The prognostic core module genes were associated with progression and prognosis of brain and glandular based cancers as well. Furthermore, these module genes were also expressed in blood, salivary gland, and spinal cord tissues. These results may provide the evidence on featured gene module which might play a prominent role in NFPA prognosis and sub-typing as effective biomarkers and therapeutic targets in the future.

Introduction

Pituitary adenomas (PAs) are the second most frequently reported primary brain tumors that are associated with increased mortality and morbidity (1, 2). Non-functioning pituitary adenomas (NFPAs) are the most common type of PAs that are not hormonally active (3). In contrast with the functioning pituitary adenomas (FPAs), which secrete excess levels of hormones that can lead easily to track endocrine syndromes, the detection of pituitary adenomas including NFPAs is clinically challenging. They are diagnosed usually in the context of mass effect (compress adjacent neurovascular tissues) leading to visual loss, dysfunction of the pituitary, and cranial neuropathies. PAs are usually non-invasive and benign neoplasms in nature, but some of them exhibit an aggressive attitude and represent invasive characteristics through cells from a neoplasm extend to the adjacent healthy tissues and infiltrating into them. Invasive PAs display local invasion, elevated risk of postoperative recurrence, and inadequacy of therapeutic response (4). Even in the perturbation of surrounding tissues, invasive PAs are not considered as malignant.

In the classification of tumors of the endocrine system, World Health Organization (WHO) offers several markers associated with invasiveness such as the Ki-67 proliferative index, the number of mitotic figures, and the expression of p53 (5). Moreover, FSH-β, LH-β, and/or α-subunit immunostain positivity and expressivity of SF-1, GATA-2, and ER-α transcription factors were proposed as the indicators of the gonadotroph lineage in NFPAs, and null cell adenomas were defined as negative immunoreactivity for both pituitary hormones and pituitary transcription factors (5, 6).

The expression profiling studies at RNA level intending the identification of molecular markers associated with NFPA invasiveness are very limited. The upregulation of MYO5A, and downregulation of E-cadherin (CDH1) and H-cadherin (CDH13) were proposed as markers of invasiveness at mRNA level (7, 8). Several microRNAs including miR-135a, miR-140-5p, miR-582-3p, miR-582-5p, and miR-938 were found as overexpressed in the invasive state, and they were also associated with tumor size and tumorigenesis (9). Moreover, the expression levels of lncRNAs, maternally-expressed gene 3 (MEG3) and Hox transcript antisense intergenic RNA (HOTAIR), were correlated with NFPA development and invasion (10).

Despite the proposed molecular markers and immunohistochemical indicators of invasiveness, early prediction of invasive NFPA remains considerably arduous since there is no consistent pattern that differentiates invasive from non-invasive PA. Thus, invasive PA needs to be properly defined with molecular markers in order to identify patients at increased risk of early recurrence or subsequent tumor progression. The present markers and classification systems are insufficient to guide precise diagnostic and therapeutic decisions in pituitary tumor invasiveness.

Differential co-expression is defined as the altered co-expression patterns of genes between two phenotypes and represents significant potential to identify gene clusters associated with the phenotype of interest (11). In several studies, co-expression analyses were performed to elucidate the PA metastasis, invasion and progression. Zhang et al. identified co-expressed gene pairs for prolactin (PRL) secreting pituitary tumor metastasis and hypothesized that SLC2A11, TENM1, IPO7, and CHGB are associated with metastasis in prolactinoma cases (12). Co-expression of somatostatin and dopamine receptors were also investigated and invariable loss in expression of both receptors in invasive growing corticotroph adenomas was reported (13).

In this study, we aimed to determine whether invasive and non-invasive PAs display differences in mRNA expression profiles through a differential co-expression network analysis framework. A novel gene module, which was differentially co-expressed in non-invasive PA when compared to invasive PA state, was presented, and its prognostic performance was tested in several tumor types through survival analyses. Topological and functional enrichment analyses were performed to elucidate the molecular mechanisms. Transcriptional regulatory elements (i.e., transcription factors and miRNAs) associated with the prognostic gene module were investigated, and expression levels of the module genes in various tissues were also screened.

Methods

Selecting Gene Expression Profiling Dataset

In order to gather NFPA associated gene expression profiling datasets, we comprehensively screened the publicly available functional genomics data repositories including ArrayExpress (14), NCBI Gene Expression Omnibus (GEO) (15) and The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/). Among three transcriptome datasets, GSE63357 (16), GSE77517 (17), and E-TABM-899 (7), carried out with NFPA samples (with no hormone secretion), we preferred to employ E-TABM-899 considering its large sample size and sampling characteristics. Expression profiles of mRNA and related clinical data of invasive and non-invasive PA were downloaded from ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/). Gene expression profiles generated using Agilent whole human genome oligo microarrays were used as a training set to construct co-expression networks, identify hub genes, and differentially expressed genes (DEGs) in this study. Transcriptome data was composed of 40 samples classified as invasive (n = 22) or non-invasive (n = 18). NFPA samples were characterized in terms of clinical and pathological features including age, sex, hormonal secretion (LH, FSH), immunohistochemical staining, tumor volume (varied between 0.5 and 77 cm³), tumor grade, follow-up time, and outcome (recurrence, stably remnant, and remnant). Immunostaining for FSH-β and LH-β and/or α-subunit was used for tumor characterization. According to WHO 2017 classification of tumors of the pituitary gland (5), the dataset was composed of 38 samples with gonadotroph lineage (positive immunostains FSH-β and LH-β and/or α-subunit) and 2 samples with null cell adenoma (negative for immunostains) (5, 6). Tumor size and invasiveness were determined with regard to preoperative magnetic resonance imaging and perioperative findings, by employing the modified Hardy criteria (18), which is grade I, if tumor < 10 mm, it is enclosed microadenoma; grade II, if tumor > 10 mm, it is enclosed macroadenoma; grade III, localized perforation of the sellar floor; and grade IV, diffusive destruction of the sellar floor. Tumor recurrence was classified according to radiological findings of re-growth of the tumor remnant or new tumor growth (in patients with no surgical remnants). Remnant tumors were accepted as stable if there was no sign of growth on two MRI within the 1-year interval and no clue of disease reactivation. Follow-up times were years that passed after surgery and varied between 3 and 10.5 years. This dataset included tumors in varied grades, 1 of grade I, 18 of grade II, 6 of grade III, and 15 of grade IV. Grade III and IV tumors, as well as one Grade II tumor exhibiting extradural extension, were accepted as invasive.

Data Processing and Identifying Differentially Expressed Genes

Pre-processing of microarray data was performed in R using expression profiles (E-TABM-899.eSet.r) downloaded from ArrayExpress database (14). The expression profiles were normalized through quantile normalization. DEGs between invasive and non-invasive PAs were identified from the normalized expression values by using Linear Models for Microarray Data (LIMMA) package (version 3.34.5) (19). The Benjamini-Hochberg method was used to control the false discovery rate. In addition, cross-validation of the DEGs was performed using the two-sample t-test. To determine the statistical significance, adjusted p < 0.05 was used, and genes that showed at least 10% change in expression were selected for further consideration in the network construction.

Pathway and Functional Enrichment Analysis for DEGs

The pathway and functional enrichment analyses were performed using ConsensusPathDB (Release 32) (20) to identify functional annotations significantly associated with the DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (21) database was preferably used as the pathway database. The Gene Ontology (GO) terminology (22) was employed as the source for annotating the molecular functions and biological processes. P-values were obtained via Fisher's Exact Test. Benjamini-Hochberg's correction was used as the multiple testing correction techniques, and enrichment results with adjusted p < 0.01 were considered statistically significant. Comparison of the significant GO terms and KEGG pathways was performed to identify the different biological functions and processes between the invasive and non-invasive samples.

Co-expression Network Construction for Invasive vs. the Non-invasive States

Co-expression networks were constructed for invasive and non-invasive PA states as previously described (11). Initially, among the 29,754 genes detected in the pituitary samples, DEGs were identified. In the case of the repetition in expression values of DEGs, the mean expression values were computed and used in further analyses. Pearson's correlation coefficient (PCC) was employed to determine the co-expression levels between DEG pairs in each state. Benjamini-Hochberg's correction was used as the multiple testing correction techniques, and the gene pairs were regarded as co-expressed if |PCC| ≥ 0.80 and adjusted p ≤ 0.05.

To identify differential co-expression profiles between two states, the following criterion was employed:

$\begin{array}{l}|({\text{PCC}}_{\text{NONINV}}-{\text{PCC}}_{\text{INV}})/{\text{PCC}}_{\text{INV}}|\ge 1\end{array}$

where PCC_NONINV and PCC_INV are the PCC of a DEG pair in non-invasive and invasive states, respectively. Resultant gene pairs were selected to construct the co-expression networks, ICON (invasive co-expression network) and NICON (non-invasive co-expression network), representing invasive and non-invasive states of PA, respectively.

Identification of Differentially Co-expressed Gene Modules

To elucidate highly connected network modules, co-expression networks were analyzed using MCODE (23) plugin of Cytoscape (v.3.5.1) (24). Modules with at least 10 nodes (genes), average connectivity ≥10 and clustering coefficient ≥ 0.5 were considered as differentially co-expressed modules between invasive and non-invasive states. Modules were further analyzed with Network Analyzer (25) and Cytohubba (26) plug-ins to determine the significantly altered hub genes with respect to local and global topological metrics (i.e., degree and betweenness centrality).

Sub-type Clustering Performance of Co-expressed Gene Modules

Considering the expression profiles of module genes, Principal component analysis (PCA) was performed to cluster invasive and non-invasive samples. The principal components explaining at least 85% of the total variance was considered in the determination of clustering performance, i.e., the ability hereupon “sub-typing” of PA.

Prognostic Power and Survival Analyses

PCA was performed based on gene expression profiles of module genes, and samples were categorized into three clusters using k-means algorithm taking into consideration the first three principal components (explaining 98% of the total variance). Cox regression was performed via MedCalc Statistical Software (v.18.2.1, Ostend, Belgium) to elucidate survival-associated genes. Analysis results were visualized and compared through Kaplan-Meier plots using the follow-up times of invasive and non-invasive groups. The prognostic power of the module was determined using the log-rank test p-value, the hazard ratio (HR), and its confidence intervals (CI).

In silico Validation in Related Tumor Types

We used SurvExpress (27) web-based biomarker validation tool in order to test the validity of the proposed prognostic gene module in various cancer types. Validation analyses were carried out through independent RNA-seq datasets obtained from TCGA (28) and microarray datasets from NCBI-GEO. Core module genes were tested through glandular based cancers including adrenocortical carcinoma (TCGA, n = 77), breast invasive carcinoma (TCGA, n = 502), and ovarian serous cystadenocarcinoma(TCGA, n = 247); brain tumors including glioblastoma (TCGA, n = 660), meningioma (NCBI-GEO, GSE16581, n = 68), and medulloblastoma (NCBI-GEO, GSE37418, n = 76) (29, 30); and other cancers including prostate adenocarcinoma (TCGA, n = 497) and lung adenocarcinoma (TCGA, n = 475). For each dataset employed, samples were partitioned into low and high-risk groups according to their prognostic index, survival multivariate analyses, and risk assessments. The prognostic performance of modules in each dataset was determined using Kaplan-Meier plots, log-rank test p-values, hazard ratios (HR), and their confidence intervals (CI). Differences in gene expression levels between high and low-risk groups were presented by box-plots and statistical significance of the difference was calculated through t-test. Heatmaps were also created to clarify the correlation of the survival analysis with gene expression levels. Samples were sorted by their prognostic index and genes were clustered by using Euclidean distance.

Screening of Gene Expressions in Varied Tissues

The expression levels of modules genes in easily collectible tissues were obtained from The Genotype-Tissue Expression (GTEx) project (31), which provides a genome-scale expression profiling (RPKM) of 51 normal human tissues, cells, and fluids.

Results

Transcriptional Profiling in Non-functioning Non-invasive Human Pituitary Adenoma

Transcriptome data (E-TABM-899) was composed of 40 NFPA samples, which were characterized in terms of clinical and pathological features including age, sex, hormonal secretion (LH, FSH), immunohistochemical staining (FSH-β and LH-β and/or α-subunit), tumor volume (0.5–77 cm³), tumor grade (I-IV), follow-up time (3–10.5 years), and outcome (recurrence, stably remnant, and remnant). Tumor size and invasiveness were determined with regard to preoperative magnetic resonance imaging and perioperative findings, by employing the modified Hardy criteria (18), Tumor recurrence was characterized according to radiological findings of re-growth of the tumor remnant or new tumor growth (in patients with no surgical remnants). As a result, the samples were classified into invasive (n = 22) and non-invasive (n = 18) groups.

Further verification of sample classification was carried out employing the expression values of MKI67 gene, which encodes Ki-67 protein accepted as a labeling index of proliferation (32). MKI67 expression levels in two groups were compared and a statistically significant difference (p = 0.02) was observed between mean expression levels of invasive (88.49 ± 2.32) and non-invasive groups (81.62 ± 1.56; Figure 1).

FIGURE 1

Figure 1. The expression values of MKI67 gene, which encodes Ki-67 protein, between tumor groups.

Identification and Functional Annotation of Differentially Expressed Genes (DEGs)

Following up a standard protocol including quantile normalization of raw data, statistical comparisons of expression levels through “limma” package, and multiple testing correction with Benjamini-Hochberg method, a total of 4,091 genes were identified as differentially expressed between two phenotypes, i.e., invasive and non-invasive pituitary adenoma, considering a statistical confidence level of p < 0.05. In addition to p-values, we also considered at least 10% change in mean expression levels (FC < 0.9 for down-regulation, FC > 1.1 for up-regulation) as a criterion in the identification of DEGs and determination of their direction of regulation. This filtering resulted in 641 DEGs of which 145 were down-regulated and 496 were up-regulated.

The functional annotations of the down-regulated DEGs were significantly enriched with transmembrane transport of organic substances (especially the cations), disassembly of protein complexes and its regulation, endocytic recycling, and response to copper ion; whereas genes associated with responses to lipids, DNA damage stimulus, and cytokines, regulation of protein localization, proteasomal protein catabolism, endoplasmic-reticulum-associated protein degradation (ERAD) pathway, regulations of telomere maintenance and oxidoreductase activity, cell aging, and oxidative demethylation were up-regulated (Figure 2A). The major KEGG pathways significantly associated with DEGs (p < 0.01) were axon guidance (17 genes), ABC family mediated transport (8 genes), SUMOylation of DNA damage response and repair (6 genes), estrogen receptor alpha signaling pathway (5 genes), imatinib pathway (3 genes), and warfarin pathway (2 genes; Figure 2B).

FIGURE 2

Figure 2. Biological process and pathway enrichments of DEGs. (A) Biological process enrichments of up- and down-regulated DEGs. (B) Pathway enrichments of DEGs.

Detection of Differentially Co-expressed Gene Modules

The possible co-expression patterns of 641 DEGs were analyzed through the employment of PCCs in invasive and non-invasive phenotypes separately to construct phenotype-specific co-expression networks. This resulted with a denser, centralized, and modular co-expression network in invasive phenotype (ICON) with 183 links among 162 genes (Figure 3A) when compared to the non-invasive co-expression network (NICON), which includes 2,385 links among 562 genes (Figure 3B). Not only the number of genes exhibiting co-expression pattern was almost 4-fold higher in the non-invasive state, but also the network centralization (0.835) and the clustering coefficient (0.348) were higher in NICON when compared to those in ICON (0.394 and 0.002, respectively).

FIGURE 3

Figure 3. Differentially co-expressed gene networks of invasive and non-invasive states. (A) The co-expression network in invasive phenotype (ICON) with 183 links among 162 and topological features of ICON (B) the non-invasive co-expression network (NICON), which includes 2,385 links among 562 genes and topological features of the network (C) Functional annotations of ICON hub genes (D) Functional annotations of NICON hub genes.

Hub genes of ICON network (ARL5, ASTN2, C11orf39, CTNNAL1, EIF2S2, NPTX2, RALGPS2, and SLC12A9) were enriched in cytolysis, Ras protein and Rho protein signal transduction and small GTPase signal transduction biological processes (Figure 3C). Hub genes of NICON network (ASTN2, C11orf39, CEACAM6, CTNNAL1, GZMM, GLT8D2, NPTX2, RALGPS2) were enriched in processes such as urogenital system development, epoxygenase P450 pathway, thyroid gland development, protein tyrosine kinase activity, drug metabolic processes, protein kinase B signaling, intracellular receptor signaling pathway (Figure 3D).

Comparative analysis of the phenotype-specific co-expression networks yielded with 10 differentially co-expressed gene modules (ranging in size from 5 to 69 genes), which were activated in non-invasive phenotype when compared to invasive. Differentially co-expressed gene modules were named as M1 (Module 1), M2, M3, M4, M5, M6, M7, M8, M9, and M10. Module sizes and their associations with biological processes and PA were tabulated (Table 1). M1 was composed of 41 genes and these genes were enriched in regulation of cell adhesion. M2 had 45 genes that were enriched in GTPase activator activity and nucleoside triphosphatase regulator activity. M3 included 14 genes that related to Actin filament binding and Ras protein signal transduction processes. M4 had 69 genes which were enriched in DNA damage checkpoint and regulation of cellular response to stress. There were 10 genes in M5 and enriched in DNA damage checkpoint process. M6 had 8 genes that were related with metal ion binding process. M7 included 8 genes that enriched in pyrophosphatase activity, hydrolase activity, and acting on acid anhydrides. M8, M9, and M10 included 5 genes that were enriched in nervous system development, Ras and GTPase mediated signal transduction, and protein complex binding, respectively. Core module included 13 genes which were enriched in ATPase activity and regulation of endocytosis. These modules have been associated with neither PA diagnosis nor PA prognosis in literature.

TABLE 1

Table 1. Identified modules and their features.

To determine the possible associations of the identified co-expression modules with clinicopathological features (i.e., age, tumor size, tumor grade, hormonal status, follow-up time, and clinical outcome), we conducted the Principal Components Analyses (PCA) as well as Kaplan Meier survival analyses, and detected significant correlations with tumor grade and sub-type for several gene modules.

Three modules, i.e., M3 with 14 genes associated with actin filament binding and Ras protein signal transduction, M6 with 8 genes associated with metal ion binding, and M7 with 8 genes associated with pyrophosphatase and hydrolase activities, presented considerable accuracy in discrimination of two sub-groups of invasive and non-invasive states in PA (Figure 4). Considering the similar discriminatory behavior of these modules, we combined them and narrowed the co-expression patterns with a higher PCC cut-off (increasing from 0.8 to 0.9) to obtain a “core module” consisting of the minimum number of genes with higher discrimination ability. This resulted with a core module of 13 genes (CEACAM6, CYP4B1, EIF2S2, HID1, IFFO1, MYO18A, PDCD2, SGIP1, SWSAP1, A_24_P127621, A_24_P255786, A_24_P683553, A_24_P916979), which indicates higher discrimination accuracy, since clusters representing invasive and non-invasive states were clearly separable and distinguishable without any overlaps in PCA plots (Figure 4). Interestingly, the genes in the core module were not enriched in any molecular pathway or biological process as a whole; but some sub-groups were enriched in ATPase activity and regulation of endocytosis (p < 0.001). Notably, four genes (A_24_P127621, A_24_P255786, A_24_P683553, A_24_P916979) located in distinct chromosomes were not annotated with any molecular function yet (Table 2).

FIGURE 4

Figure 4. Principle Components Analysis (PCA) plots of differentially co-expressed gene modules (M3, M6, M7, and core module) that have significant prognostic performance.

TABLE 2

Table 2. Differentially co-expressed genes represented in the core module.

The Prognostic Performance of the Core Module Between Different PA States

The genes in the core module indicated significantly different co-expression patterns between invasive and non-invasive phenotypes of PA (Figure 5). The clustering of samples based on the expression profiles of core module genes via the first principle component (PC1) of the PCA analysis (describing 89.8% of the total variance) and k-means algorithm resulted with three distinct sub-groups of patients (Figure 6A). The cluster representing lower loadings on PC1 consisted of only invasive and grade IV tumor samples; whereas the cluster with higher loadings on PC1 included non-invasive and grade II tumors. Furthermore, survival analysis with the Cox regression model and Kaplan-Meier estimates indicated the predictive power of patient survival among the sub-types (Figure 6B). There was a significant difference among groups in terms of survival probability estimates (p = 0.0024, log-rank test), and the average survival probability of non-invasive and grade II tumors was significantly higher than that of invasive and grade IV tumors with a hazard ratio of 2.96 (p = 0.011).

FIGURE 5

Figure 5. The degree of co-expression of core module genes between invasive and non-invasive states. Upper-triangle represents invasive state whereas lower-triangle shows non-invasive state.

FIGURE 6

Figure 6. Clustering of sub-types of NFPA with PCA analysis. (A) The clustering of samples via first principle component of expression profiles of core module genes (B) survival analysis by Cox regression model and Kaplan-Meier estimates indicated the predictive power of patient survival among the sub-types.

Transcriptional Regulators of the Core Module Genes

In order to elucidate the regulatory mechanism behind the co-expression pattern of module genes and to evaluate the condition-specific expression pattern alterations, we performed serial analyses to link the key regulators of transcriptional control, i.e., transcription factors (TFs) and microRNAs (miRNAs), to genes in the core module. A total of 10 TFs were found as regulators of the module genes and GATA1, GATA2, ETS1, ESR1, and PRDM14 were top five regulators with highest degree values (Figure 7A). The majority of the core module genes were co-regulated by ETS1 and GATA2. Furthermore, among module genes EIF2S2, IFFO1, and PDCD2 were co-regulated by the same TFs. In contrast to TFs, core module genes were regulated by distinct miRNAs (Figure 7B). Among those, mir-335-5p was the only common regulator between two gene pairs (HID1 and MYO18A). Members of the core module were found in varied chromosomal locations. CEACAM6 and SWSAP1 were located on chromosome 19 whereas HID1 and MYO18A were on chromosome 17, CYP4B1, and SGIP1 were on chromosome 1, PDCD2 was on chromosome 6, A_24_P127621 was on chromosome 7, A_24_P683553 was on chromosome 9, A_24_P255786 was on chromosome 10, IFFO1 was on chromosome 12, EIF2S2 was on chromosome 20, and A_24_P916979 located on chromosome X.

FIGURE 7

Figure 7. Transcriptional regulators of prognostic genes. (A) Transcription factors that regulate prognostic core module genes (B) microRNAs (miRNAs) regulating the core module genes.

The Specificity of the Core Module to NFPA

In order to analyze NFPA-specific expression and prognostic performance of the core module, we carried out several analyses over various human cancers with or without tissue similarity with PA. First of all, we tested its prognostic performance in other glandular based cancers including adrenocortical carcinoma, breast invasive carcinoma and ovarian serous cystadenocarcinoma (Figure 8), and observed that the core module presents high prognostic performance in adrenocortical carcinoma (HR = 7.22, p = 4.69 × 10⁻⁶), breast invasive carcinoma (HR = 2.60, p = 1.84 × 10⁻⁴), and ovarian serous cystadenocarcinoma (HR = 2.02, p = 4.56 × 10⁻³). Secondly, we tested the performance of core module in brain tumors (Figure 9), and significant results were obtained in glioblastoma (HR = 3.54, p < 10⁻²⁰), meningioma (HR = 10.54, p = 8.90 × 10⁻⁴), and medulloblastoma (HR = 14.25, p = 6.02 × 10⁻³). Lastly, we tested the prognostic performance of the core module in two other unrelated tissue but highly prevalent cancers, namely lung, and prostate adenocarcinoma. This time survival analyses pointed out insignificant performance in both lung (p = 0.08) and prostate (p = 0.15) adenocarcinomas (Figure S1). As a consequence, our core module could be considered as prognostic for both glandular-origin cancers and brain cancers, but not for lung and prostate cancers.

FIGURE 8

Figure 8. Prognostic power of core module through related tissue type cancers. Core module was also prognostic for (A) Adrenocortical carcinoma, (B) breast adenocarcinoma, (C) ovarian serous adenocarcinoma.

FIGURE 9

Figure 9. Prognostic power of core module through 3 types of brain tumors. (A) glioblastoma, (B) meningioma, (C) medulloblastoma.

Tissue-Specific Expressions of Core Module Genes

We searched for the expression levels of module genes in easily collectible tissues (especially not need to biopsy or taking by minimally invasive operations) in order to decipher whether these genes can be detectable at mRNA level. Core module genes have detectable expression values in varied tissue types like blood, kidney, lung, spleen, spinal cord, saliva, pancreas, and uterus. It was also difficult to obtain a cumulative behavior in expression patterns of the genes in the same tissue; however, we found that every core module genes have considerable expression values in three tissue types, which are blood, salivary gland, and spinal cord (Figure 10).

FIGURE 10

Figure 10. Tissue-specific expressions of core module genes.

Discussion

There is a great need for interpretation of omics data to elucidate the multi-layered molecular mechanisms lay behind the invasion of PA. However, omics datasets associated with PA invasiveness is scarce in the publicly available databases. Even in TCGA, which is the most comprehensive and multi-dimensional database for cancer genomics with 33 cancer types, genomics and/or transcriptomics of PA were not considered. The scarcity of omics datasets indicates the urgent need for further efforts on mRNA profiling in NFPA.

Despite designated markers associated with invasive behavior of PA tumors such as the Ki-67 proliferative index, the number of mitotic figures, and the expression of p53 (48, 49), the early prediction of PA invasiveness is still considerably hard, since the patient's state tends to be delayed due to the clinical symptoms of NFPA which are not obvious in the early stage. In the era of precision medicine, accurate biomarkers for cancer-specific prognosis and early prediction were urgently required and this convincingly enhances decision making for patient management.

Instead of the construction of the commonly used biological networks in literature like PPIs and metabolic networks, we constructed gene co-expression networks since it has several advantages such as almost complete coverage of human genes, little knowledge bias obtained from the published literature, and the ability to create cancer type–specific networks. Identification of changes in co-expression patterns of genes among invasive and non-invasive samples provides information about the invasiveness of PA-specific gene modules. Previous studies in literature objected to predict candidate prognostic biomarkers or protein-protein interaction networks constructed around differentially expressed genes in PA, not co-expressed gene modules. The altered co-expression schemes of module genes among invasive to the non-invasive state were neglected in these studies. Here, we made a differential co-expression analysis to elucidate PA genes and their co-expression pairs. As a result, the novel co-expressed gene module presented here may be regarded as “systems biomarkers” that lead to the design of effective therapeutic strategies in PA prognosis.

In the present study, co-expressed gene networks were constructed for invasive and non-invasive states (ICON and NICON, respectively) separately. ICON hub genes were enriched mostly in signaling processes like Ras, Rho, and small-GTPase signaling. It was reported that Ras small GTPase had a key signal transduction role in activating MAPK cascade (50). Rho GTPases are activated by G-protein-coupled receptors (GPCRs) and play critical roles in the invasion and metastasis of cancer cells (51). Hub genes of NICON were enriched in biological processes like urogenital system development, thyroid gland development, tyrosine kinase activity, protein kinase B signaling, and intracellular receptor signaling. The pituitary is part of the endocrine system and it was expected to found urogenital system development and thyroid gland development processes considering the secretion of TSH and FSH in the pituitary. Protein kinase B (Akt) signaling was found as over-expressed and over-phosphorylated in human PAs and proposed a potential role for Akt with respect to p27 deregulation (opposite direction to each other) (52).

The differential co-expression profiling followed by clustering and survival analyses resulted with a prognostic core module composed of 13 genes; namely, CEACAM6, CYP4B1, EIF2S2, HID1, IFFO1, MYO18A, PDCD2, SGIP1, SWSAP1, and 4 unknown genes (A_24_P127621, A_24_P683553, A_24_P255786, A_24_P916979). The unknown genes, A_24_P127621, A_24_P683553, A_24_P255786, and A_24_P916979, have specific gene sequences and distinct chromosomal locations (chromosome 7, chromosome 9, chromosome 10, and chromosome X, respectively); however, there was no known functional annotation of them in any of the annotation sources. Further in-depth experimental research studies might be performed to elucidate their function and association with PA tumorigenesis and/or invasiveness. The high prognostic power of the core module was not specific to NFPA but also observed in cancers with the same tissue origin such as glandular-origin cancers and brain cancers.

Among the module genes, CEACAM6 encodes a protein that is a member of the carcinoembryonic antigen (CEA) family whose crew are cell surface glycoproteins (33). Members of this family have a function in cell adhesion and they are widely employed as serum tumor markers for determinations of carcinoma (53). Moreover, over-expression of CEACAM6 modulates cancer progression through aberrant cell differentiation, anti-apoptosis, cell growth, and resistance to therapeutic agents (54). In addition, CEACAM6 over-expression in varied malignancies promoting cell invasion and metastasis, therefore, it represents a characteristic advantage of tumor cells which are responsible for an invasive phenotype (55). It is also a major determinant for malignant phenotype of pancreatic cancer through over-expression of the positive regulation of epithelial–mesenchymal transition (35, 36, 56).

CYP4B1 gene encodes a member of the cytochrome P450 superfamily of enzymes which is responsible for the oxidation of steroids, fatty acids, and xenobiotics therefore involved in reactions in drug metabolism (37). CYP4B1 was reported as procarcinogen for some chemicals in bladder carcinoma, down-regulated in esophageal squamous cell carcinoma and up-regulated in breast tumors when comparing the surrounding healthy tissues (38, 57, 58).

EIF2S2 encodes eukaryotic translation initiation factor 2 and functions as a vehicle in the early phases of protein synthesis. The deletion of EIF2S2 gene has been reported with suppression of testicular germ cell tumor incidence and recessive lethality in mice (39). It is also reported that EIF2S2 is up-regulated and involved in major enriched genes in epithelial ovarian carcinoma (EOC) samples and cervical cancer vs. normal group (53, 59, 60). In ER-negative breast cancer, EIF2S2 gene was found as a member of novel loci (61).

HID1 gene encodes a protein functioning in the trafficking of cargos that are important for the sorting/biogenesis/maturation of dense core vesicles and plays an important role in the development of cancers in a broad range of tissues (40, 62). It is alternatively named as DMC1 (down-regulated in multiple cancer 1). The loss of expression of HID1 was reported in breast, cervical, lung, thyroid, renal, and gastrointestinal cancer cell lines (62).

IFFO1 gene is a member of the intermediate filament family which include essential components of the cytoskeleton and nuclear envelope (41, 42). In ovarian cancer, its promoter was hypermethylated and IFFO1 was proposed as a biomarker (63). In addition, it has been identified as a strong prognostic indicator of breast cancer (63). The hypomethylation and up-regulation of IFFO1 have been also reported in endometrioid endometrial adenocarcinoma samples (64).

MYO18A encodes a protein which is part of a complex that unifies lamellar actomyosin bundles and needed for cell migration. MYO18A has been found highly expressed in metastatic prostate cancer and its knockdown affects the cytoskeleton and cell migration (65). MYO18A gene fusion has been identified in an acute cell leukemia patient (66).

PDCD2 (programmed cell death 2) gene encodes a nuclear protein expressed in a variety of tissues and its expression is controlled by transcriptional repressor BCL6. It plays an important role in cell death and/or in the regulation of cell proliferation. PCDC2 has been found as a tumor suppressor and involved in the pathogenesis of osteosarcoma. Aberrant expression of PDCD2 is associated with many tumors, such as leukemia and gastric cancer (44, 67). PDCD2 is an important predictor of clinical relapse in acute leukemia patients (67). Plus, loss of PDCD2 expression could induce gastric cancer development and progression through cell growth arrest at the early S phase of the cell cycle and reported as a putative tumor suppressor in gastric stromal tumors (68).

SGIP1 (SH3 domain GRB2 like endophilin interacting protein) functions as an endocytic protein that has effects on signaling in neuronal systems including energy homeostasis (45, 46). Hypomethylation and retrotransposition of SGIP1 have been reported in colorectal cancer samples (69). SGIP1 has been found as a potential therapeutic target for obesity- and diabetes-related symptoms, since the selective reduction of the expression of SGIP1 consequenced with inhibition of food intake and the decrement of body weight in rat models (70). SGIP1 has been identified as the second most significantly upregulated gene in the colon and rectal cancers samples (71). The expression level of SGIP1 has been shown significantly decreased in gastric cancer when compared to control samples (72).

SWSAP1 (SWS1-Associated Protein 1) gene is involved in homologous recombination repair by binding single-stranded DNA activity and it also has ATPase activity. Depletion of SWSAP1 contributes to defects in homologous recombination repair and Knockdown of SWSAP1 ends up with increased cell sensitivity to the DNA damaging (47).

Remaining core module gene members (A_24_P127621, A_24_P255786, A_24_P683553, A_24_P916979) were unknown with their unknown functions. The present study may contribute to assigning function/ functions to these genes. Known module genes were briefly explained with their functions and roles in various types of cancers. But there was no research found that could report module members for PA progression and prognosis. Therefore, identified module genes and module itself are novel for PA prognosis and they might be considered as the indicator of the non-invasiveness.

While our study provides some insight into the prognosis of PA in different states, more efforts should be performed to validate them clinically and extend our findings. We examined the prognosis of PA through core module genes in the context of gene co-expression networks. Results were in silico validated through various cancer types and it can be concluded that the core module is novel for PA prognosis. In addition, core module has also found as prognostic for brain tumors (glioblastoma, meningioma, and medulloblastoma) and glandular- originated tumors (adrenocortical carcinoma, breast invasive carcinoma, and ovarian serous cystadenocarcinoma). Our findings are based on TCGA, ArrayExpress, and GEO data (module investigation and validation), so a critical extension of this work would be to learn whether the patterns can be recapitulated by clinical trials. All in all, future efforts should be carried out to incorporate this systems-level understanding of prognostic genes into the practice of construction of effective clinical prognostic vehicles.

Conclusion

Pituitary adenomas are the most common intracranial tumors in the central nervous system. They are accepted as benign in general but somehow tumors can exhibit gross invasion into surrounding tissues rarely. Invasion results in resistance to conventional treatment methods and leading to early and frequent recurrences. In order to shed light on the multi-layered molecular mechanism lay behind the invasion of PA, there is a great need for omics-level data and their integration into meta-analyses. Our study used gene co-expression analysis to construct a gene co-expression network, elucidate and validate network modules associated with the invasiveness and prognosis of PA. Eventually, a core module with 13 genes including CEACAM6, CYP4B1, EIF2S2, HID1, IFFO1,MYO18A, PDCD2, SGIP1, SWSAP1, and 4 unknown genes (A_24_P127621, A_24_P255786, A_24_P683553, A_24_P916979) were identified and in silico validated in association with the indicator of invasiveness and prognosis of PA, plus some related cancer types. The research not only increases the theoretical knowledge but also provides a novel prognostic module and therapeutic strategy for PA.

Data Availability

Publicly available datasets were analyzed in this study. This data can be found here: https://www.ebi.ac.uk/arrayexpress/experiments/E-TABM-899/.

Author Contributions

BA and KYA designed the data analysis framework. BA performed the data analyses and evaluated the results. KYA conceived and directed the study. BA drafted the manuscript. KYA revised the manuscript. All authors read and approved the final manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2019.00361/full#supplementary-material

Figure S1. Prognostic performance of core module genes through lung and prostate cancers. Survival analyses indicated insignificant prognostic performance in both lung (p = 0.08) and prostate (p = 0.15) carcinomas.

References

1. Ostrom QT, Gittleman H, Fulop J, Liu M, Blanda R, Kromer C, et al. CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2008-2012. Neuro-oncol. (2015) 17(suppl_4):iv1-iv62. doi: 10.1093/neuonc/nov189

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Fernandez A, Karavitaki N, Wass JA. Prevalence of pituitary adenomas: a community-based, cross-sectional study in Banbury (Oxfordshire, UK). Clin Endocrinol. (2010) 72:377–82. doi: 10.1111/j.1365-2265.2009.03667.x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Mercado M, Melgar V, Salame L, Cuenca D. Clinically non-functioning pituitary adenomas: pathogenic, diagnostic and therapeutic aspects. Endocrinol Diabetes Nutr. (2017) 64:384–95. doi: 10.1016/j.endien.2017.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Landeiro JA, Fonseca EO, Monnerat ALC, Taboada GF, Cabral GAPS, Antunes F. Nonfunctioning giant pituitary adenomas: invasiveness and recurrence. Surg Neurol Int. (2015) 6:179–89. doi: 10.4103/2152-7806.170536

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Llyod RV, Osamura RY, Klöppel G, Rosai J (eds.). WHO Classification of Tumours of Endocrine Organs. Lyon: International Agency for Research on Cancer (2017). p. 355.

6. Lopes MBS. The 2017 World Health Organization classification of tumors of the pituitary gland: a summary. Acta Neuropathol. (2017) 134:521–35. doi: 10.1007/s00401-017-1769-8

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Galland F, Lacroix L, Saulnier P, Dessen P, Meduri G, Bernier M, et al. Differential gene expression profiles of invasive and non-invasive non-functioning pituitary adenomas based on microarray analysis. Endocr Relat Cancer. (2010) 17:361–71. doi: 10.1677/ERC-10-0018

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Qian ZR, Sano T, Yoshimoto K, Asa SL, Yamada S, Mizusawa N, et al. Tumor-specific downregulation and methylation of the CDH13 (H-cadherin) and CDH1 (E-cadherin) genes correlate with aggressiveness of human pituitary adenomas. Modern Pathol. (2007) 20:1269. doi: 10.1038/modpathol.3800965

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Butz H, Likó I, Czirják S, Igaz P, Korbonits M, Rácz K, et al. MicroRNA profile indicates downregulation of the TGFβ pathway in sporadic non-functioning pituitary adenomas. Pituitary. (2010) 14:112–24. doi: 10.1007/s11102-010-0268-x

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li Z, Li C, Liu C, Yu S, Zhang Y. Expression of the long non-coding RNAs MEG3, HOTAIR, and MALAT-1 in non-functioning pituitary adenomas and their relationship to tumor behavior. Pituitary. (2014) 18:42–7. doi: 10.1007/s11102-014-0554-0

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Gov E, Arga KY. Differential co-expression analysis reveals a novel prognostic gene module in ovarian cancer. Sci Rep. (2017) 7:4996. doi: 10.1038/s41598-017-05298-w

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zhang W, Zang Z, Song Y, Yang H, Yin Q. Co-expression network analysis of differentially expressed genes associated with metastasis in prolactin pituitary tumors. Mol Med Rep. (2014) 10:113–8. doi: 10.3892/mmr.2014.2152

PubMed Abstract | CrossRef Full Text | Google Scholar

13. de Bruin C, Pereira A, Feelders R, Romijn J, Roelfsema F, Sprij-Mooij D, et al. Coexpression of dopamine and somatostatin receptor subtypes in corticotroph adenomas. Endocrinology. (2009) 150:1558–9. doi: 10.1210/endo.150.3.9998

CrossRef Full Text | Google Scholar

14. Kolesnikov N, Hastings E, Keays M, Melnichuk O, Tang Y, Williams E, et al. ArrayExpress update—simplifying data submissions. Nucleic Acids Res. (2014) 43:D1113–6. doi: 10.1093/nar/gku1057

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Barrett T, Wilhite S, Ledoux P, Evangelista C, Kim I, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets—update. Nucleic Acids Res. (2012) 41:D991–5. doi: 10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Hernández-Ramírez LC, Morgan RM, Barry S, D'Acquisto F, Prodromou C, Korbonits M. Multi-chaperone function modulation and association with cytoskeletal proteins are key features of the function of AIP in the pituitary gland. Oncotarget. (2018) 9:9177. doi: 10.18632/oncotarget.24183

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Xing W, Qi Z, Huang C, Zhang N, Zhang W, Li Y, et al. Genome-wide identification of lncRNAs and mRNAs differentially expressed in non-functioning pituitary adenoma and construction of an lncRNA-mRNA co-expression network. Biol Open. (2019) 8:bio037127. doi: 10.1242/bio.037127

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Hardy J. Functional microsurgical treatment of hypersecreting pituitary disorders. Stereotact Funct Neurosurg. (1978) 41:217–8. doi: 10.1159/000102419

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Smyth GK. Limma: linear models for microarray data. In: Gentleman R, Carey V, Dudoit S, Irizarry R, Huber W, editors. Bioinformatics and Computational Biology Solutions Using R and Bioconductor. New York, NY: Springer (2005). p. 397–420. doi: 10.1007/0-387-29362-0_23

CrossRef Full Text | Google Scholar

20. Kamburov A, Stelzl U, Lehrach H, Herwig R. The ConsensusPathDB interaction database: 2013 update. Nucleic Acids Res. (2012) 41:D793–800. doi: 10.1093/nar/gks1055

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. (2016) 45:D353–61. doi: 10.1093/nar/gkw1092

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. (2000) 25:25. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. (2003) 4:2. doi: 10.1186/1471-2105-4-2

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Assenov Y, Ramírez F, Schelhorn S, Lengauer T, Albrecht M. Computing topological parameters of biological networks. Bioinformatics. (2007) 24:282–4. doi: 10.1093/bioinformatics/btm554

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Chin C, Chen S, Wu H, Ho C, Ko M, Lin C. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. (2014) 8:S11. doi: 10.1186/1752-0509-8-S4-S11

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Aguirre-Gamboa R, Gomez-Rueda H, Martínez-Ledesma E, Martínez-Torteya A, Chacolla-Huaringa R, Rodriguez-Barrientos A, et al. SurvExpress: an online biomarker validation tool and database for cancer gene expression data using survival analysis. PLoS ONE. (2013) 8:e74250. doi: 10.1371/journal.pone.0074250

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Weinstein JN, Collisson EA, Mills GB, Shaw KR, Ozenberger BA, Ellrott K, et al. The cancer genome atlas pan-cancer analysis project. Nat Genet. (2013) 45:1113. doi: 10.1038/ng.2764

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Lee Y, Liu J, Patel S, Cloughesy T, Lai A, Farooqi H, et al. Genomic landscape of meningiomas. Brain Pathol. (2010) 20:751–62. doi: 10.1111/j.1750-3639.2009.00356.x

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Robinson G, Parker M, Kranenburg TA, Lu C, Chen X, Ding L, et al. Novel mutations target distinct subgroups of medulloblastoma. Nature. (2012) 488 43–8. doi: 10.1038/nature11213

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Lonsdale J, Thomas J, Salvatore M, Phillips R, Lo E, Shad S, et al. The genotype-tissue expression (GTEx) project. Nat Genet. (2013) 45:580. doi: 10.1038/ng.2653

CrossRef Full Text | Google Scholar

32. Thapar K, Kovacs K, Scheithauer BW, Stefaneanu L, Horvath E, Peter J P, et al. Proliferative activity and invasiveness among pituitary adenomas and carcinomas: an analysis using the MIB-1 antibody. Neurosurgery. (1996) 38:99–107. doi: 10.1097/00006123-199601000-00024

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Thompson JA. Molecular cloning and expression of carcinoembryonic antigen gene family members. Tumor Biol. (1995) 16:10–6. doi: 10.1159/000217923

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Kuroki M, Matsuo Y, Kinugasa T, Matsuoka Y. Augmented expression and release of nonspecific cross-reacting antigens (NCAs), members of the CEA family, by human neutrophils during cell activation. J Leukoc Biol. (1992) 52:551–7. doi: 10.1002/jlb.52.5.551

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Jantscheff P, Terracciano L, Lowy A, Glatz-Krieger K, Grunert F, Micheel B, et al. Expression of CEACAM6 in resectable colorectal cancer: a factor of independent prognostic significance. J Clin Oncol. (2003) 21:3638–46. doi: 10.1200/JCO.2003.55.135

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Duxbury MS, Matros E, Clancy T, Bailey G, Doff M, Zinner MJ, et al. CEACAM6 is a novel biomarker in pancreatic adenocarcinoma and PanIN lesions. Ann Surg. (2005) 241:491–6. doi: 10.1097/01.sla.0000154455.86404.e9

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Baer BR, Rettie AE. CYP4B1: an enigmatic P450 at the interface between xenobiotic and endobiotic metabolism. Drug Metabol Rev. (2006) 38:451–76. doi: 10.1080/03602530600688503

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Imaoka S, Yoneda Y, Sugimoto T, Ikemoto S, Hiroi T, Yamamoto K, et al. Androgen regulation of CYP4B1 responsible for mutagenic activation of bladder carcinogens in the rat bladder: detection of CYP4B1 mRNA by competitive reverse transcription-polymerase chain reaction. Cancer Lett. (2001) 166:119–23. doi: 10.1016/S0304-3835(00)00572-3

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Heaney JD, Michelson MV, Youngren KK, Lam M-YJ, Nadeau JH. Deletion of eIF2beta suppresses testicular cancer incidence and causes recessive lethality in agouti-yellow mice. Hum Mol Genet. (2009) 18:1395–404. doi: 10.1093/hmg/ddp045

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang L, Zhan Y, Song E, Yu Y, Jiu Y, Du W, et al. HID-1 is a peripheral membrane protein primarily associated with the medial- and trans- Golgi apparatus. Protein Cell. (2011) 2:74–85. doi: 10.1007/s13238-011-1008-3

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Quinlan R, Hutchison C, Lane B. Intermediate filament proteins. Protein Profile. (1994) 1:779–911.

PubMed Abstract | Google Scholar

42. Steinert P. Molecular and cellular biology of intermediate filaments. Annu Rev Biochem. (1988) 57:593–625. doi: 10.1146/annurev.biochem.57.1.593

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Hsu RM, Hsieh YJ, Yang TH, Chiang YC, Kan CY, Lin YT, et al. Binding of the extreme carboxyl-terminus of PAK-interacting exchange factor β (βPIX) to myosin 18A (MYO18A) is required for epithelial cell migration. BBA- Mol Cell Res. (2014) 1843:2513–27. doi: 10.1016/j.bbamcr.2014.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhang J, Wei W, Jin H-C, Ying R-C, Zhu A-K, Zhang F-J. Programmed cell death 2 protein induces gastric cancer cell growth arrest at the early S phase of the cell cycle and apoptosis in a p53-dependent manner. Oncol Rep. (2014) 33:103–10. doi: 10.3892/or.2014.3551

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Shimada A, Yamaguchi A, Kohda D. Structural basis for the recognition of two consecutive mutually interacting DPF motifs by the SGIP1 μ homology domain. Sci Rep. (2016) 6:19565. doi: 10.1038/srep19565

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Chwedorowicz R, Raszewski G, Kapka-Skrzypczak L, Sawicki K, Studzinski T. Event-related potentials (ERP) and SGIP1 gene polymorphisms in alcoholics: relation to family history of alcoholism and drug usage. Ann Agric Environ Med. (2016) 23:618–24. doi: 10.5604/12321966.1226856

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Liu T, Wan L, Wu Y, Chen J, Huang J. hSWS1·SWSAP1 is an evolutionarily conserved complex required for efficient homologous recombination repair. J Biol Chem. (2011) 286:41758–66. doi: 10.1074/jbc.M111.271080

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Trouillas J, Roy P, Sturm N, Dantony E, Cortet-Rudelli C, Viennet G, et al. A new prognostic clinicopathological classification of pituitary adenomas: a multicentric case–control study of 410 patients with 8 years post-operative follow-up. Acta Neuropathol. (2013) 26:123–35. doi: 10.1007/s00401-013-1084-y

CrossRef Full Text | Google Scholar

49. Raverot G, Dantony E, Beauvy J, Vasiljevic A, Mikolasek S, Borson-Chazot F, et al. Risk of recurrence in pituitary neuroendocrine tumors: a prospective study using a five-tiered classification. J Clin Endocrinol Metab. (2017) 102:3368–74. doi: 10.1210/jc.2017-00773

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Suojun Z, Feng W, Dongsheng G, Ting L. Targeting Raf/MEK/ERK pathway in pituitary adenomas. Eur J Cancer. (2012) 48:389–95. doi: 10.1016/j.ejca.2011.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

51. O'Hayre M, Degese MS, Gutkind JS. Novel insights into G protein and G protein-coupled receptor signaling in cancer. Curr Opin Cell Biol. (2014) 27:126–35. doi: 10.1016/j.ceb.2014.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Musat M. Enhanced protein kinase B/Akt signalling in pituitary tumours. Endocr Relat Cancer. (2005) 12:423–33. doi: 10.1677/erc.1.00949

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Kuroki M, Matsushita H, Matsumoto H, Hirose Y, Senba T, Yamamoto T. Nonspecific cross-reacting antigen-50/90 (NCA-50/90) as a new tumor marker. Anticancer Res. (1999) 19:5599–606.

PubMed Abstract | Google Scholar

54. Blumenthal RD, Leon E, Hansen HJ, Goldenberg DM. Expression patterns of CEACAM5 and CEACAM6 in primary and metastatic cancers. BMC Cancer. (2007) 7:2. doi: 10.1186/1471-2407-7-2

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Rizeq B, Zakaria Z, Ouhtit A. Towards understanding the mechanisms of actions of carcinoembryonic antigen-related cell adhesion molecule 6 in cancer progression. Cancer Sci. (2018) 109:33–42. doi: 10.1111/cas.13437

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Chen J, Li Q, An Y, Lv N, Xue X, Wei J, et al. CEACAM6 induces epithelial-mesenchymal transition and mediates invasion and metastasis in pancreatic cancer. Int J Oncol. (2013) 43:877–85. doi: 10.3892/ijo.2013.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Karagoz K, Lehman HL, Stairs DB, Sinha R, Arga KY. Proteomic and metabolic signatures of esophageal squamous cell carcinoma. Curr Cancer Drug Tar. (2016) 16:721–36. doi: 10.2174/1568009616666160203113721

CrossRef Full Text | Google Scholar

58. Iscan M, Klaavuniemi T, Çoban T, Kapucuoglu N, Pelkonen O, Raunio H. The expression of cytochrome P450 enzymes in human breast tumours and normal breast tissue. Breast Cancer Res Treat. (2001) 70:47–54. doi: 10.1023/A:1012526406741

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Wang Q, Tang Y, Yu H, Yin Q, Li M, Shi L, et al. CCL18 from tumor-cells promotes epithelial ovarian cancer metastasis via mTOR signaling pathway. Mol Carcinog. (2015) 55:1688–99. doi: 10.1002/mc.22419

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Long NP, Jung KH, Yoon SJ, Anh NH, Nghi TD, Kang YP, et al. Systematic assessment of cervical cancer initiation and progression uncovers genetic panels for deep learning-based early diagnosis and proposes novel diagnostic and prognostic biomarkers. Oncotarget. (2017) 8:22689. doi: 10.18632/oncotarget.22689

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Siddiq A, Couch FJ, Chen GK, Lindström S, Eccles D, Millikan RC, et al. A meta-analysis of genome-wide association studies of breast cancer identifies two novel susceptibility loci at 6q14 and 20q11. Hum Mol Genet. (2012) 21:5373–84. doi: 10.1093/hmg/dds381

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Harada H, Nagai H, Tsuneizumi M, Mikami I, Sugano S, Emi M. Identification of DMC1, a novel gene in the TOC region on 17q25.1 that shows loss of expression in multiple human cancers. J Hum Genet. (2001) 46:90–5. doi: 10.1007/s100380170115

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Campan M, Moffitt M, Houshdaran S, Shen H, Widschwendter M, Daxenbichler G, et al. Genome-scale screen for DNA methylation-based detection markers for ovarian cancer. PLoS ONE. (2011) 6:e0028141. doi: 10.1371/journal.pone.0028141

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Liu Y, Nan F, Lu K, Wang Y, Liu Y, Wei S, et al. Identification of key genes in endometrioid endometrial adenocarcinoma via TCGA database. Cancer Biomarkers. (2017) 21:11–21. doi: 10.3233/CBM-170164

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Makowska KA, Hughes RE, White KJ, Wells CM, Peckham M. Specific myosins control actin organization, cell morphology, and migration in prostate cancer cells. Cell Rep. (2015) 13:2118–25. doi: 10.1016/j.celrep.2015.11.012

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Ussowicz M, Jaśkowiec A, Meyer C, Marschalek R, Chybicka A, Szczepanski T, et al. A three-way translocation of MLL, MLLT11, and the novel reciprocal partner gene MYO18A in a child with acute myeloid leukemia. Cancer Genet. (2012) 205:261–5. doi: 10.1016/j.cancergen.2012.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Barboza N, Minakhina S, Medina DJ, Balsara B, Greenwood S, Huzzy L, et al. PDCD2 functions in cancer cell proliferation and predicts relapsed leukemia. Cancer Biol Ther. (2013) 14:546–55. doi: 10.4161/cbt.24484

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Wang W, Song XW, Bu XM, Zhang N, Zhao C-H. PDCD2 and NCoR1 as putative tumor suppressors in gastric gastrointestinal stromal tumors. Cell Oncol. (2015) 39:129–37. doi: 10.1007/s13402-015-0258-0

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Pitkänen E, Cajuso T, Katainen R, Kaasinen E, Välimäki N, Palin K, et al. Frequent L1 retrotranspositions originating from TTC28 in colorectal cancer. Oncotarget. (2014) 5:1781. doi: 10.18632/oncotarget.1781

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Trevaskis J, Walder K, Foletta V, Kerr-Bayles L, Mcmillan J, Cooper A, et al. Src homology 3-domain growth factor receptor-bound 2-like (Endophilin) interacting protein 1, a novel neuronal protein that regulates energy balance. Endocrinology. (2005) 146:3757–64. doi: 10.1210/en.2005-0282

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Li JN, Zhao L, Wu J, Wu B, Yang H, Zhang HH, et al. Differences in gene expression profiles and carcinogenesis pathways between colon and rectal cancer. J Dig Dis. (2011) 13:24–32. doi: 10.1111/j.1751-2980.2011.00551.x

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Bae J-H, Kang MJ, Yang K-M, Kim T-O, Yi JM. Epigenetically silenced microRNAs in gastric cancer: functional analysis and identification of their target genes. Oncol Rep. (2015) 34:1017–26. doi: 10.3892/or.2015.4036

PubMed Abstract | CrossRef Full Text | Google Scholar