Candidate genes for predicting the survival of patients with gastric cancer: a study based on The Cancer Genome Atlas (TCGA) database.

BACKGROUND: Gastric cancer (GC) is the second most frequent cause of cancer-related mortality in the world, and the five-year survival rate for GC remains very low universally. In recent years, it has become a consensus that genetic changes are associated with carcinogenesis of GC, and precision medicine based on genetic changes is one of the most popular treatments for GC patients. However, the association between some genes and GC-related protein signaling pathways is still not well understood. This study revealed that seven genes were closely related to the survival probability in GC patients.
METHODS: We downloaded the gene expression data of GC patients from The Cancer Genome Atlas (TCGA) databases, and integrated bioinformatic analysis was performed, such as differential gene expression analysis, including Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathways analyses, as well as survival analysis. The r package "survival" was used to analyze the Kaplan-Meier survival analysis, which showed the associations between specific gene expressions and the outcomes of patients with GC to identify which genes could be potential prognostic biomarkers.
RESULTS: This study revealed that seven genes: alcohol dehydrogenase 4 (ADH4), histamine receptor H3 (HRH3), neuropeptide Y2 receptor (NPY2R), apolipoprotein AI (APOA1), N-acetylgalactosaminyltransferase 14 (GALNT14), leucine-rich repeats and IQ motif containing 1 (LRRIQ1), and coiled-coil-domain-containing 57 (CCDC57). These seven genes were closely related to the survival probability of GC patients (P<0.05).
CONCLUSIONS: Our study found seven genes which could be considered as candidate prognostic biomarkers and therapeutic targets. 2020 Translational Cancer Research. All rights reserved.

Introduction

Gastric cancer (GC) is the second most frequent cause of cancer-related mortality in the world. Despite the developments in endoscopic technology, the great progression made in early cancer screening, and the achievements made in relation to Helicobacter pylori eradication, the 5-year survival rate for GC remains very low worldwide (1).

Dynamic changes in the genome play an essential role in the progress of carcinogenesis (2). The Cancer Genome Atlas (TCGA) provides a comprehensive overview of gene expression, RNA-seq, DNA copy-number, somatic mutations, and DNA methylation profiles in tumors, as well as providing the matched clinical information of patients with cancer (3). This publicly available cancer genomics data set allows for improved diagnostic methods, treatment criteria and, ultimately, cancer prevention (4).

Many studies have proved that, in GC patients, the TNM stage is not the only factor impacting survival (5); gene expression also bears a strong association. Previous studies have revealed that the overexpression of tumor protein 53 (p53) and Mucin 1 (MUC1), and the decrease of expressions of phosphatase and tension homolog gene (PTEN), E-cadherin gene, and SMAD family member 4 (SMAD4), were found to be associated with poor prognosis of GC patients (6). Recent studies have also found that people with high expression of LncRNA AL139147 show a tendency towards poor prognosis (7). Competing endogenous RNAs (ceRNA) analysis has also shown that the complex mechanisms of the ceRNA network are essential in the progression of GC (8). Various genes that could be considered as candidate prognostic biomarkers and therapeutic targets are yet to be revealed and comprehensively understood.

In this study, we obtained the gene expression profiles of 375 gastric tumors and 32 adjacent non-tumor samples from TCGA database. A Gene expression matrix was obtained, and R package “edgeR” was used to examine differentially expressed genes (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) (9). The gene expression profiles were combined with clinical survival information. Integrated bioinformatic analyses were performed using “R”, including differential gene expression analysis, The Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis (10,11), as well as survival analysis. To identify which gene could predict the outcomes of GC patients, Kaplan-Meier survival analysis was performed.

Methods

Gene expression profile and patient clinical data

We downloaded the gene expression profiles and clinical data of GC patients from TCGA database (https://cancergenome.nih.gov/) in November 2018 and analyzed the statistics between December 2018 and May 2019. The exclusion criteria were as follows: (I) without clinical information or prognostic statistics like survival time; (II) without matched adjacent non-tumor tissues; and (III) not stomach adenocarcinoma. Ultimately, 407 samples including 375 GC tumor tissues and 32 adjacent non-tumor samples, were collected for integrated bioinformatics analysis. There was no need for ethical approval as all data in this study were downloaded from public databases (TCGA), and the data processing met the TCGA publication guidelines (https://cancergenome.nih.gov/publications/guidelines).

Differential gene expression in GC

The gene expression profiles of tumor tissue and adjacent non-tumor tissue samples from GC patients were analyzed in R using the “edgeR” package and were normalized by log2 transformation. We used fold change (FC) to characterize the expression differences. Each gene has its associated P values. The “edgeR” package was used to determine the differentially expressed genes with a cutoff of P<0.05 and |logFC |>2 to define the differential expression of genes in GC patients. The unbiased t-test provided by the “Limma” package in “R” was used to evaluate the significant P value of differences in gene expression (12). All the genes were tested by t-test to determine their corresponding P value. A heat map of the top 30 differentially expressed genes were drawn by the “pheatmap” package in “R” (13). The heat map was divided into two categories, the tumor tissue group and the adjacent nontumor tissue group. Red represents the up-regulation of gene expression, and the green represents the down-regulation of gene expression.

Functional enrichment analysis

To better understand the biological functions of the dysregulated genes, GO biological enrichment and KEGG pathways analysis were performed through the “ggplot2” and “clusterProfiler” package in “R” (14). DAVID database was used to carry out functional enrichment analysis (https://david-d.ncifcrf.gov/) (15). GO analysis results included three parts, biological process (BP), cellular component (CC), and molecular function (MF). P<0.05 was considered significant.

Survival analysis to search for the candidate genes

GC samples were divided into two groups according to gene expression: the high expression group and the low expression group. Kaplan-Meier survival analysis was conducted using the “survival” package in “R” to explore the associations between the expression of a specific gene and prognosis of GC patients. We analyzed the top 30 differentially expressed genes from 1,313 genes with expression differences, as well as all the genes enriched in the top 29 KEGG pathways, ranked by P value. The log-rank test was used to determine significant differences in survival curves (16), and P value <0.05 was considered as statistically significant.

Results

Identification of mRNAs in GC

A total of 407 samples, including 375 GC tumor tissue samples and 32 adjacent nontumor tissue samples, were collected for this study. There were 1,313 differentially expressed genes in total, including 781 up-regulated and 532 down-regulated genes identified in GC and matched normal tissues. The cut-off criteria of differentially expressed genes was P<0.05 and |logFC| >2. The volcano plot of the differentially expressed genes is presented in . The red dots represent the up-regulated genes, while the green dots represent the down-regulated genes. The heat map of the top 30 differentially expressed genes ranked according to the fold change was conducted in R with the package “pheatmap” (https://cran.Rproject.org/web/packages/pheatm/apindex.html) and is shown in .

Figure 1

A volcano plot of differentially expression genes in gastric cancer patients. The red dot represents up-regulated genes. The red dot represents up-regulated genes, The green dot represents down-regulated genes.

Figure 2

Heatmap of the top 30 differentially expressed genes ranked according to the fold change, the right side of the sample is the tumor group, the left side of the heatmap is the matched normal tissue group.

Functional analysis

To better understand the genes’ function, Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomics (KEGG) analysis were performed in “R”. The up-regulated and down-regulated genes were separately analyzed by KEGG analysis. The top 29 terms with the lowest P value were selected. The results (, ) showed that the down-regulated genes were significantly enriched in pathways such as the cGMP-PKG, estrogen, and cAMP signaling pathways, as well as the PPAR signaling pathway. The down-regulated genes also interacted with gastric acid secretion, protein digestion and absorption, and insulin secretion. The up-regulated genes were primarily enriched in cytokine-cytokine receptor interaction. The top six GO terms (, ) were “digestion”, “peptide cross-linking”, “keratinocyte differentiation”, “erythrocyte differentiation”, “proteolysis”, and “detection of chemical stimulus involved in sensory perception of the bitter taste”. GO analysis results included BP, CC, and MF, and P<0.05 was considered as statistically significant.

Figure 3

Pathway enrichment map of 1,313 differentially expressed genes, the up-regulated and down-regulated genes were separately analyzed in KEGG analysis. The left side is the down-regulated genes, the right side is the up-regulated genes.

Table 1

Pathway enrichment analysis of the 1,313 differentially expressed genes

Cluster	Description	P value	Count	Gene
Down	Protein digestion and absorption	3.11E−09	16	MEP1B/PGA3/CELA3A/COL4A6/CELA2A/MME/ATP1A2/CPB1/PGA4/SLC15A1/MEP1A/PGA5/XPNPEP2/SLC6A19/SLC7A9/SLC8A2
Down	Metabolism of xenobiotics by cytochrome P450	1.58E−07	13	CYP3A4/AKR1C1/ALDH3A1/GSTA2/ADH7/UGT1A7/GSTA1/GSTM5/ADH4/UGT2B17/UGT1A5/ADH1B/CBR1
Down	Chemical carcinogenesis	3.97E−07	13	CYP3A4/ALDH3A1/GSTA2/ADH7/UGT1A7/GSTA1/GSTM5/ADH4/UGT2B17/UGT1A5/AKR1C2/ADH1B/CBR1
Down	Drug metabolism-cytochrome P450	4.62E−06	11	CYP3A4/ALDH3A1/GSTA2/ADH7/UGT1A7/GSTA1/GSTM5/ADH4/UGT2B17/UGT1A5/ADH1B
Down	Neuroactive ligand-receptor interaction	8.53E−06	25	GRIK5/GIP/CORT/LPAR3/GRIA2/GABRA5/TACR2/HRH3/GHRL/PENK/GRIA4/GALR1/VIP/NPY2R/AVPR1B/ADCYAP1/SST/GRPR/GHR/CCKAR/LEP/GCG/P2RY14/GRP/PTH2R
Down	Retinol metabolism	9.82E−05	9	CYP3A4/ADH7/UGT1A7/RDH12/BCO1/ADH4/UGT2B17/UGT1A5/ADH1B
Down	Fat digestion and absorption	0.000127444	7	MTTP/APOB/LIPF/MOGAT2/FABP1/APOA1/ABCG5
Down	Pancreatic secretion	0.000412655	10	CELA3A/CELA2A/ATP1A2/CPB1/ATP2B3/CLCA1/CLCA4/CCKAR/SLC26A3/KCNMA1
Down	Insulin secretion	0.000665041	9	GIP/CREB3L3/ATP1A2/CACNA1S/ADCYAP1/CCKAR/GCG/KCNMB1/KCNMA1
Down	Gastric acid secretion	0.001163971	8	KCNE2/ATP4B/ATP4A/ATP1A2/SST/SLC9A4/MYLK/GAST
Down	cAMP signaling pathway	0.002601643	14	GIP/CREB3L3/GRIA2/ATP1A2/PLN/GHRL/ATP2B3/GRIA4/VIP/CACNA1S/MYL9/ADCYAP1/SST/GCG
Down	Glycolysis/gluconeogenesis	0.002889032	7	ALDOB/ALDH3A1/ADH7/PCK1/G6PC/ADH4/ADH1B
Down	Starch and sucrose metabolism	0.003138338	5	MGAM/TREH/ENPP3/MGAM2/G6PC
Down	Vitamin digestion and absorption	0.004184031	4	RBP2/APOB/GIF/APOA1
Down	Glutathione metabolism	0.004686566	6	NAT8/GSTA2/GSTA1/GPX3/GSTM5/ANPEP
Down	Estrogen signaling pathway	0.005488775	10	KRT37/ESR1/TFF1/CREB3L3/GNAO1/KRT15/HSPA2/PGR/KRT20/KRT13
Down	Steroid hormone biosynthesis	0.006585195	6	CYP3A4/AKR1C1/UGT1A7/UGT2B17/UGT1A5/AKR1C2
Down	Drug metabolism-other enzymes	0.00666097	7	CYP3A4/GSTA2/UGT1A7/GSTA1/GSTM5/UGT2B17/UGT1A5
Down	cGMP-PKG signaling pathway	0.007128674	11	CREB3L3/ATP1A2/PLN/ATP2B3/CACNA1S/MYL9/SLC8A2/MYLK/KCNMB1/RGS2/KCNMA1
Down	Carbohydrate digestion and absorption	0.007537726	5	MGAM/MGAM2/SLC2A5/ATP1A2/G6PC
Down	Galactose metabolism	0.01062647	4	MGAM/MGAM2/G6PC/AKR1B10
Down	Fructose and mannose metabolism	0.013229258	4	ALDOB/PFKFB1/KHK/AKR1B10
Down	Pentose and glucuronate interconversions	0.014669893	4	UGT1A7/UGT2B17/UGT1A5/AKR1B10
Down	PPAR signaling pathway	0.017632254	6	RXRG/APOC3/PLIN4/PCK1/FABP1/APOA1
Down	Tyrosine metabolism	0.017839499	4	ALDH3A1/ADH7/ADH4/ADH1B
Down	Cell adhesion molecules (CAMs)	0.021986616	9	CD276/CLDN14/NRXN1/CLDN23/MPZ/CLDN17/NEGR1/CNTN1/NRXN3
Down	Glycerolipid metabolism	0.028159835	5	LIPF/MOGAT2/DGKB/AKR1B10/PNLIPRP3
Down	Porphyrin and chlorophyll metabolism	0.029764972	4	UGT1A7/COX10/UGT2B17/UGT1A5
Up	Cytokine-cytokine receptor interaction	0.000755068	22	CXCL8/IL13RA2/IFNL2/IL15RA/CXCL10/INHBA/CXCL9/PF4V1/CXCL6/GDF1/CXCL11/TNFSF11/TNFRSF11B/PPBP/THPO/CXCL5/BMP8A/IL17F/CCL26/GDF15/CCL3/CSH2

Count: the number of enriched genes in each term.

Figure 4

The top 6 GO terms. Count: the number of enriched genes in each term. The blue box represent biological process (BP), the green box represents cellular component (CC), the red box represents molecular function (MF).

Table 2

The 6 GO terms

Category	Term	Name	Count	Percent	P value	FDR
MF	GO:0007586	Digestion	18	1.417322835	0.000000231	0.00042
MF	GO:0018149	Peptide cross-linking	16	1.25984252	0.00000026	0.000472
CC	GO:0030216	Keratinocyte differentiation	18	1.417322835	0.00000399	0.007254384
CC	GO:0030218	Erythrocyte differentiation	18	1.417322835	0.00000399	0.007254384
BP	GO:0006508	Proteolysis	58	4.566929134	0.0000111	0.020147299
BP	GO:0001580	Detection of chemical stimulus involved in sensory perception of bitter taste	12	0.94488189	0.0000259	0.047106424

Count: the number of enriched genes in each term. BP, biological process; CC, cellular component; MF, molecular function.

Survival analysis

To ascertain which candidate genes may influence survival outcomes in GC patients, survival analysis was performed using the Kaplan-Meier method with a log-rank statistical test. The patients with GC were categorized into a high-expression group and a low-expression group according to its median gene expression level. We downloaded the survival information of the samples from the TCGA database and obtained the matched survival status of each sample. We analyzed the top 30 differentially expressed genes using the “survival” package in “R”, ranked according to the fold change, we also analyzed all the genes enriched in the top 29 KEGG pathway, ranked by P value.

Seven genes were found to be associated with survival: ADH4, HRH3, NPY2R, APOA1, GALNT14, LRRIQ1, and CCDC57. A significance level of P<0.05 was set as the cut-off criteria, and the results are shown in . GALNT14, LRRIQ1, and CCDC57 were selected as candidate genes from the top 30 differentially expressed genes, and ADH4, HRH3, NPY2R, and APOA1 were associated with the top 29 KEGG pathways: “Metabolism of xenobiotics by cytochrome P450”, “Chemical carcinogenesis”, “Drug metabolism-cytochrome P450” “Retinol metabolism”, “Glycolysis/Gluconeogenesis” and “Tyrosine metabolism”, HRH3 and NPY2R were found to be linked with the neuroactive ligand-receptor interaction pathway. APOA1 interacted with the fat digestion and absorption, vitamin digestion and absorption, and PPAR signaling KEGG pathways. GC patients who had high expression of LRRIQ1, GALNT14, APOA1, NPY2R, HRH3, and ADH4 had a better prognosis than GC patients with low expressions of these genes, while the GC patients with low expression of CCDC57 often had poor survival outcomes.

Figure 5

The survival analysis of dysregulated genes, red line represents the high expression group of gastric cancer patients while the green line represents the low expression group of patients.

Discussion

GC is one of the most malignant cancers worldwide, although great progress has been made in endoscopic surveillance for early GC, and many new molecular targeted drugs have been invented and clinically applied, such as the human epidermal growth factor receptor 2 (HER-2)-targeted drug trastuzumab (17). In spite of this, the 5-year survival rate remains low (29.6%) for GC patients around the world (1). Many genes are overexpressed in GC, and some of these genes could be potential prognosis predictors and/or therapeutic targets. It has been proved that the accumulation of mutations in crucial genes may cause cancer by altering normal programs of differentiation and cell proliferation and death (18). Genetic changes often lead to the alteration of biological processes. TCGA project wants to identify dysregulated pathways and candidate driver genes in GC (19).

In this study, we conducted some bioinformatic analyses to determine the candidate genes which can predict survival in GC patients. At first, we found a total of 1,313 differentially expressed genes, including 781 up-regulated and 532 down-regulated genes. Among the top 30 differentially expressed genes and all the differentially expressed genes enriched in the top 29 KEGG pathways, 7 genes (ADH4, HRH3, NPY2R, APOA1, GALNT14, LRRIQ1, and CCDC57) were selected as the candidate genes. GC patients with low expression of CCDC57 often had poor survival outcomes. GC patients with low expression of any one of the other six genes (ADH4, HRH3, NPY2R, APOA1, GALNT14, and LRRIQ1) often had a good survival outcome.

acetylgalactosaminyltransferase 14 (GALNT14) belongs to the polypeptide N-acetylgalactosaminyltransferase family. Previous studies found that the loss function of GALNTs can result in altered glycoproteins and can cause tumor aggressiveness in various kinds of cancer (20). The genotype of polypeptide GALNT14 has also be put forward as a potential prognostic predictor for patients undergoing chemotherapy for hepatocellular carcinoma (21). The human alcohol dehydrogenase 4 gene (ADH4) is a member of the human alcohol dehydrogenase (ADH) family, which plays a role in the process of ethanol metabolism (22). Neuropeptide Y (NPY) is an appetite hormone that has been reported to be a candidate gene associated with the development of obesity and control of food intake (23-27). Apolipoprotein AI (APOA1) belongs to the apolipoprotein family (28). By using gene expression array analysis, it has been found that APOA1 mRNA expression in ovarian serous is a marker of longer survival (29). A retrospective study involving 1,201 GC patients who received surgery showed that patients with high ApoB1/ApoA1 (≥1) had shorter overall survival (30). The histamine receptor H3 (HRH3) has been identified as an important molecule in inflammation and carcinogenesis. Recent studies have found that HRH4 is involved in inflammation-related colorectal carcinogenesis (31). Coiled-coil-domain-containing 57 (CCDC57) has been found to be slightly higher in uterine leiomyomata (32). Previous studies have found that the methylation of CCDC57 is related to age, tumor location, and Helicobacter infection in early gastric carcinogenesis (33). Our study shows that these genes are associated with the outcomes of GC patients, but their molecular mechanisms are still poorly understood.

In summary, we selected seven candidate genes that could be considered as candidate prognostic biomarkers in GC patients. These seven genes may become future therapeutic targets in GC. However, our study needs another validation cohort to verify our results, and further investigation and molecular experiments are required to explore the roles of these genes in GC better.