Identification of key pathways regulated by a set of competitive long non-coding RNAs in oral squamous cell carcinoma

Abstract

Objective

The aim of this study was to identify important pathways regulated by a set of long non-coding RNAs (lncRNAs) in oral squamous cell carcinoma (OSCC).

Methods

A lncRNA-mediated competitive endogenous RNA network (LMCN) was constructed using information on microRNA (miRNA)–mRNA interactions and lncRNA–miRNA intersections from the E-GEOD-37991 transcription profiling data in the ArrayExpress database. A random walk with restart ranking algorithm was then applied to evaluate the influences of protein-coding genes regulated by competitive lncRNAs. Pathway enrichment scores were calculated based on the propagation scores of protein-coding genes. Finally, permutation tests were used to estimate the significance of the pathways.

Results

We obtained lncRNA–mRNA interactions based on miRNAs common to both miRNA–mRNA interactions and lncRNA–miRNA intersections, and used interactions with a z-score > 0.7 to construct a LMCN. Ten lncRNAs were identified as source nodes in the LMCN, and nine pathways with enrichment scores >0.8, including ‘Cell cycle’, ‘Endocytosis’, and ‘Pathways in cancer’, were significantly enriched by these source nodes.

Conclusions

Nine significant pathways regulated by a set of competitive lncRNAs were identified in OSCC, which may play important roles in the development of OSCC via the cell cycle and endocytosis.

Introduction

Oral cancer may occur in the lips, floor of the mouth, upper and lower alveolar ridges, buccal mucosa, hard palate, sublingual region, retromolar trigone, and anterior two-thirds of the tongue. ¹ Oral cancer is the third most common carcinoma in developing countries and the sixth most common carcinoma worldwide, ² with oral squamous cell carcinoma (OSCC) accounting for more than 95% of all malignancies in the oral cavity. ³ Despite significant progress in the diagnosis and clinical treatment of OSCC in recent decades, the overall survival rate has not improved significantly. ⁴

Long non-coding RNAs (lncRNAs) are a class of non-protein coding RNAs of more than 200 nucleotides in length, which have been discovered in a wide range of biological processes to date. ⁵ Many of these lncRNAs have been shown to be uniquely expressed in specific cancer types or differentiated tissues, and increasing evidence suggests that lncRNAs perform key roles in promoting cancer development and progression. ⁶ LncRNA-UCA1 has been proposed to play an oncogenetic role in OSCC through regulating the Wnt/b-catenin signaling pathway, ⁷ while lncRNA-MALAT1 was found to promote tumor metastasis and growth through inducing epithelial-mesenchymal transition in OSCC. ⁸

Although dysregulation of lncRNAs has been shown in OSCC, lncRNA profiles in this type of cancer remain largely unknown. In this study, we applied a novel method, LncRNAs2Pathways, ⁹ to identify important lncRNA-related pathways in OSCC. LncRNAs2Pathways identifies pathways influenced by the set of lncRNAs of interest based on a global network propagation algorithm. ⁹ This method has been successfully applied in glioma, and in pancreatic and prostate cancers, with most of the identified pathways proposed to be related to the progression of prostate cancer. ⁹

The aim of this study was to characterize key pathways regulated by a set of competitive lncRNAs in OSCC. We collected microRNA (miRNA)–mRNA interactions and lncRNA–miRNA intersections from the starBase v2.0 database, and transcription profiling data from the ArrayExpress database and processed them using the LncRNAs2Pathways method to identify significant pathways likely to be regulated by a set of competitive lncRNAs.

Materials and methods

Data

Transcription profiling data in E-GEOD-37991 ¹⁰ , which includes transcriptional filing data for tumor and paired non-tumor samples, were downloaded from the ArrayExpress database. Control (n = 40) and OSCC data (n = 40) were collected. According to the probe ID and gene symbol in the annotation file, the probe ID was changed to its gene symbol using the feature function Filter. ¹¹ A total of 18,631 genes were finally obtained.

A total of 10,212 pairs of lncRNA–miRNA intersections and 423,975 pairs of miRNA–mRNA interactions were downloaded from the starBase v2.0 database.

The 18,631 genes with transcription data were mapped to the lncRNA–miRNA intersections and miRNA–mRNA interactions, and mRNAs (n = 9125) and lncRNAs (n = 10) common to both were collected as new profiling data (n = 9135). Interactions containing genes from the new profiling data were then collected. A total of 221 pairs of lncRNA–miRNA intersections and 293,185 pairs of miRNA–mRNA interactions were finally obtained.

Construction of competitive endogenous RNA (ceRNA) network

We evaluated the significance of the miRNAs common to both mRNA–miRNA interactions and lncRNA–miRNA interactions using a hypergeometric test, as follows: ¹²

$$P=1-{\displaystyle \sum _{t=0}^x}\frac{\left(\begin{array}{c}K\\ t\end{array}\right)\left(\begin{array}{l}N-K\hfill \\ M-t\hfill \end{array}\right)}{\left(\begin{array}{c}N\\ M\end{array}\right)}$$

(1)

The P-values were subjected to false discovery rate correction. An adjusted P < 0.05 was used as the threshold. Finally, 10 lncRNAs, 6206 mRNAs, and 10,772 interactions were obtained.

We characterized the lncRNA–mRNA interactions by calculating Pearson’s correlation coeffcients according to the expression of the competing lncRNA–mRNAs interactions. Correlation coefficients were then converted to the z-statistic according to Fisher’s z transform method, ¹³ as follows:

$$\text{z}=0.\text{5}\left[\text{ln}\left(\text{1}+\text{r}\right)-\text{ln}\left(\text{1}-\text{r}\right)\right]$$

(2)

where r is Pearson’s correlation coefficient. LncRNA–mRNA interactions with z > 0.7 were considered as co-expression networks and used to construct a lncRNA-mediated ceRNA (LMCN) network, including 10 lncRNAs and 481 mRNAs. The lncRNAs in the LMCN network were considered as source nodes.

Evaluation of influence of protein-coding genes regulated by competitive lncRNAs

Protein-coding genes are more likely to be regulated if they are located close to the source nodes. We evaluated the influence of protein-coding genes regulated by competitive lncRNAs using the random walk with restart ranking algorithm, ¹⁴ as follows:

$${\text{p}}^{\text{t}+\text{1}}=\left(\text{1}-\text{r}\right){\text{Mp}}^{\text{t}}+{\text{rp}}^0$$

(3)

where r indicates the restart probability of the walk in every time step at the source node. In this study, r was set at 0.7.

The probabilities of protein-coding genes were then normalized and defined as propagation scores. All the protein-coding genes were ranked in the gene list L, based on their propagation scores.

Calculating pathway enrichment scores (ES) to evaluate biological pathways

Pathway ES were calculated based on propagation scores of protein-coding genes. Human-related pathways were downloaded from the Kyoto Encyclopedia of Genes and Genomes database. Pathways with > 500 or < 15 protein-coding genes were deleted.

Protein-coding genes in pathways were mapped to the gene list L. The mapped pathway was more likely to be influenced by the combinatorial effects of the competitive lncRNAs if the genes in that pathway occurred near the top of the list L. We identified ES according to the method proposed by Han et al., ⁹ as follows:

$$F_{\mathrm{In}P}\left(i\right)=\sum _{\begin{array}{c}{g}_j\in P\\ j\le i\end{array}}\frac{{\left|tj\right|}^p}{N_R}$$

(4)

$$N_R=\sum _{g_j\in P}{\left|tj\right|}^p$$

(5)

$$F_{\mathrm{Not}P}\left(i\right)=\mathrm{ }\sum _{\begin{array}{c}{g}_j\notin P\\ j\le i\end{array}}\frac{1}{N_{\mathit{NotP}}}$$

(6)

$$\mathrm{ES}\left(\mathrm{P}\right)=\underset{i\in L}{\max }\left\{F_{\mathrm{In}P}\left(i\right)-F_{\mathrm{Not}P}\left(i\right)\right\}$$

(7)

where P was set at 1.

Pathway significance analysis

We estimated the significance of the pathway ES using permutation tests. ⁹ After 1000 permutation tests, the P-value of the observed ES was calculated by comparison with the set of scores in ES_NULL, shown as P-value = M/N, where N indicates the number of permutation tests and M indicates the number of ES_NULL values that were greater than the observed ES. The P-value was then adjusted by the false discovery rate, ¹⁵ and pathways with adjusted P < 0.05 were selected as target pathways regulated by the set of competitive lncRNAs.

No animal or human tissue or cells were directly involved in this study. The transcription profiling data were obtained from the ArrayExpress database and ethical permission was therefore unnecessary.

Results

Construction of LMCN network

We identified lncRNA–mRNA interactions based on common miRNAs between lncRNA–miRNA intersections and miRNA–mRNA interactions, and used interactions with a z-score > 0.7 to construct the LMCN network. The network included 10 lncRNAs and 481 mRNAs (Figure 1).

Figure 1.

Long non-coding RNA (lncRNA)-mediated competitive endogenous RNA network. The network included 10 lncRNAs and 481 mRNAs. Competitive lncRNAs were considered as source nodes, including MALAT1, SNHG10, SNHG5, SNHG8, HCP5, SNHG3, H19, MEG3, DGCR5, and CASC2

Evaluation of influence of protein-coding genes regulated by competitive lncRNAs

We evaluated the influence (propagation scores) of protein-coding genes regulated by competitive lncRNAs and ranked the protein-coding genes in the gene list L, according to their propagation scores. We considered that the greater the probability of a protein-coding gene being located closer to the source nodes, the more this gene was likely to be influenced by the source nodes.

Identification of significant pathways

Based on pathway ES, we performed 1000 permutation tests to estimate the significance of the pathways. Nine significant pathways with ES values >0.8 were finally identified, including ‘Cell cycle’, ‘Endocytosis’, ‘Pathways in cancer’, ‘Chemokine signaling pathway’, ‘Axon guidance’, ‘Spliceosome’, ‘Focal adhesion’, ‘MAPK signaling pathway’, and ‘Regulation of actin cytoskeleton’ (Table 1).

Table 1.

Pathways regulated by the set of competitive long non-coding RNAs

Pathway	Enrichment score	P-value	Adj. P-value
Cell cycle	0.9	<0.01	<0.01
Endocytosis	0.89	<0.01	<0.01
Pathways in cancer	0.89	<0.01	<0.01
Chemokine signaling pathway	0.85	<0.01	<0.01
Axon guidance	0.85	<0.01	<0.01
Spliceosome	0.84	<0.01	<0.01
Focal adhesion	0.83	<0.01	<0.01
MAPK signaling pathway	0.83	<0.01	<0.01
Regulation of actin cytoskeleton	0.81	<0.01	<0.01

Adj. P-value: P-value adjusted by false discovery rate

Discussion

In this study, we applied the LncRNAs2Pathways method to identify important pathways regulated by a set of lncRNAs of interest. Nine pathways were identified, of which the top three pathways were ‘Cell cycle pathway’, ‘Endocytosis pathway’, and ‘Pathways in cancer’.

‘Cell cycle pathway’ plays a vital role in cell differentiation, angiogenesis, and cell proliferation and survival. Cell cycle progression of mitotic cells is accomplished via a sequence of events, including DNA replication and mitosis, with temporal gaps known as G1 and G2 phases. Several genes in this pathway have been reported to be closely related to cancers, including the genes encoding cyclin-dependent kinases, p53, and retinoblastoma protein (Rb). Cyclin-dependent kinases are critical regulatory enzymes affecting all cell cycle transitions. ¹⁶ p53 is a tumor suppressor protein that has been widely reported in many kinds of human cancers,^17,18 and p53 mutation is generally accepted to be a common event in human cancers, including OSCC. ¹⁹ Rb is also a well-known tumor-suppressor, which functions through transition of the restriction point in progression of the cell cycle.

‘Endocytosis’ is responsible for bringing ligands, nutrients, plasma membrane proteins, and lipids into the cell interior and removing them from the cell surface. Oncogenic mutants of growth factor receptors show impaired endocytosis, and this peculiar endocytic system in cancer could be exploited by pharmacological therapies, to deliver effective patient treatment. ²⁰ Although endocytosis was characterized as an important pathway in OSCC, more studies are needed to clarify the mechanisms of this pathway in OSCC.

‘Pathways in cancer’ involves several signaling pathways and tumor-related genes, including Wnt, Ras, p53, and cell apoptosis. Aberrant activation of the Wnt pathway is highly associated with tumorigenesis and tumor metastasis. ²¹ Ras is a small G protein that cycles between inactive GDP-bound and active GTP-bound forms, and functions as a molecular switch for signal transduction initiated in the cell membrane. Mutationally activated Ras proteins have been shown to induce nearly 30% of all human cancers. ²² Cell apoptosis plays an important role in cancer, and de-regulation of apoptotic signaling and activation of anti-apoptotic systems allow cancer cells to evade apoptosis, leading to uncontrolled cell proliferation, tumor progression, therapeutic resistance, and cancer recurrence. ²³

Although the functions of some of the nine significant identified pathways have been investigated in OSCC, most have not. Further studies are therefore needed to clarify the roles of these other pathways, to provide new insights into the development and progression of OSCC.

This study had some limitations. Some of the lncRNAs identified by this method as being associated with OSCC remain to be validated. Furthermore, this method regarded all lncRNAs as equally important, and the extent of the association between individual lncRNAs and cancer were not considered.

In conclusion, we identified nine significant pathways that were regulated by a set of competitive lncRNAs in OSCC, suggesting that these lncRNAs may play important roles in the development of OSCC via the cell cycle and endocytosis.

Declaration of conflicting interest

The authors declare that there is no conflict of interest.