Preprocessing of expression profile data
RNA expression profile datasets that contained no less than 10 pairs of human pancreatic cancer and normal control samples were retrieved from NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/) and EBI ArrayExpress (Https://www.ebi.ac.uk/arrayexpress/). A total of 2 circRNA, 5 miRNA and 7 mRNA expression profile chip datasets were thus obtained, as shown in Table 1.
Table 1
RNA datasets used for analysis, including 2 circRNA datasets, 5 miRNA datasets and 7 mRNA datasets.
| Dataset | Platform | Normal | Disease | Total samples |
circRNA | GSE79634 | GPL19978 | 20 | 20 | 40 |
GSE69362 | GPL19978 | 6 | 6 | 12 |
miRNA | E-MTAB-753 | A-AFFY-184 | 17 | 17 | 34 |
GSE71533 | GPL18058 | 16 | 36 | 52 |
GSE43796 | GPL15159 | 5 | 6 | 11 |
GSE41369 | GPL16142 | 9 | 9 | 18 |
GSE32678 | GPL7723 | 7 | 25 | 32 |
mRNA | GSE62452 | GPL6244 | 61 | 69 | 130 |
GSE71989 | GPL570 | 8 | 14 | 22 |
GSE55643 | GPL6480 | 8 | 45 | 53 |
GSE28735 | GPL6244 | 45 | 45 | 90 |
GSE32676 | GPL570 | 7 | 25 | 32 |
GSE16515 | GPL570 | 16 | 36 | 52 |
GSE15471 | GPL570 | 39 | 39 | 78 |
There were two original types of data: TXT and CEL. Log2 transformation was performed using the Limma package (version 3.32.5, http://bioconductor.org/packages/release/bioc/html/limma.html)[3] in R3.4.1 to transform the TXT format data from a skewed distribution to an approximately normal distribution. The CEL format data were normalized by the Oligo package (version 1.41.1, http://www.bioconductor.org/packages/release/bioc/html/oligo.html)[4] in R3.4.1.
Screening Of Differentially Expressed Genes
Flow chat of the bioinformatics analysis was shown in Fig. 1. Significantly differentially expressed circRNAs in the two circRNA datasets were screened with the Limma package[3], with a p value less than 0.05 and | logFC |> 0.585 (multiplied by 1.5) as the thresholds. All the circRNA names were converted into those commonly used in the circRNA database (Http://www.circbase.org/)[5](see Table 1). The differentially expressed miRNAs and mRNAs were screened by MetaDE (version 1.0.5, https://CRAN.R_project.org/package=MetaDE)[6] in R3.4.1, with tau2 = 0 and Q Pval > 0.05 as the homogeneity test parameters and FDR < 0.05 as the threshold for significant expression differences between genomes. The MetaQC package (version 0.1.13, https://cran.r_project.org/web/pathe ckages/MetaQC/index.html)[7] in R3.4.1 combined with standardized mean rank (SMR) and principal component analysis (PCA) was used as an objective data quality control measure and confirmed that the data in the 5 miRNA datasets and 7 mRNA datasets were well distributed and suitable for further analysis (supplementary Fig. 1). Bidirectional hierarchical clustering of the above differentially expressed RNAs was performed based on Euclidean distance with the Pheatmap package (version 1.0.8, https://cran.r_project.org/package=pheatmap)[8] in R3.4.1.
Regulatory Network Construction
The sequences of differentially expressed circRNAs were downloaded from circBase, and those of miRNAs were downloaded from miRBase Release 21 (http://www.mirbase.org/)[9]. CircRNA-miRNA interactions were predicted through the miRanda program (http://cbio.mskcc.org/miRNA2003/miranda.html)[10]. The target genes of differentially expressed miRNAs were predicted with TargetScan (Release 7.1, http://www.targetscan.org/vert_71/)[11], MicroCosm (version 5, https: // www .ebi.ac.uk / enright-srv / microcosm / htdocs / targets / v5 /)[12], miRecords (http://mirecords.umn.edu/miRecords)[13], and miRNAMap (version 2, http: //mirnamap.mbc.nctu .edu.tw /)[14]. Genes that were identified by at least two of the above databases and were negatively correlated with the corresponding miRNAs were considered potential target genes, and thus, the miRNA-mRNA regulatory relationship was confirmed. The target gene interactions were also explored with BioGRID (version 3.4.153, http://thebiogrid.org/)[15, 16], HPRD Release 9 (http://www.hprd.org/) [17] and STRING (version 10.5, https://string-db.org/)[18].
The data of pancreatic cancer samples with survival information were obtained from the TCGA database (https://gdc_portal.nci.nih.gov/). The miRNAs and mRNAs in the above regulatory relationships that were significantly related to prognosis were screened by single-factor Cox regression analysis with the survival package (version: 2.40-1, https://cran.r-project) in the R3.4.1 language .org / package = survival)[19] to construct a ceRNA network related to prognosis. The miRNA-associated pathways involved in the ceRNA network were downloaded from miRWalk 2.0 (http://zmf.umm.uni_heidelberg.de/apps/zmf/mirwalk2/index.html)[20].
Finally, the regulatory relationships were visualized with Cytoscape 3.3 (http://www.cytoscape.org/). GO and KEGG pathway analyses were performed by DAVID 6.8 (https://david.ncifcrf.gov/)[21, 22] online, and a p value less than 0.05 was considered the significance threshold.
Verification of the 8 most differentially expressed circRNAs in pancreatic cancer cell lines
The pancreatic cancer cell lines AsPC-1, PANC-1, BxPC-3, CFPAC-1, and HPAF-II and the normal human pancreatic cell line HPDE6-C7 were purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences, thawed and cultured in DMEM (Gibco, C11995500BT) or RPMI-1640 medium (Gibco, C11875500BT) containing 10% fetal bovine serum, 1% penicillin/streptomycin (complete medium, Gibco) with or without 1% glutamine and 1% sodium pyruvate (Gibco) according to the manufacturer’s protocols. The cells were cultured in the recommended incubator at 37 °C supplied with 5% CO2. Total RNA was isolated from the above cell lines with RNAiso Plus (TRIzol, TaKaRa) according to the manufacturer’s instructions. The RNA samples (0.5 µg each) were reverse transcribed into cDNA by the reverse transcription kit PrimeScript™ RT Master Mix (Perfect Real-Time, TaKaRa, Dalian, China). The relative expression levels of 8 differentially expressed circRNAs, hsa_circ_0079385, hsa_circ_0002191, hsa_circ_0006502, hsa_circ_0092367, hsa_circ_0000919, hsa_circ_0013587, hsa_circ_0000376 and hsa_circ_0008240, were determined by qRT-PCR by using Power SYBR Green PCR Master (Thermo, 4367659), and the primer sequence of each circular RNA is presented in supplementary table 1. All experiments were repeated 3 times in every cell line. Amplification plots and dissociation curves were generated to confirm the quality of the data (supplementary Fig. 2).
Functional studies of hsa_circ_0006502 in CFPAC-1 and AsPC-1 cell lines
The hsa_circ_0006502 cDNA with BamHI (GGATCC) and EcoRI (GAATTC) restriction sites was synthesized by Yuanzi Biotech (Shanghai, China) and inserted into the pLCDH_ciR vector to obtain the hsa_circ_0006502-overexpressing plasmid pLCDH_ciR_hsa_circ_0006502. Lentiviral particles were obtained by cotransfecting HEK293T cells with pLCDH_ciR_hsa_circ_0006502 and packaging plasmids (pMD2 and pAX2). CFPAC-1 and AsPC-1 cells were infected with concentrated pLCDH_ciR_hsa_circ_0006502 lentiviral particles and then screened by puromycin. qRT-PCR was performed 72 h later to detect the expression level of hsa_circ_0006502. Cell proliferation was tested with a cell counting kit-8 (CCK-8) assay at 24 h, 48 h, and 72 h. The apoptosis assay was performed with flow cytometry and the migratory and invasive abilities were assessed in CFPAC-1 and AsPC-1 cells 48 h after transduction of hsa_circ_0006502 by Transwell assays according to the manufacturer’s protocols.
Statistical analysis
All results are presented as the mean ± standard deviation (SD). Statistical analyses were performed with GraphPad Prism 7 (GraphPad Software, San Diego, CA). Student’s t test and one-way ANOVA were utilized to analyze significant differences, and P < 0.05 and P < 0.01 were considered to represent statistically significant and extremely significant differences, respectively.