Identification of metastasis-associated genes in colorectal cancer using metaDE and survival analysis.

The aim of the present study was to detect the candidate genes involved in the metastasis of colorectal cancer (CRC). Gene expression profiles of primary and metastatic CRC samples in the GSE14297 and GSE49355 datasets were downloaded from the Gene Expression Omnibus database. Subsequent to processing, Fishers exact test and the metaDE package in R language were applied to screen the differentially expressed genes (DEGs) between primary and metastatic CRC samples. In addition, function and pathway enrichment analysis was performed using online tools in the Database for Annotation, Visualization, and Integrated Discovery resource and common DEGs in GSE14297 and GSE49355 were identified. Their expression values in another dataset, GSE29621, were then collected in order to screen the genes with high standard deviations between primary and metastatic samples, which were considered as candidate metastasis-associated genes. Candidate genes were finally verified by performing survival analysis via the log-rank test. A total of 370 DEGs were screened in GSE14297 and GSE49355, and 77 common DEGs were identified. Upregulated DEGs were mainly enriched in the immune, energy metabolism and drug metabolism-associated functions. Downregulated DEGs were mainly enriched in cell adhesion-associated functions. A total of 12 genes, including the carbonic anhydrase II (CA2), carcinoembryonic antigen-related cell adhesion molecule 7 (CEACAM7), Fc fragment of immunoglobulin G binding protein (FCGBP), and placenta-specific 8 (PLAC8), were the candidate metastasis-associated genes, among which FCGBP expression significantly decreased the overall survival time of patients. The selected candidate metastasis-associated gene, FCGBP, may be used as a potential therapeutic target in patients with metastatic CRC.

Introduction

Colorectal cancer (CRC) is a type of malignant tumor developed in the rectum or colon, and it is estimated that there are currently >1 million people who are living with CRC in the USA (1). It is the second leading cause of cancer-associated mortalities among all cancer types that affect both females and males (2). An association has been identified between CRC and a relatively high 5-year survival rate at early stages of the disease, while the survival is poor when distant and nodal metastasis is exhibited (3,4). The number of patients with metastatic status is up to 30%, and these patients present a 5-year survival rate of <10% (5,6). In China, the mortality associated with CRC increased by 17.9% between 1973–1975 and 2004–2005 (7). Thus, it is of great importance to identify the genetic changes and potential underlying mechanisms in CRC.

Genes involved in the metastasis process in CRC have been extensively studied. Kirsten rat sarcoma viral oncogene homolog (KRAS) mutations have been reported to be a predictor of resistance to the monoclonal antibodies targeting the epidermal growth factor receptor (EGFR) in metastatic CRC (8). In addition, B-Raf proto-oncogene, serine/threonine kinase (BRAF), and phosphatase and tensin homolog (PTEN) mutations are considered to be implications for the targeted treatment in metastatic CRC (9). Furthermore, a previous study revealed that the forced expression of chemokine (C-X-C motif) receptor 3 (CXCR3) promoted preferential metastasis to draining lymph nodes in CRC, while CXCR3 knockdown significantly decreased dissemination of cancer cells to the lungs and liver; thus, CXCR3 may be a potential therapy against metastatic CRC (10). The receptor of advanced glycation end products is a prognostic biomarker of CRC metastasis (11).

Furthermore, numerous metastasis-associated genes have been screened by gene expression profiling. For instance, Stange et al (12) used microarray data to identify genes associated with CRC metastasis to the liver. As a result, 163 unique genes were identified to be significantly overexpressed, whereas 15 genes were significantly downregulated (12). These genes, including CYP2E1, CYP4A11, CRP, ORM1, SAA1, APOA1, FGA and FGB, may be associated with metabolism and inflammatory response. Del Rio et al (13) identified the 33-gene signature to classify the hepatic metastases, primary tumors, normal colon mucosa and normal liver tissues, and indicated that these genes may influence the CRC metastasis to the liver by involving the extracellular matrix remodeling. However, the genes identified in a single dataset may be limited if they have not been confirmed in other datasets.

The present study aimed to use the data of Stange et al (12) and Del Rio et al (13) together, to further detect the candidate metastasis-associated genes in CRC. The metaDE package in R language, which implements 12 major meta-analyses in differential expression screening (14), was used to screen the differentially expressed genes (DEGs) between primary and metastatic cancer samples in the two datasets. Functional enrichment was also conducted for the significantly associated functions and pathways. By calculating the standard deviations of common DEGs in another dataset with clinical data, candidate metastasis-associated genes were collected, followed by survival analysis.

Materials and methods

Microarray data

The Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database was retrieved for obtaining the microarray data with the following accession numbers: GSE14297 (12), GSE49355 (13,15) and GSE29621 (16). A total of 18 primary CRC and 18 matched liver metastasis samples were available in the GSE14297 dataset, based on the GPL6370 Illumina human-6 v2.0 expression beadchip (extended). Similarly, the expression data of 20 primary CRC and 19 matched liver metastasis samples in the GSE49335 dataset, based on GPL10430 Rosenstiel Fundulus heteroclitus 7K array, were downloaded. In addition, the GSE29621 dataset included 46 primary CRC samples, 18 metastatic samples and 1 sample with unknown metastatic status, based on the platform of GPL570 (HG_U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array.

Data processing

The gene expression profiles of GSE14297 were prenormalized, while the raw CEL data (media.affymetrix.com/support/developer/powertools/changelog/gcos-agcc/cel.html) in GSE49335 and GSE29621 were initially normalized using the robust multi-array average method in R software, version 2.6.0 (R Foundation for Statistical Computing, Vienna, Austria; www.R-project.org/). The median value of multiple probes corresponding to a same gene was used as the expression value.

Microarray meta-analysis for DEGs

Initially, the expression values obtained from the GSE14297 and GSE49335 datasets were used to screen DEGs between primary and metastatic CRC samples. To eliminate discrepancies, the metaDE package in R language (14) in R was used as a supplement for Fishers exact test. Next, clustering analysis was performed to detect the distinguishing effect of metaDE on differential expression in different sample groups. The threshold for DEGs was a false discovery rate (FDR) of <0.05 in metaDE, combined with a P value of <0.05 in Fisher's exact test. P<0.05 and FDR<0.05 were considered to indicate a statistically significant difference.

DEG enrichment analysis

In order to identify the significantly altered functions and pathways during the metastasis of CRC, the online tools in Database for Annotation, Visualization, and Integrated Discovery (DAVID version 6.7; http://david.abcc.ncifcrf.gov/) were used for the enrichment analysis. FDR<0.01 was set as the cut-off value for the enrichment process.

Metastasis-associated genes in CRC

Fold-changes (increase or decreased) in the expression of the selected genes were investigated, and DEGs with a fold-change of >2 were considered to be the metastasis-associated genes in CRC. Common metastasis-associated genes were obtained by comparing data in the profiles of GSE14297 and GSE49335. Next, the expression values of these common metastasis-associated genes recorded in GSE29621 were obtained, in order to calculate the standard deviation of common genes in all the 65 samples. Genes with a high standard deviation in their expression were defined as the candidate metastatic genes in CRC.

Survival analysis

The clinical data of CRC patients were also obtained from the GSE29621 dataset, and were subjected to survival analysis. Using the log-rank test, the effects of metastasis and adjuvant chemotherapy on survival were detected. Additionally, all the 65 samples were divided into the high expression and low expression groups, according to the expression levels of the candidate metastatic genes. Subsequently, the log-rank test was conducted to identify the effect of these candidate genes on the survival of CRC patients.

Results

Data processing and DEG screening

Initially, the expression profiles were normalized, and the results are shown in Fig. 1. Using the threshold of FDR<0.05 in metaDE and P<0.05 in Fisher's exact test, a total of 370 genes were identified to be differentially expressed between the primary and metastatic CRC samples in GSE14297 and GSE49335, among which 153 upregulated and 217 downregulated genes were detected. Fig. 2 shows the number of genes that are differentially expressed against various FDR values. The clustering analysis results indicated that the DEGs selected using metaDE analysis may be able to distinguish metastatic cancer samples from non-metastatic samples (Fig. 3).

Figure 1.

Box plots for the (A) GSE14927 and (B) GSE49355 datasets following normalization. The x-axis represents the samples, and the y-axis represents the distribution of expression levels.

Figure 2.

Plots of numbers of DEGs against FDR obtained from the metaDE analysis. DEGs, differentially expressed genes; FDR, false discovery rate.

Figure 3.

Hierarchical clustering of DEGs identified using metaDE analysis in primary and metastatic CRC samples. Values 0 and 1 at the top of the graph represent the primary and metastatic samples, respectively. Dataset 1 refers to the profiles of GSE14297, and database 2 refers to the profiles of GSE9355. Green indicates the downregulated genes, while red indicates the upregulated genes in metastatic samples. DEGs, differentially expressed genes; CRC, colorectal cancer.

DAVID was applied to investigate the function and pathway enrichment. Upregulated genes were mainly enriched in immune-associated functions, such as repose to wounding and inflammatory response, as well as in metabolic process-associated functions, including peptidase inhibitor activity and carbohydrate binding. Complement and coagulation cascades, and drug metabolism were among the significantly enriched pathways (Table I). Downregulated genes were significantly enriched in the functions of vasculature development, cell adhesion, biological adhesion and extracellular matrix structure; however, no pathways were enriched (Table II).

Table I.

Enriched functions and pathways of upregulated genes in the GSE14927 and GSE40355 datasets.

Category	Term	Count	FDR
GOTERM_BP_FAT	GO:0009611 - response to wounding	62	6.60×10−35
GOTERM_BP_FAT	GO:0002526 - acute inflammatory response	31	1.14×10−28
GOTERM_BP_FAT	GO:0006954 - inflammatory response	40	1.52×10−21
GOTERM_BP_FAT	GO:0006952 - defense response	47	2.73×10−17
GOTERM_BP_FAT	GO:0006956 - complement activation	17	5.18×10−16
GOTERM_CC_FAT	GO:0005615 - extracellular space	72	3.85×10−38
GOTERM_CC_FAT	GO:0005576 - extracellular region	106	1.02×10−32
GOTERM_CC_FAT	GO:0044421 - extracellular region part	73	2.28×10−29
GOTERM_CC_FAT	GO:0034364 - high-density lipoprotein particle	13	5.44×10−13
GOTERM_CC_FAT	GO:0032994 - protein-lipid complex	14	1.92×10−12
GOTERM_MF_FAT	GO:0004866 - endopeptidase inhibitor activity	25	1.14×10−15
GOTERM_MF_FAT	GO:0030414 - peptidase inhibitor activity	25	4.15×10−15
GOTERM_MF_FAT	GO:0004857 - enzyme inhibitor activity	30	4.14×10−14
GOTERM_MF_FAT	GO:0004867 - serine-type endopeptidase inhibitor activity	18	2.68×10−11
GOTERM_MF_FAT	GO:0030246 - carbohydrate binding	26	1.05×10−7
KEGG_PATHWAY	hsa04610: Complement and coagulation cascades	31	1.22×10−28
KEGG_PATHWAY	hsa00982: Drug metabolism	15	1.31×10−7
KEGG_PATHWAY	hsa00980: Metabolism of xenobiotics by cytochrome P450	13	1.32×10−5
KEGG_PATHWAY	hsa00830: Retinol metabolism	11	4.99×10−4

GO, gene oncology; BP, biological process; CC, cellular component; MF, molecule function; FDR, false discovery rate.

Table II.

Enriched functions of downregulated genes in datasets of GSE14927 and GSE40355.

Category	Term	Count	FDR
GOTERM_BP_FAT	GO:0001944 - vasculature development	17	3.30×10−6
GOTERM_BP_FAT	GO:0007155 - cell adhesion	26	1.42×10−5
GOTERM_BP_FAT	GO:0022610 - biological adhesion	26	1.46×10−5
GOTERM_BP_FAT	GO:0001568 - blood vessel development	15	1.42×10−4
GOTERM_BP_FAT	GO:0048514 - blood vessel morphogenesis	13	1.31×10−3
GOTERM_BP_FAT	GO:0042127 - regulation of cell proliferation	23	7.84×10−3
GOTERM_CC_FAT	GO:0044421 - extracellular region part	53	2.54×10−22
GOTERM_CC_FAT	GO:0005576 - extracellular region	72	1.22×10−21
GOTERM_CC_FAT	GO:0031012 - extracellular matrix	29	8.38×10−15
GOTERM_CC_FAT	GO:0005578 - proteinaceous extracellular matrix	28	1.31×10−14
GOTERM_CC_FAT	GO:0005615 - extracellular space	31	8.12×10−9
GOTERM_CC_FAT	GO:0044420 - extracellular matrix part	11	3.75×10−4
GOTERM_CC_FAT	GO:0030934 - anchoring collagen	5	1.56×10−3
GOTERM_CC_FAT	GO:0005581 - collagen	7	1.59×10−3
GOTERM_MF_FAT	GO:0005201 - extracellular matrix structural constituent	11	1.11×10−5

GO, gene oncology; BP, biological process; CC, cellular component; MF, molecule function; FDR, false discovery rate.

Screening of candidate metastasis-associated genes

Following the intersection of DEGs in the GSE49355 and GSE14927 datasets, a total of 77 common DEGs were obtained, consisting of 17 downregulated and 59 upregulated genes. The expression values of these 77 common DEGs in the GSE29621 dataset were used to calculate the standard deviation, and 12 genes with a standard deviation >1 were identified: carbonic anhydrase II (CA2); carcinoembryonic antigen-related cell adhesion molecule 7 (CEACAM7); chloride channel accessory 4 (CLCA4); CLCA1; CXC ligand 14 (CXCL14); Fc fragment of immunoglobulin G binding protein (FCGBP); immunoglobulin J polypeptide, linker protein for immunoglobulin α and µ polypeptides (IGJ); lipocalin 2 (LCN2); matrix metallopeptidase 1 (MMP1); MMP3; peptidase inhibitor 3, skin-derived (PI3); and placenta-specific 8 (PLAC8). A larger deviation value indicated great changes in the samples; thus, these 12 genes were regarded as candidate metastasis-associated genes, and may exert an intensive effect on the survival of patients.

Univariate survival analysis was conducted via the log-rank test. As shown in Fig. 4, a significantly decreased overall survival time was observed in patients with metastasis (P<0.05), while the disease-free survival times were not changed according to the metastatic status of the cancer. Furthermore, adjuvant chemotherapy demonstrated no statistically significant effect on the overall and disease-free survival times (P>0.05). Analysis of the 12 candidate metastasis-associated genes revealed that downregulation of FCGBP markedly decreased the patients' overall survival time (P=0.184); however, no effect was observed on the disease-free survival time (Fig. 5).

Figure 4.

Survival curves demonstrating the effects of metastasis and adjuvant chemotherapy on survival time in colorectal cancer patients. chemo, chemotherapy.

Figure 5.

Survival curve of the effect of the FCGBP gene on the survival time of colorectal cancer patients. FCGBP, Fc fragment of immunoglobulin G binding protein.

Discussion

At present, metastasis remains an important factor contributing towards the majority of cancer-associated mortalities, as metastatic cancer is often resistant to conventional therapies (17,18). However, reliable and sensitive methods for detecting early metastasis in CRC are not currently available. In the present study, gene expression data were initially analyzed using the metaDE package in R language, in order to identify the DEGs in CRC samples, and a total of 370 genes were identified as DEGs between primary and metastatic cancer samples. In addition, clustering analysis demonstrated the reliability of the identified DEGs. By comparing the expression levels in different datasets, 12 genes were identified among the DEGs as the metastasis-associated genes in CRC, including CA2, CEACAM7, CLCA4, CLCA1, CXCL14, FCGBP, IGJ, LCN2, MMP1, MMP3, PI3 and PLAC8. Subsequently, FCGBP was demonstrated to influence the patients' survival.

The carbonic anhydrase family is a set of proteins that are important to the pH regulatory system. Carbonic anhydrase IX (CA9) is overexpressed in numerous types of cancer cells, and various agents targeting CA9 are currently in preclinical or clinical development for cancer therapy (19,20). CA9 is also a therapeutic target for metastasis (21). CA2 and CA9 were demonstrated to be overexpressed in endometrial adenocarcinoma, regulating the pH in the tumor microenvironment (22). Similarly, CA2 was identified as a metastasis-associated gene in CRC in the present study. Proteomics analysis identified CA2 as a potential biomarker, which exerted significant inhibitory effects on cell growth in CRC (23). CA2 expression was also associated with lymph node metastasis in gastric cancer, and a previous study reported that its reduction may contribute to tumor metastasis (24). Therefore, CA2 may also be used as a therapeutical target for metastasis.

CEACAM7 is expressed in normal colon tissues, but is downregulated in CRC (25). By contrast, its expression is upregulated in liver metastases (26). CEACAM6 also belongs to the CEACAM subgroup, which is mainly associated with the cell membrane. A previous study suggested that CEACAM6 plays a role in tumor cell migration, adhesion and invasion, as well as the formation of distant metastases (27). In addition, the expression of CEACAM6 in CRC is an independent prognostic factor (28). Furthermore, downregulation of CEACAM7 in early-stage adenomas represents certain observable molecular events that may lead to CRC (29). In the present study, CEACAM7 was also identified as a metastasis-associated gene.

FCGBP has been reported to be implicated in ulcerative colitis, which is a chronic inflammatory disease predisposing to CRC (30,31). This gene may play a role in anti-inflammation and cell protection in tissues (32). The mutation of FCGBP has been previously identified in CRC (33), while it was also identified to be a DEG in metastatic CRC in the current study. Notably, in gallbladder adenocarcinoma, FCGBP functions as an important maker for clinical prognosis (34). Using the univariate Cox model and Kaplan-Meier survival curve, a correlation between FCGBP expression and overall survival in ovarian adenocarcinoma was observed (35). The survival curve of CRC patients with liver metastasis in the present study also detected an association of FCGBP with prognosis. Furthermore, FCGBP was significantly enriched in the function of cell adhesion. Metastasis is facilitated by the cell-cell interactions between the endothelium and tumor cells in distant tissues. Cell adhesion occurring in vasculature of specific organs is an essential step in cancer metastasis (36). Thus, we infer that FCGBP may exert an important role during metastasis in CRC by participating in cell adhesion.

In conclusion, using the metaDE package in R language, DEGs were screened in two sets of gene expression profiles, which may be able to distinguish the metastasis samples from the primary cancer samples. Candidate metastasis-associated genes identified by comparing with a third expression dataset were further subjected to survival analysis. Therefore, CA2 has the potential to be used as a therapeutic target for metastatic CRC, while FCGBP may affect patients' survival by participating in cell adhesion.