MicroRNA (miR) signatures may aid the diagnosis and prediction of cancer; therefore, miRs associated with the prognosis of esophageal squamous cell carcinoma (ESCC) were screened. miR-sequencing (seq) and mRNA-seq data from early-stage ESCC samples were downloaded from The Cancer Genome Atlas (TCGA) database, and samples from subjects with a >6-month survival time were assessed with Cox regression analysis for prognosis-associated miRs. A further two miR expression datasets of ESCC samples, GSE43732 and GSE13937, were downloaded from the Gene Expression Omnibus database. Common miRs between prognosis-associated miRs, and miRs in the GSE43732 and GSE13937, datasets were used for risk score calculations for each sample, and median risk scores were applied for the stratification of low- and high-risk samples. A prognostic scoring system of signature miRs was subsequently constructed and used for survival analysis for low- and high-risk samples. Differentially-expressed genes (DEGs) corresponding to all miRs were screened and functional annotation was performed. A total of 34 prognostic miRs were screened and a scoring system was created using 10 signature miRs (hsa-miR-140, −33b, −34b, −144, −486, −214, −129-2, −374a and −412). Using this system, low-risk samples were identified to be associated with longer survival compared with high-risk samples in the TCGA and GSE43732 datasets. Age, alcohol and tobacco use, and radiotherapy were prognostic factors for samples with different risk scores and the same clinical features. There were 168 DEGs, and the top 20 risk scores positively-correlated and the top 20 risk scores negatively-correlated DEGs were significantly enriched for six and 10 functional terms, respectively. ‘Tight junction’ and ‘melanogenesis’ were two significantly enriched pathways of DEGs. miR-214, miR-129-2, miR-37a and miR-486 may predict ESCC patient survival, although further studies to validate this hypothesis are required.
With respect to prognosis and mortality, esophageal squamous cell carcinoma (ESCC) is the 8th most common type of cancer and the 6th most common cancer-associated cause of premature mortality (1). Globally, ~450,000 people are affected by ESCC and this incidence is growing (2), with ~500,000 new cases diagnosed each year (3). The 5-year survival for patients with ESCC remains low (15–20%) (4). The cure rate of early stage ESCC is as high as 50% following surgical resection (5), although a number of patients with ESCC are not candidates for surgery due to comorbid conditions, including advanced age. In these cases, the 30-day mortality is 2–10% (6).
Numerous studies have revealed that smoking and pre-diagnosis alcohol consumption are risk factors for ESCC, and the surgical technique, biological behavior, postoperative treatment and response to chemoradiotherapies contribute to improving prognosis (7,8). There are additional genetic alterations that contribute to the prognosis of ESCC, including somatic mutations, copy number variations and gene expression alterations (9). MicroRNAs (miRs) are useful diagnostic and prognostic indicators for human cancer (10), and miR-377 suppresses the initiation and progression of ESCC by inhibiting cluster of differentiation 133 and vascular endothelial growth factor (11). miR-1290 and miR-613 are prognostic factors for patients with ESCC (12,13), and high expression of miR-103/107 is associated with poor survival in patients with ESCC (14). Nevertheless, miRs may cooperate to drive the progression and prognosis of esophageal carcinoma.
miR signatures may aid in the diagnosis and prognosis of cancer (15). Feber et al (16) assessed the association of miR expression with patient survival and lymph node metastasis by evaluating miR expression in 45 primary tumors. This previous study identified that miR profiles have prognostic value for staging patients with ESCC. The present study screened signature miRs involved in predicting ESCC using miR-sequencing (seq) and mRNA datasets from The Cancer Genome Atlas (TCGA; gdc-portal.nci.nih.gov) and the Gene Expression Omnibus (GEO; www.ncbi.nlm.nih.gov) database. Subsequently, a prognostic scoring system was created to identify predictive miRs using sample risk scores. All cancer samples were divided into high- and low-risk categories and validated using the scoring system, and the differentially-expressed genes (DEGs) associated with miRs were functionally annotated.
Materials and methodsMicroarray data
miR-seq and mRNA-seq data from early stage ESCC samples were downloaded from TCGA on March 18, 2017 and 89 samples with miR and mRNA expression data were obtained by matching barcodes. These were early stage (stage I and II) cancer samples. This dataset was used as a test dataset.
A further two miR expression datasets of ESCC samples, GSE43732 and GSE13937, were downloaded from the GEO database. The GSE43732 dataset was based on the platform of GPL16543, and contained 53 early stage cancer samples. The GSE13937 dataset was based on the platform of GPL8835, and contained 31 early stage cancer samples. These two datasets were used as validation datasets. Clinical feature data for all downloaded datasets were also collected (Table I).
Prognostic miRs
The overall prognosis of patients with early stage ESCC is comparatively good. Samples with a <6-month censor time are not representative samples for analyzing prognostic factors. Therefore, miR-seq data samples from TCGA with a survival time of <6 months were removed to avoid introducing more mixed factors, and the remaining 77 samples assessed with Cox regression analysis using the survival package in R (17) to identify prognostic miRs (threshold of P<0.01 for the log rank test).
Prognostic scoring system
Prognostic miRs were matched with miRs in the GSE43732 and GSE13937 datasets, and common ones were collected. Selected miRs were ranked according to log rank P-values to construct a prognosis scoring system. miRs were added singly subsequent to the first three, until the highest P-value representing correlation significance between samples and overall survival time was obtained. When the P-value was greatest, miRs were considered to be signature miRs, and the scoring system was created using these miRs.
Risk scores are used to assess risk factors for large samples (18). Signature miRs were used to calculate risk scores for samples in the TCGA dataset using the following formula:
Risk score = β gene 1 × expr gene 1 + β gene 2 × expr gene 2 + … + β gene n × expr gene n, where β gene indicates the regression coefficients of the gene, and the exp gene indicates its expression levels.
The risk scores of validation samples (GSE43732 and GSE13937) were computed, and a median risk score was applied to stratify low- and high-risk samples. Subsequently, survival correlation coefficients between low- and high-risk samples in the TCGA and GEO datasets, and correlations among risk scores, were assessed. In addition, correlations between clinical features and sample prognosis were analyzed via Cox regression.
Functional annotation of samples with different prognosis risks
The matched RNA-seq data was downloaded from TCGA according to the barcodes of the samples used in the prognostic miRNA analysis. The RNA-seq data was used to screen the DEGs between high- and low-risk samples using the limma package in R (bioconductor.org/packages/release/bioc/html/limma.html) (19). A false discovery rate (FDR) of <0.05 was set as the threshold. Correlation coefficients for gene expression and risk scores were computed, and positively and negatively-correlated genes were annotated with respect to significant functional terms, and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg) pathway terms, using DAVID (david.ncifcrf.gov) (20).
ResultsPrognostic miRs
Using Cox regression analysis on samples that indicated a survival time of >6 months, 34 prognostic miRs from the miR-seq dataset were screened and 16 common miRs were identified between the GSE43732 and GSE13937 datasets (Table II).
Prognostic scoring system
To create a prognostic scoring system, common miRs between prognostic miRs and miRs in the GEO datasets were added singly following the first three, until the highest P-value representing connection significance between samples and overall survival time was obtained.
A prognostic scoring system was created using the 10 signature miRs with the greatest P-values, and low-risk samples had greater survival in the TCGA and GSE43732 datasets. These data appear in Fig. 1A and B. Differences in the GSE13937 dataset were not notable (Fig. 1C). Regression analysis revealed that risk scores were correlated with prognosis (P=0.0141; Table III). Differences in expression among 10 signature genes in samples stratified by clinical features were noted, and Table IV shows the risk factors that were prognostic for samples with different risk scores (P<0.05). Survival curves are presented in Figs. 2–5. Risk scores for samples, survival time and expression clustering heatmaps of the 10 signature miRs from the TCGA, GSE13937 and GSE43732 datasets are in Fig. 6.
Functional annotation of samples with different prognosis risks
In total, 168 DEGs were identified, and 58 were negatively-associated with risk scores, with 110 positively-associated with risk scores. The expression pattern of the top 20 DEGs positively- and negatively-associated with risk scores differed significantly between low and high-risk samples (Fig. 7A). The GO enrichment of the DEGs is presented in Fig. 7B. The top 20 positively-associated DEGs were significantly enriched in six KEGG pathways, including: hsa05217-Basal cell carcinoma, hsa04916-Melanogenesis, hsa04610-Complement and coagulation cascades, hsa04530-Tight junction, hsa04340-Hedgehog signaling pathway and hsa03320-PPAR signaling pathway (Fig. 7C).
Discussion
In order to screen miRs involved in the prognosis of ESCC, miR-seq and mRNA-seq data for early stage ESCC samples were downloaded from TCGA, with a further two miR expression datasets, GSE43732 and GSE13937, downloaded from the GEO database. miR-seq data samples with a survival time of >6 months were subjected to Cox regression analysis to assess prognostic value. Common prognostic miRs, and miRs in the GSE43732 and GSE13937 datasets, were used for risk score calculations, and a median risk score was used to stratify low- and high-risk samples. A prognostic scoring system of 10 signature miRs was made according to survival analysis between low- and high-risk samples. It was noted that low-risk samples had greater survival compared with high-risk samples in the TCGA and GSE43732 datasets. Age, alcohol and tobacco use, and radiotherapy were prognostic factors for samples with different risk scores. The present study identified 168 DEGs for all miRs, 110 of which were positively correlated with risk scores. The top 20 positively-correlated and top 20 negatively-correlated DEGs were significantly enriched in six and 10 functional terms, respectively. There were six significantly enriched KEGG pathways, including ‘tight junction’ and ‘melanogenesis’.
Prognostic scoring is used to predict survival and disease recurrence for a number of types of cancer (21). Wang et al (17) established a 53-gene expression system to be used to predict overall survival for gastric cancer. Mao et al (22) created a 12-gene prognostic scoring system to guide adjuvant therapy for breast cancer. Yang et al (23) created a miR signature to stratify patients with Barrett's esophagus with different prognostic risks for targeted chemoprevention.
A number of miRs in the prognostic system used in the present study have been previously implicated in ESCC or some other malignant tumors. miR-214, a miR that regulates cancer cell proliferation, migration and invasion by targeting phosphatase and tensin homolog in gastric cancer, has been reported to reduce cell survival via downregulation of Bcl2l2 in cervical cancer cells (24,25). The predictive value of miR-214 for prognosis and multidrug resistance has been implicated in ESCC (26). Overexpression has been reported to enhance cisplatin sensitivity in ESCC by directly targeting surviving, and indirectly through CUG triplet repeat RNA binding protein 1 (27). miR-129-2 suppresses the proliferation and migration of ESCC via downregulation of SRY-related HMG box 4, and miR-129 is hypothesized to be a novel therapeutic target and biomarker in gastrointestinal cancer (28,29). miR-37a is a prognostic marker for patient survival in early-stage non-small cell lung cancer (30). miR-39a has been implicated in cell proliferation, migration and invasion in gastric cancer by targeting SRC kinase signaling inhibitor 1 (31). miR-486-5p expression is frequently decreased in human cancer. Low or unaltered expression of miR-486-5p compared with neighboring normal tissues has been demonstrated to be associated with a poor prognosis, and high expression with a good prognosis, in gastric cancer (32). miR-486 was observed to be downregulated in ESCC tissues (33). In patients with ESCC, miR-486-3p was highly expressed following chemotherapy treatment (34). In conclusion, miR-214, miR-129-2, miR-37a and miR-486 may predict survival in patients with ESCC, although these data require validation with larger studies.
Acknowledgements
The present study was funded by the Natural Science Foundation of China Youth Fund Project of Xinjiang (grant no. 2016D01C322).
ReferencesJemalASiegelRXuJWardECancer Statistics, 2010: CA Cancer J Clin6013002010HerszényiLTulassayZEpidemiology of gastrointestinal and liver tumorsEur Rev Med Pharmacol Sci14249258201020496531FerlayJSoerjomataramIDikshitREserSMathersCRebeloMParkinDMFormanDBrayFCancer incidenceand mortality worldwide: Sources, methods and major patternsin GLOBOCAN 2012Int J Cancer136E359E386201510.1002/ijc.2921025220842PennathurAGibsonMKJobeBALuketichJDOesophageal carcinomaLancet381400412201310.1016/S0140-6736(12)60643-623374478RiceTWAdelsteinDJZuccaroGFalkGWGoldblumJRAdvances in the treatment of esophageal carcinomaGastroenterologist527829419979436004BennettCVakilNBergmanJHarrisonROdzeRViethMSandersSGayLPechOLongcroft-WheatonGConsensus statements for management of Barrett's dysplasia and early-stage esophageal adenocarcinoma, based on a Delphi processGastroenterology143336346201210.1053/j.gastro.2012.04.032225376135538857KimTGrobmyerSRSmithRBen-DavidKAngDVogelSBHochwaldSNEsophageal cancer-the five year survivorsJ Surg Oncol103179183201110.1002/jso.2178421259254ChristeinJDHollingerEFMillikanKWPrognostic factors associated with resectable carcinoma of the esophagusAm Surg68258263200211893104WangVEGrandisJRKoAHNew strategies in esophageal carcinoma: Translational insights from signaling pathways and immune checkpointsClin Cancer Res2242834290201610.1158/1078-0432.CCR-16-029227370606FangYXGaoWQRoles of microRNAs during prostatic tumorigenesis and tumor progressionOncogene33135147201410.1038/onc.2013.5423455326LiBXuWWHanLChanKTTsaoSWLeeNPYLawSXuLYLiEMChanKWMicroRNA-377 suppresses initiation and progression of esophageal cancer by inhibiting CD133 and VEGFOncogene3639864000201710.1038/onc.2017.29282881405511242XieRWuSNGaoCCYangXZWangHGZhangJLYanWMaTHPrognostic value of combined and individual expression of microRNA-1290 and its target gene nuclear factor I/X in human esophageal squamous cell carcinomaCancer Biomark20325331201710.3233/CBM-17002928800311GuanSWangCChenXLiuBTanBLiuFWangDHanLWangLHuangXMiR-613: A novel diagnostic and prognostic biomarker for patients with esophageal squamous cell carcinomaTumor Biol3743834391201610.1007/s13277-015-4271-8GuoYChenZZhangLZhouFShiSFengXLiBMengXMaXLuoMDistinctive microRNA profiles relating to patient survival in esophageal squamous cell carcinomaCancer Res682633200810.1158/0008-5472.CAN-06-441818172293HirokiEAkahiraJSuzukiFNagaseSItoKSuzukiTSasanoHYaegashiNChanges in microRNA expression levels correlate with clinicopathological features and prognoses in endometrial serous adenocarcinomasCancer Sci101241249201010.1111/j.1349-7006.2009.01385.x19891660FeberAXiLPennathurAGoodingWEBandlaSWuMLuketichJDGodfreyTELitleVRMicroRNA prognostic signature for nodal metastases and survival in esophageal adenocarcinomaAnn Thorac Surg9115231530201110.1016/j.athoracsur.2011.01.056214200703399250WangPWangYHangBZouXMaoJHA novel gene expression-based prognostic scoring system to predict survival in gastric cancerOncotarget755343553512016274193735342421ZhangCBZhuPYangPCaiJQWangZLLiQBBaoZSZhangWJiangTIdentification of high risk anaplastic gliomas by a diagnostic and prognostic signature derived from mRNA expression profilingOncotarget636643366512015264366994742201RitchieMEPhipsonBWuDHuYLawCWShiWSmythGKlimma powers differential expression analyses for RNA-sequencing and microarray studiesNucleic Acids Res43e47201510.1093/nar/gkv007256057924402510Huang daWShermanBTLempickiRASystematic and integrative analysis of large gene lists using DAVID bioinformatics resourcesNat Protoc44457200910.1038/nprot.2008.21119131956WangLMKevansDMulcahyHO'SullivanJFennellyDHylandJO'DonoghueDSheahanKTumor budding is a strong and reproducible prognostic marker in T3N0 colorectal cancerAm J Surg Pathol33134141200910.1097/PAS.0b013e318184cd5518971777MaoXYLeeMJZhuJZhuCLawSMSnijdersAMGenome-wide screen identifies a novel prognostic signature for breast cancer survivalOncotarget814003140162017281223285355157YangHGuJWangKKZhangWXingJChenZAjaniJAWuXMicroRNA expression signatures in barrett's esophagus and esophageal adenocarcinomaClin Cancer Res1557445752200910.1158/1078-0432.CCR-09-0385197379492745487YangTSYangXHWangXDWangYLZhouBSongZSMiR-214 regulate gastric cancer cell proliferation, migration and invasion by targeting PTENCancer Cell Int1368201310.1186/1475-2867-13-68238349023716801WangFLiuMLiXTangHMiR-214 reduces cell survival and enhances cisplatin-induced cytotoxicity via down-regulation of Bcl2l2 in cervical cancer cellsFEBS Lett587488495201310.1016/j.febslet.2013.01.01623337879ZhouYHongLPrediction value of miR-483 and miR-214 in prognosis and multidrug resistance of esophageal squamous cell carcinomaGenet Test Mol Biomarkers17470474201310.1089/gtmb.2012.0518237213453668496PhatakPByrnesKAMansourDLiuLCaoSLiRRaoJNTurnerDJWangJYDonahueJMOverexpression of miR-214-3p in esophageal squamous cancer cells enhances sensitivity to cisplatin by targeting survivin directly and indirectly through CUG-BP1Oncogene3520872097201610.1038/onc.2015.27126234674KangMLiYLiuWWangRTangAHaoHLiuZOuHmiR-129-2 suppresses proliferation and migration of esophageal carcinoma cells through downregulation of SOX4 expressionInt J Mol Med325158201310.3892/ijmm.2013.138423677061FeslerAZhaiHJuJmiR-129 as a novel therapeutic target and biomarker in gastrointestinal cancerOnco Targets Ther7148114852014251877284149397VõsaUVooderTKoldeRFischerKVälkKTõnissonNRoosipuuRViloJMetspaluAAnniloTIdentification of miR-374a as a prognostic marker for survival in patients with early-stage nonsmall cell lung cancerGenes Chromosomes Cancer50812822201110.1002/gcc.2090221748820XuXWangWSuNZhuXYaoJGaoWHuZSunYmiR-374a promotes cell proliferation, migration and invasion by targeting SRCIN1 in gastric cancerFEBS Lett589407413201510.1016/j.febslet.2014.12.02725554419ChenHRenCHanCWangDChenYFuDExpression and prognostic value of miR-486-5p in patients with gastric adenocarcinomaPLoS One10e0119384201510.1371/journal.pone.0119384257933944368750LiuSGQinXGZhaoBSQiBYaoWJWangTYLiHCWuXNDifferential expression of miRNAs in esophageal cancer tissueOncol Lett516391642201310.3892/ol.2013.1251237618283678876HummelRWangTWatsonDIMichaelMZvan der HoekMHaierJHusseyDJChemotherapy-induced modification of microRNA expression in esophageal cancerOncol Rep2610111017201121743970
Survival curves for patients with early stage esophageal carcinoma stratified by low- and high-risk. Samples from (A) The Cancer Genome Atlas, and (B) GSE43732 and (C) GSE13937 datasets. **P<0.05.
Survival curves of high- and low-risk samples of different ages. (A) Samples <60 years of age. High-risk samples are red and low-risk samples are black. (B) Samples > 60 years of age. High-risk samples are purple and low-risk samples are blue. (C) Combined survival curves of samples with age groups above and below the median age. Curves crossed with P>0.05 represent different samples which cannot be distinguished by risk score, while curves with P<0.05 represent samples that may be distinguished by risk score. **P<0.05.
Survival curves of high- and low-risk samples with different alcohol consumption. (A) Samples from non-drinkers. High-risk samples are red and low-risk samples are black. (B) Samples from drinkers. High-risk samples are purple and low-risk samples are blue. (C) Combied survival curves of samples from drinkers and non-drinkers. Curves crossed with P>0.05 represent different samples which cannot be distinguished by risk score, while curves with P<0.05 represent samples that may be distinguished by risk score. **P<0.05.
Survival curves of high- and low-risk samples from smokers and non-smokers. (A) Samples from non-smokers. High-risk samples are red, and low-risk samples are black. (B) Samples from smokers. High-risk samples are purple, and low-risk samples are blue. (C) Combined survival curves of samples from smokers and non-smokers. (D) Survival curves of non-smokers, smokers and the combination. **P<0.05.
Survival curves of high- and low-risk samples with/without radiation therapy. (A) Samples with no radiation therapy. High-risk samples are red, and low-risk samples are black. (B) Samples with radiation therapy. High-risk samples are purple, and low-risk samples are blue. (C) Combined survival curves from those with/without radiation therapy. **P<0.05.
Risk scores, survival and expression clustering heatmap of the 10 signature microRNAs of all early stage esophageal carcinoma samples. Samples from (A) The Cancer Genome Atlas, and (B) GSE13937 and (C) GSE43732 datasets.
Expression pattern and functional annotation of the DEGs positively- and negatively-associated with risk scores. (A) Expression pattern of the top 20 DEGs positively- and negatively-associated with risk scores. X-axis represents the samples in TCGA dataset, wich risk scores increase from left to right. Y-axis represents the DEGs expression levels. (B) The GO enrichment of the DEGs. (C) KEGG pathway enrichment of the top 20 positively-associated DEGs.
Clinical features of cancer samples downloaded from TCGA and the Gene Expression Omnibus.
Clinical feature
TCGA (n=89)
GSE43732 (n=53)
GSE13937 (n=31)
Age, mean ± standard deviation
63.02±12.44
59.21±9.26
–
Gender
Male
62
43
–
Female
17
10
–
Pathologic_M
M0
79
–
–
M1
0
–
–
Pathologic_N (/N1)
N0
64
43
–
N1
24
10
–
Pathologic_T
T0
1
–
–
T1
24
7
–
T2
33
15
–
T3
31
31
–
Alcohol
Yes
66
32
23
No
22
21
7
Smoking
Yes
17
34
23
No
31
19
7
Reformed
35
–
–
New tumor
Yes
27
–
–
No
60
–
–
Radiation therapy
Yes
18
–
14
No
65
–
17
Mortality
Succumbed
27
24
14
Survived
67
29
17
Overall survival time (months) (mean ± standard)
20.47±20.59
44.7±24.05
29.78±20.89
TCGA, The Cancer Genome Atlas.
Common miRs between prognosis-associated miRs, and miRs in GSE43732 and GSE13937.
miR
P-value
hsa-miR-129-2
3.29 ×10−05
hsa-miR-34b
1.86 ×10−04
hsa-miR-374a
1.92 ×10−04
hsa-miR-412
1.99 ×10−04
hsa-miR-140
4.66 ×10−04
hsa-miR-214
5.15 ×10−04
hsa-miR-144
1.57 ×10−03
hsa-miR-376b
1.59 ×10−03
hsa-miR-486
1.67 ×10−03
hsa-miR-33b
3.99 ×10−03
hsa-let-7f-1
6.22 ×10−03
hsa-miR-494
6.24 ×10−03
hsa-miR-33a
6.37 ×10−03
hsa-miR-432
6.73 ×10−03
hsa-miR-219-1
7.88 ×10−03
hsa-miR-188
9.87 ×10−03
miR, microRNA.
Cox regression results for the prognosis-associated clinical features.
Clinical feature
P-value
Hazards regression (confidence interval)
Age, >60 years vs. <60 years
0.97
1.016 (0.438–2.356)
Sex, male vs. female
0.615
1.325 (0.442–3.97)
Alcohol, yes vs. no
0.916
0.943 (0.318–2.793)
Tobacco, yes vs. no vs. reformed
0.56
0.872 (0.551–1.382)
New tumor, yes vs. no
0.726
1.168 (0.491–2.778)
Radiation therapy, yes vs. no
0.9302
0.951 (0.3113–2.907)
Risk score, high vs. low
0.0141
1.21 (1.005–1.458)
Prognostic factors in high- and low-risk samples under the same clinical features.