Integrated analysis of circular RNA-associated ceRNA network in pancreatic ductal adenocarcinoma.

Circular RNAs (circRNAs) have displayed dysregulated expression in several types of cancer. However, the functions of the majority of circRNAs in pancreatic ductal adenocarcinoma (PDAC) remain unknown. The present study aimed to investigate the expression, functions and molecular mechanisms of circRNAs in PDAC. The circRNAs, mRNAs and the microRNA (miRNAs) expression profiles were obtained from three Gene Expression Omnibus microarray datasets, and a circRNA-miRNA-mRNA and circRNA-miRNA-hubgene network was established. The interactions between proteins were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins database, and hubgenes were identified using the MCODE plugin. A total of eight differentially expressed circRNAs (DEcircRNAs), 44 differentially expressed miRNAs (DEmiRNAs), and 2,052 differentially expressed mRNAs (DEmRNAs) were identified. The present study successfully constructed a circRNA-miRNA-mRNA competing endogenous RNA (ceRNA) network based on four circRNAs, six miRNAs and 111 mRNAs in PDAC. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathways analyses indicated that DEmRNAs may participate in the pathogenesis and progression of PDAC. The protein-protein interaction network and module analysis identified six hubgenes (THBS1, FN1, TIMP3, TGFB2, ITGA1 and ITGA3). Furthermore, the circRNA-miRNA-hubgene regulatory modules were constructed based on the three DEcircRNAs, one DEmiRNAs and five DEmRNAs. In conclusion, the results of the present study improve the current understanding of the pathogenesis of PDAC. Copyright: © Song et al.

Introduction

As the incidence and mortality rates of pancreatic ductal adenocarcinoma (PDAC) continue to increase annually, it has been estimated to become the second leading cause of cancer-associated mortality in Europe and the USA by 2030 (1,2). To date, surgery is the only treatment option available; however, the 5-year overall survival (OS) time remains unsatisfactory (1,3). This is largely due to a low early diagnostic rate and the fact that the majority of patients exhibit local invasion or distant metastasis (4). In addition, systemic chemotherapy has a limited impact and significant toxicity on the treatment of patients with advanced-stage PDAC. In the majority of cases, these patients are largely resistant to molecularly-targeted agents and immunotherapy (5). Therefore, it is essential to understand the potential mechanism of the carcinogenesis of PDAC and to identify novel markers for developing more effective therapeutic approaches.

Circular RNAs (circRNAs) are a class of noncoding RNAs with continuous and covalently closed circular structures (6). These molecules are not easily degraded by nucleases in the absence of free 3′ and 5′ends, which makes them more stable than the majority of linear RNAs (7,8). The continuous development of high-throughput sequencing technologies and analysis has allowed for the identification of numerous circRNAs that are abnormally expressed in tumor tissues and have an important influence on tumor progression (9–11). CircRNAs can function as a microRNA (miRNA) sponge to repress miRNA function using miRNA response elements (MREs). This inhibits the activity of miRNAs and regulates the expression of their downstream target genes in numerous types of malignancies (12). Rao et al (13) demonstrated that circRNA-0007874 (circMTO1) expression was decreased in glioblastoma tissues; moreover, elevated circMTO1 expression is known to inhibit cell proliferation and promote apoptosis in both in vivo and in vitro conditions. miR-630 is a targeted miRNA of circMTO1. Therefore, Rao et al (13) established a circMTO1/miR-630/temozolomide (TMZ) competing endogenous RNA (ceRNA) network, suggesting that circMTO1 could reverse chemical resistance to TMZ by regulating miR-630. Furthermore, circZFR interacts with C8orf4 through sponging miR-1261 in papillary thyroid carcinoma (14).

In the present study, differentially expressed gene (DEG) expression profiles were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The flowchart for this procedure is presented in Fig. 1. After predicting the sponge miRNA of circRNA and miRNA-mRNA pairs, the present study successfully established the circRNA-miRNA-mRNA network. Subsequently, the present study performed a series of analyses, including functional enrichment analyses and the interactions between proteins. The circRNA-miRNA-hubgene regulatory modules were constructed in order to better understand the pathogenesis of PDAC.

Figure 1.

Flowchart of ceRNA network analysis. circRNA, circular RNA; MREs, miRNA response elements; miRNAs, micro RNAs; ceRNA, competing endogenous RNA; PPI network, protein-protein interaction network; GO, gene ontology; KEGG, kyto encyclopedia of genes and genomes.

Materials and methods

Microarray data processing

The present study downloaded two circRNA expression profiles [GSE69362 (15) and GSE79634 (16)] from the GEO database based on the GPL19978 platform (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi). The present study also downloaded the expression profiles of mRNA and miRNA [GSE60980 (17)] from GEO database based on GPL14550 and GPL15159 platforms (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi). The GSE69362 dataset included 31 normal pancreatic tissues and six PDAC tissues. The datasets of GSE79634 included 20 PDAC tissues and 20 paracancerous tissues. The array data for GSE60980 included the miRNA expression profiles of 51 PDAC tissues and six normal tissues, and mRNA expression profiles of 49 PDAC tissues and 12 normal tissues. No ethical approval nor informed consent was required in the present study due to the data being publicly available from the GEO.

Screening for DEGs

The raw data from the microarray datasets were normalized and log2-transformed. The DEGs of each dataset were identified using the Limma package (version 3.40.6) in the Bioconductor package (18). Subsequently, the present study integrated and ranked the differentially expressed circRNAs (DEcircRNAs) with a robust rank aggregation method (19). The FDR <0.05 and |log2 fold change (FC)| >1 were considered as the threshold values for DEGs selection.

Construction of the ceRNA network

The Circular RNA Interactome (https://circinteractome.nia.nih.gov/) and Cancer-Specific CircRNA databases (http://gb.whu.edu.cn/CSCD/) were used to predict the regulatory relationships between circRNAs and miRNAs. Only the overlapping genes were selected as candidate target miRNAs, which were further screened by the differentially expressed miRNAs (DEmiRNAs). Subsequently, the present study used miRTarBase (20) and TargetScan (21) databases to identify miRNA targeted mRNAs. Only the mRNAs recognized by both databases were considered as candidate mRNAs and were subsequently intersected with the differentially expressed mRNAs (DEmRNAs) in order to determine the DEmRNAs that were targeted by the DEmiRNAs. The circRNA-miRNA-mRNA network was established using a combination based on circRNA-miRNA pairs and miRNA-mRNA pairs and was visualized using Cytoscape software (version 3.7.0; http://cytoscape.org/).

Functional enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted using the Cluster Profiler package (version 3.12.0) of R software (version 3.6.1; http://www.r-project.org) (22), in order to assess the primary function of the DEmRNAs in the ceRNA network in tumorigenesis.

Construction of the protein-protein interaction (PPI) network

The present study established a PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) in order to assess the interactions between DEmRNAs. Cytoscape 3.7.0 was used for visualization. The MCODE plugin was then used to extract hub genes from the PPI network (23).

Statistical analysis

All data were analyzed using R version 3.6.1. The paired Student's t-test was performed to compare the DEGs between the PDAC tissues and paracancerous tissues and FDR filtering was used for comparative analysis. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of differentially expressed circRNA, mRNA and miRNA

In total, 282 and 174 DEcircRNAs were identified from the GSE79634 and GSE69362 datasets, respectively. Among these, 120 upregulated and 162 downregulated circRNAs in the GSE79634 dataset (Fig. 2A), and 116 upregulated and 58 downregulated circRNAs in the GSE69362 dataset (Fig. 2B), were identified. All DEcircRNAs in the GSE79634 and GSE69362 datasets are presented in Fig. 3A. The GSE60980 dataset included mRNA and miRNA expression profiles. Based on this dataset, a total of 44 DEmiRNAs (17 upregulated and 27 downregulated miRNAs) and 2,052 DEmRNAs (1,036 upregulated and 1,016 downregulated mRNAs) were identified in PDAC (Fig. 2C and D). The basic information of the three datasets is listed in Table I. The DEcircRNAs of GSE79634 and GSE69362 datasets were integrated and ranked using a robust method (Table II). In total, eight DEcircRNA (six upregulated and two downregulated circRNAs) were identified (P<0.05; Fig. 3B). The basic characteristics of the eight circRNAs are presented in Table III. The basic structural patterns of the six circRNAs are presented in Fig. 4.

Figure 2.

Volcano plot of DEGs of the three microarray datasets. (A) GSE79634 (circRNAs). (B) GSE69362 (circRNAs). (C) GSE60980 (miRNAs). (D) GSE60980 (mRNAs). Red indicates upregulated DEGs and green indicates downregulated DEGs. DEGs, differentially expressed genes.

Figure 3.

Heatmap of the DEcircRNAs on the GSE79634 and GSE69362 microarray datasets. (A) All DEcircRNAs. (B) Eight DEcircRNAs with robust rank aggregation method. The node color changes gradually from green to red in ascending order according to the log2(foldchange) of genes. DEcircRNAs, differentially expressed circRNAs.

Table I.

Basic information of the three microarray datasets from Gene Expression Omnibus.

Author	Year	Data source	Platform	Geographical location	Sample size (T/N)	No. of genes	(Refs.)
Li et al	2015	GSE69362	GPL19978	China	6/6	4094 circRNAs	(15)
Guo et al	2016	GSE79634	GPL19978	China	20/20	1836 circRNAs	(16)
Sandhu et al	2015	GSE60980	GPL15159	Norway	51/6	1368 miRNAs	(17)
Sandhu et al	2015	GSE60980	GPL14550	Norway	49/12	42545 mRNAs	(17)

circRNAs, circular RNAs; miRNAs, micro RNAs.

Table II.

Total of 8 differentially expressed circRNAs using robust rank aggregation method.

circRNA ID	logFC	P-value	FDR
hsa_circ_0013912	2.751753	1.58×10−11	7.76×10−10
hsa_circ_0092314	2.01195	1.66×10−10	4.63×10−09
hsa_circ_0006220	3.736558	0.000103	0.00307
hsa_circ_0043278	3.426752	0.000141	0.00405
hsa_circ_0000977	3.517108	0.000201	0.00555
hsa_circ_0001666	3.124223	0.013511	0.023253
hsa_circ_0013587	−2.15798	1.51×10−09	2.44×10−08
hsa_circ_0092367	−1.81243	4.40×10−07	2.72×10−06

circRNA, circular RNA; FC, fold change; FDR, false discovery rate.

Table III.

Basic characteristics of the eight differently expressed circRNAs.

circRNA ID	Position	Genomic length	Strand	Best transcript	Gene symbol	Regulation
hsa_circ_0013912	chr1:145601529-145601852	323	Antisense	NM_006468	POLR3C	Up
hsa_circ_0092314	chr22:20113099-20113439	340	Sense	NM_002882	RANBP1	Up
hsa_circ_0006220	chr17:35800605-35800763	158	Sense	NM_001488	TADA2A	Up
hsa_circ_0043278	chr17:35797838-35800763	2925	Sense	NM_001488	TADA2A	Up
hsa_circ_0000977	chr2:10784445-10808849	24404	Antisense	NM_024894	NOL10	Up
hsa_circ_0001666	chr6:170626457-170639638	13181	Sense	NM_032448	FAM120B	Up
hsa_circ_0013587	chr1:113661854-113662145	291	Sense	NM_014813	LRIG2	Down
hsa_circ_0092367	chr15:25325262-25326442	1180	Sense	NR_003329	SNORD116-14	Down

circRNA, circular RNA.

Figure 4.

Structural patterns of the seven circRNAs. (A) hsa_circ_0013912, (B) hsa_circ_0006220, (C) hsa_circ_0043278, (D) hsa_circ_0000977, (E) hsa_circ_0001666, (F) hsa_circ_0013587. Different colors in the circular structure of the circRNA represents the position of the exon.

The potential miRNAs targets of the eight DEcircRNAs were retrieved from the CSCD and CircInteractome online database. A total of 409 circRNA-miRNA pairs were identified. After intersecting with the DEmiRNAs, only 10 circRNA-miRNA pairs, including four circRNAs (hsa_circ_0006220, hsa_circ_0043278, hsa_circ_0001666 and hsa_circ_0092367) and six DEmiRNAs (hsa-mir-1, hsa-mir-214, hsa-mir-224, hsa-mir-223, hsa-mir-1305 and hsa-mir-375), remained. Subsequently, the miRTarBase and TargetScan databases were used to identify target mRNAs of six DEmiRNAs. After the targeted mRNAs intersected with DEmRNAs, the remaining DEmRNAs were used as candidate genes. The results indicate that the ceRNA network included 111 DEmRNAs. Finally, the present study constructed a ceRNA network based on the four circRNAs, six miRNAs and 111 mRNAs (Fig. 5).

Figure 5.

ceRNA network of circRNA-miRNA-mRNA in pancreatic ductal adenocarcinoma. The network consisting of four circRNA nodes, six miRNA nodes, 111 mRNA nodes. Triangles indicate circRNAs, rounded rectangles indicate miRNA, and octagons indicate mRNA. The nodes highlighted in red represent upregulation and the nodes in blue represent downregulation. circRNA, circular RNA; miRNA, micro RNA.

Functional enrichment analysis of DEmRNA

GO and KEGG pathways analyses were performed in order to investigate the biological function of the 111 DEmRNAs. Among the 129 biological process terms, the most enriched GO terms were ‘extracellular matrix organization’ and ‘extracellular structure organization’ (P<0.05). The mRNAs associated with cellular components were most relevant to the extracellular matrix (P<0.05). In terms of molecular function, DEmRNAs were primarily enriched in ‘cell adhesion molecule binding’ (P<0.05). The GO terms are presented in Table IV. Furthermore, the KEGG pathway analysis indicated that the majority of the DEmRNAs were involved in ‘focal adhesion’ and ‘microRNAs in cancer’. The KEGG pathways are presented in Table V.

Table IV.

GO terms enriched by DEmRNA involved in the ceRNA network.

A, Biological process
Terms	Functional description	P-value	Gene count
GO:0030198	‘Extracellular matrix organization’	5.58×10−08	17
GO:0043062	‘Extracellular structure organization’	2.77×10−07	17
GO:0048732	‘Gland development’	3.19×10−05	15
GO:0001655	‘Urogenital system development’	3.19×10−05	13
GO:2000826	‘Regulation of heart morphogenesis’	3.19×10−05	6
B, Cellular component
Terms	Functional description	P-value	Gene count
GO:0031012	‘Extracellular matrix’	3.73×10−06	16
GO:0005925	‘Focal adhesion’	3.24×10−05	13
GO:0005924	‘Cell-substrate adherens junction’	3.24×10−05	13
GO:0030055	‘Cell-substrate junction’	3.24×10−05	13
GO:0044420	‘Extracellular matrix component’	0.000276	7
C, Molecular function
Terms	Functional description	P-value	Gene count
GO:0001968	‘Fibronectin binding’	0.000146	5
GO:0050839	‘Cell adhesion molecule binding’	0.000275	14
GO:0005201	‘Extracellular matrix structural constituent’	0.000802	6
GO:0019838	‘Growth factor binding’	0.001385	7
GO:0005518	‘Collagen binding’	0.002286	5

GO, Gene Ontology.

Table V.

Kyoto Encyclopedia of Genes and Genomes pathways enriched by differentially expressed mRNA involved in the competing endogenous RNA network.

Pathway ID	Functional description	P-value	Genes	Count
hsa04510	Focal adhesion	2.25×10−05	BIRC3/COL4A1/FN1/ITGA1/ITGA3/LAMA4/MET/THBS1/VASP	9
hsa04512	ECM-receptor interaction	4.13×10−05	COL4A1/FN1/ITGA1/ITGA3/LAMA4/THBS1	6
hsa05222	Small cell lung cancer	8.41×10−05	BIRC3/COL4A1/E2F1/FN1/ITGA3/LAMA4	6
hsa05206	MicroRNAs in cancer	0.0001	E2F1/FSCN1/MET/NOTCH3/PDCD4/PIM1/SERPINB5/TGFB2/THBS1/TIMP3	10
hsa05410	Hypertrophic cardiomyopathy (HCM)	0.00047	ITGA1/ITGA3/TGFB2/TPM2/TPM4	5
hsa05414	Dilated cardiomyopathy (DCM)	0.000683	ITGA1/ITGA3/TGFB2/TPM2/TPM4	5
hsa05205	Proteoglycans in cancer	0.000961	FN1/MET/PDCD4/TGFB2/THBS1/TIMP3/TWIST1	7
hsa04933	AGE-RAGE signaling pathway in diabetic complications	0.001052	COL4A1/F3/FN1/PIM1/TGFB2	5

Construction of PPI network and module analysis

The PPI network was constructed, and included 37 nodes and 36 edges (Fig. 6A), following the removal of unconnected nodes. The degree, betweenness centrality and key circRNA-miRNA-mRNA regulatory modules were extracted using the MCODE approach (24) from the PPI network in order to investigate the hubgenes in the process of PDAC carcinogenesis. The significant module contained six nodes and 11 edges. The hubgenes included THBS1, FN1, TGFB2, ITGA1, ITGA3 and TIMP3 (Fig. 6B). Subsequently, the present study established a circRNA-miRNA-hubgene subnetwork (Fig. 7), including 16 subnetwork regulatory modules. Since the expression levels of hsa_circ_0092367 and TGFB2 were inconsistent, the hsa_circ_0092367/hsa-mir-375/TGFB2 regulatory module was excluded.

Figure 6.

Identification of hubgenes from the PPI network with the MCODE algorithm. This network consists of 37 nodes and 36 edges. (A) The PPI network of 111 genes. (B) The PPI network of six hubgenes that were extracted from (A).

Figure 7.

circRNA-miRNA-hubgene network. The network consisting of three circRNAs, one miRNAs, and five hubgenes. Green indicates circRNAs, yellow indicates miRNA, and blue indicates mRNA. Diamonds indicate circRNAs, rounded rectangles indicate miRNA, and ellipses indicate mRNA.

Discussion

circRNAs have become popular topics for research over previous years. circRNAs differ structurally to the well-known linear RNA due to the absence of the 3′ or 5′ polarities or polyadenylated tails (8). This increases their stability and protects them against degradation by the RNase-R enzyme (25). circRNA is abundant in eukaryotic cells where it is found to be structurally stable, highly conserved, and with tissue, timing and disease specificity (26,27). These characteristics make circRNAs potential biomarkers for several types of tumor (28,29). Although the exact functions of the majority of circRNAs remain unclear, a number of recent studies have suggested that circRNAs affect the initiation and development of different types of malignancies (30–32). Cao et al (32) reported that circRNA_100876 was highly expressed in esophageal squamous cell carcinoma. Furthermore, knockdown of circRNA_100876 was demonstrated to inhibit proliferation, migration and progression of the epithelial-mesenchymal transition (EMT). In addition, Chen et al (30) revealed that circPRMT5 promotes the EMT through sponging miR-30c to promote bladder cancer metastasis. Nevertheless, the exact role of circRNAs in PDAC remains undefined.

The present study first performed microarray analysis to identify DEcircRNAs in PDAC samples and associated normal samples. In order to increase the accuracy of the results, the present study used two online databases to predict their MREs. Only the genes that overlapped in both algorithms were selected as candidate miRNAs. The results identified 10 DEcircRNA-DEmiRNA pairs by intersecting with the DEmiRNAs. Using the same technique, the present study identified 120 DEmiRNA-DEmRNA pairs. Subsequently, a circRNA-miRNA-mRNA regulatory network was constructed, including four circRNAs, six miRNAs and 111 mRNAs.

Several studies have indicated that circRNA expression is dysregulated in PDAC, resulting in its pathogenesis and prognosis, and so it is considered to be a tumor-associated biomarker (33–35). Huang et al (34) demonstrated that hsa_circ_0000977 was significantly upregulated in PDAC tissues. Silencing hsa_circ_0000977 in vitro was revealed to decrease cell proliferation and induce cell cycle arrest. Furthermore, the authors identified hsa_circ_0000977 can regulate the expression of PLK1 by sponging hsa-miR-874-3p in the cytoplasm. Similarly, elevated circ_0030235 was observed in PDAC cell lines. Overexpression of circ_0030235 was associated with low survival rates and advanced clinicopathological features. As such, circ_0030235 is expected to be a prognostic factor for PDAC (35). In the present study, a total of four DEcircRNAs (hsa_circ_0006220, hsa_circ_0043278, hsa_circ_0001666 and hsa_circ_0092367) were identified to be involved in the ceRNA network. To the best of our knowledge, these four circRNAs have not previously been reported.

miRNAs are a large class of endogenous noncoding RNAs, 19–25 nucleotides in length, that are involved in the regulation of cell proliferation, differentiation, apoptosis and migration (36,37). Previous studies have investigated the binding of circRNAs to miRNAs and their interactions in pancreatic cancer (38,39). Hao et al (39) indicated that circ_0007534 promoted PDAC cell proliferation, apoptosis and invasion by sponging miR-892b and miR-625. An et al (38) reported that circZMYM2 participates in progression of pancreatic cancer by sponging miR-335-5p. In the present study, we predicted the correlation between four circRNAs and six miRNAs involved in the ceRNA network. Of these six miRNAs, four have previously been reported in PDAC (40–44). Zhu et al (42) demonstrated that miRNA-224 promotes PDAC cell proliferation and migration.

The GO and KEGG enrichment analyses suggest that these DEmRNAs have a significant effect on tumor-associated biological functions. Among the 12 total pathways, ‘focal adhesion’, ‘ECM-receptor interaction’, ‘microRNAs in cancer’ and ‘proteoglycans in cancer’ are associated with the progression of PDAC (45–48). The PPI network was established, including the six hubgenes (THBS1, FN1, TGFB2, ITGA1, ITGA3, and TIMP3), to further identify the key circRNAs participating in the regulatory network. Among these hubgenes, four genes (THBS1, TGFB2, ITGA1 and ITGA3) have been identified to play critical roles in the carcinogenesis and development of PDAC (49–52). However, to the best of our knowledge, the association between these genes and circRNAs has not yet been investigated. The present study established 16 circRNA-miRNA-hubgene axes in PDAC. Since the expression levels of hsa_circ_0092367 and TGFB2 were inconsistent, the hsa_circ_0092367/hsa-mir-375/TGFB2 regulatory module was excluded, leaving 15 circRNA-miRNA-hubgene axes. However, given that the results are based on bioinformatics, further studies are required in order to verify the potential role of the 15 axes in PDAC.

The present study presents several limitations. First, the number of samples is small. Future studies should include larger sample sizes in order to establish more accurate results. As neither of the two GSE datasets used in the present study provided patient survival information, the prognostic value of DEGs was not able to be investigated. In addition, although the patients' clinicopathological parameters were provided in each citation of their dataset, they do not provide the corresponding GEO sample ID number for each patient. Therefore, it was not possible to assess the association between DEGs and the clinicopathological parameters. Furthermore, the conclusions of the present study are only based on the current existing tools and databases, and thus lack in vitro analyses.

The present study successfully established a circRNA-associated ceRNA network and circRNA-miRNA-hubgenes regulatory modules via bioinformatics analysis. The results demonstrated that three significant circRNAs (hsa_circ_0006220, hsa_circ_0043278 and hsa_circ_0001666) may play important roles in PDAC progression, which provides new insight into the pathogenesis for patients with PDAC.