Overexpression of YTHDF1 is associated with poor prognosis in patients with hepatocellular carcinoma

2.1 Patients and clinicopathological data

The level-3 expression data (RNA-seqV2) and clinicopathological data of 373 HCC patients and 50 normal samples were downloaded from The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/) data portal. The data on patients’ clinical features including age at diagnosis, days to last follow-up, pathology T stage and pathology N stage were retrospectively abstracted from patient records. All the patients were staged using the 2009 TNM classification of American Joint Committee on Cancer/International Union Against Cancer.

A total of 372 patients with survival data in the TCGA database were given a follow-up exam ranging from 0 to 120.73 months, and the median survival was 19.66 months. For the analysis of survival and follow-up, the date of surgery was used as the beginning of the follow-up period. Overall survival (OS) and disease-free survival (DFS) were calculated in months from the data of diagnosis to the time of death and to the time of tumor progression or recurrence, or death of the patient from HCC, respectively. All patients that died from diseases other than HCC or unexpected events were excluded from the cohort. In this study, we focused on exploring whether m6A related protein were associated with 5-years OS and DFS. In present study, the upper quarter m6A related mRNA expression in all HCC tissues was used as the cutoff point to divide all cases into gene-high and gene-low groups.

2.2 GO and KEGG pathway analysis

To identify functions of m6A related proteins in HCC, we performed GO function enrichment analysis by using MAS 3.0 (http://bioinfo.capitalbio.com/mas3/project/index). KEGG pathway enrichment analysis was also performed to identify pathways enriched in HCC using DAVID system (https://david.ncifcrf.gov/home.jsp) [29]. The q-value 0.05 was considered as significantly.

2.3 Co-expression network construction and analysis

In present study, we performed co-expression analysis of candidate genes by calculating the Pearson correlation coefficient of paired genes. The co-expressed gene pairs with the absolute value of Pearson correlation coefficient ⩾ 0.3 were selected for further study.

OPEN IN VIEWER

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; string-db.org/) database is a pre-computed global resource for the investigation and analysis of associations between proteins. We used STRING to analyze the interactions between DEGs and constructed a PPI network. The database reveals protein interactions, including experimental and predicted protein interaction information. In this study, only experimentally validated interactions with a combined score ⩾ 0.4 were selected as significant. Then, PPI networks were constructed using the Cytoscape software.

OPEN IN VIEWER

Numerical data were presented as the mean ± standard deviation of at least three determinations. Statistical comparisons between groups of normalized data were performed using Student’s t-test using SPSS 13.0 software (SPSS, Inc., Chicago, IL, USA) or Mann-Whitney U-test according to the test condition. Statistical comparisons among multiple groups of normalized data were performed using analysis of v one-way ANOVA according to the test condition. The log-rank test was used to compare differences in survival times. Kaplan Meier and Cox regression analyses were utilized to assess the association between candidate genes and overall survival as well as the prognosis of HCC. 0.05 was considered to indicate a statistically significant difference with a 95% confidence level.

OPEN IN VIEWER

OPEN IN VIEWER

3.1 m6A related proteins were dysregulated in hepatocellular carcinoma

In present study, we analyzed TCGA data to explore the expression patterns of m6A related proteins in hepatocellular carcinoma. A total of 373 HCC patients and 50 normal samples were included in TCGA data. As shown in Fig. 1, we observed seven m6A related genes (YTHDF1, YTHDF2, ALKBH1, ALKBH5, FTO, METTL3, and METTL14) were differentially expressed in HCC compared to match normal patients. Of these genes, YTHDF1 (Fig. 1A, 0.001) and METTL3 (Fig. 1G, 0.001) were up-regulated and YTHDF2 (Fig. 1B, 0.05), ALKBH1 (Fig. 1D, 0.001), ALKBH5 (Fig. 1E, 0.001), FTO (Fig. 1F, 0.05) and METTL14 (Fig. 1H, 0.001) were down-regulated in hepatocellular carcinoma (Fig. 1).

3.2 YTHDF1 was positive correlated with pathology stage in hepatocellular carcinoma

We next analyzed the correlation between the expression levels of m6A related proteins and clinicopathological characteristics. As shown in Fig. 2, we found that YTHDF1 expression was significantly up-regulated in Stage III/IV HCC compared to Stage I (Fig. 2A, 0.001) and Stage II (Fig. 2A, 0.05) HCC. However, we did not observe the difference of expression levels of other m6A related genes among different stages of HCC.

3.3 Overexpression of YTHDF1 is associated with poor prognosis in patients with hepatocellular carcinoma

OS curves were plotted according to candidate gene expression level by the Kaplan-Meier method and log-rank tests. As shown in Fig. 3, overexpression of YTHDF1 (Fig. 3A, p = 0.013) and YTHDF2 (Fig. 3B, p = 0.0235) and down-regulation of FTO (Fig. 3F, p = 0.0038) was correlated with significantly shorter 5-year survival rate. In a multivariate analysis (Table 1), YTHDF1 expression remained an independent prognostic variable for

overall survival. We also explored whether these m6A related genes were associated with disease free survival of HCC patients. Similar to OS analysis, we also found the DFS rates were lower in YTHDF1-high patients in both datasets comparing to YTHDF1-low patients (Fig. 4A, p = 0.0131). Interestingly, we found down-regulation of YTHDF3 (Fig. 4C, p = 0.038), ALKBH5 (Fig. 4E, p = 0.0132), FTO (Fig. 4F, p = 0.0450) and METTL14 (Fig. 4H, p = 0.0306) were correlated with shorter DFS rate (Fig. 4).

OPEN IN VIEWER

Collectively, our results had showed overexpression of YTHDF1 was associated with the progression and could predict poor prognosis of HCC. However, the potential roles of YTHDF1 in HCC remains unclear. To predict the functions of the YTHBF1, we first constructed gene co-expression networks to identify interactions among differentially expressed mRNAs and YTHBF1 by using TCGA. The co-expressed YTHBF1-gene pairs with the absolute value of Pearson correlation coefficient ⩾ 0.3 were selected and the co-expression network was established. Based on the information in the STRING database, we constructed a protein-protein interaction network of YTHBF1 co-expressed genes in the HCC (Fig. 5A). We could observe three PPI hubs were collected using these genes.

MAS 3.0 system was used to perform GO and pathway analysis on the co-expressed genes. A total of 413 differentially expressed genes were classified according to GO term. Our analysis revealed that these co-expressed genes were associated with transcription, cell cycle, cell division, mitosis, oxidation reduction, response to DNA damage stimulus, DNA repair, DNA replication, RNA splicing, chromatin modification, microtubule-based movement, protein amino acid phosphorylation, modification-dependent protein catabolism, proteolysis, lipid metabolism (Fig. 5B). KEGG pathway enrichment analyses showed co-expre- ssed genes of YTHBF1 mainly participate in regulating Valine, leucine and isoleucine degradation, Cell cycle, Complement and coagulation cascades, Propanoate metabolism, Fatty acid metabolism, Butanoate metabolism, Lysine degradation, beta-Alanine metabolism, Tryptophan metabolism, Caprolactam de- gradation, DNA polymerase, PPAR signaling pathway, Homologous recombination, Limonene and pinene degradation, Fatty acid elongation in mitochondria, etc. (Fig. 5C).