Fifteen pairs (10 for RT-PCR and 5 for RNA sequencing) of non-neoplastic gastric tissues (NT) and GCs tissues from Shenzhen Second People’s Hospital (The First affiliated Hospital of Shenzhen University) were examined in the study. All tissues received no radiotherapy or chemotherapy before surgery and stored in RNAlater immediately after surgery. All patients provided written informed consent, approved by the ethics committee of Shenzhen University School of medicine.

Total RNA was extracted from 5 paired GC samples. Ribosomal RNA was digested using ribo zero kit and linear RNA was removed using rnase R. CircRNA was fragmented and cDNA was synthesized with six base random primers and purified. Strand-specific library is prepared via incorporating a chemical label deoxy-UTP (dUTP) during production of the second-strand cDNA. The second-strand cDNA was specifically digested by UNG enzyme. The constructed library was analyzed by Agilent 2,100 Bioanalyzer and was sequenced on an Illumina Hiseq 2,500 (Chi Biotech, Shenzhen).

We used 10 samples from 5 gastric cancer patients to investigate potential circRNA biomarker. From each patients, two samples were collected: one from the cancer tissue and other from the adjacent normal tissue. After getting the raw data in fastq format, the quality control analysis was performed. The quality control was done using the NGSQCToolkit (Patel and Jain, 2012) and the quality score 20 was used as the cutoff point. Filtering/trimming was carried out if the quality of the reads failed to reach the cutoff point. The fastq reads were aligned to the human reference genome (version GRCh38) using the BWA (Li and Durbin, 2009) aligner. CircRNAs were identified using the software CIRI (version 2) (Gao et al., 2015). The identified circRNAs were then annotated with the gene annotation file corresponding to the reference genome and the full length circRNA sequences were extracted. The full length circRNA sequences for all circRNAs were considered as the reference genome and the fastq reads were mapped using the bowtie2 (Langmead and Salzberg, 2012) aligner. Then count data was generated using bedtools (Quinlan and Hall, 2010) multiBamCov with the output of bowtie2 (converted to bam, sorted and indexed). R package DESeq (Anders and Huber, 2010) and DESeq2 (Love et al., 2014) were used to identify differentially expressed (DE) circRNAs for individual patient and combining all patients together respectively. |Log2FoldChange|>1 and p-value<0.05 were considered as the cut-off for defining significant DE circRNAs. For the DE analysis of combination of all patients, we also considered false discovery rate (FDR) at 10%.

According to the manufacturer’s protocol, total RNA was extracted using Trizol reagent (251,808, Invitrogen). Reverse transcription and real-time PCR were performed using GoScript™ Reverse Transcription Mix (A2800 and A6002, Promega). Divergent primers were designed. All primers were synthesized by Sangon Biotech. Results were normalized to the 18S mRNA in each sample. We used the following primer sequences.

We were interested with those circRNAs which were DE in four or more patients (out of five patients). Then, some top DE circRNAs were selected for miRNA interaction analysis. The circRNA-miRNA interaction was predicted using miRanda (Enright et al., 2003) software. The miRNA sequences were downloaded from mirbase (Kozomara et al., 2019). From the circRNA-miRNA interaction, a sub-network was constructed keeping only the gastric cancer related miRNAs. Then a circRNA-miRNA-gene network was constructed using Cytoscape (Shannon et al., 2003). The genes were the host genes of the circRNAs. Finally, a sub-network was constructed with the top 5 hub circRNAs.

Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/cgi/input.pl) is a database for obtaining the protein-protein interaction of the provided genes. The host genes of the DE circRNAs were mapped to STRING database, and the protein-protein interaction among the genes was obtained. Then, a PPI network was constructed using the software Cytoscape, and from the network, the top 50 hub genes were selected for GO term and KEGG pathway enrichment analyses. With the selected hub genes, GO term and KEGG pathway enrichment analyses were performed using DAVID (Huang et al., 2009) to know the function of the circRNAs. The threshold p-value<0.05 was used for the significance of the enrichment analysis.

A total of 45,783 circRNAs were identified in all samples and 21,652 of them were co-expressed in cancer and normal samples (Figure 1A). Among the total circRNAs, 79% (36,218/45,783) were exonic, 1% (462/45,783) were intronic, 4% (1848/45,783) were intergenic, 15% (6,744/45,783) were sense overlapping and 1% (511/45,783) were antisense (Figure 1C). The length distribution of the circRNAs was shown in Figure 1B and observed that the length of 74% (44 + 30) of the circRNAs was less than 1,000 nucleotides (nt). The length of only 13% circRNAs was greater than 2000 nt. The number of back-spliced reads for the 97% of the circRNAs were less than 20 (Figure 1D). The average number of back spliced reads per circRNA was 6 whereas the minimum and the maximum number of back spliced reads were 2 and 898 respectively. The distribution of circRNAs in different chromosomes were heterogeneous, most circRNAs were originated from chromosome NC_000,001.11 (chromosome 1) (Figure 1E). The average number of circRNAs per chromosome was 1,405. The length of the circRNAs ranged from 48 to 199,579 nt while median and average lengths were 563 nt and 3,042 nt respectively (Table 1).

The expression values of the 10 samples from the five patients were shown in Figure 2A. The expression values were almost homogeneous across the samples. We also performed the Principal Component Analysis (PCA) of the expression values and the loading of the PC1 and PC2 were shown in Figure 2B. The distribution of circRNAs across the number of exons was shown in Figure 2C. Majority of the circRNAs were composed of 2 or 3 exons. About 76% of the circRNAs were derived from less than or equal to 5 exons. The distribution of the number of circRNAs per gene was given in Figure 2D. Total 45,783 circRNAs were originated from 8,095 genes. Majority of the cases (about 30%), the number of circRNA per gene was 1. About 55% of the cases, the number of circRNAs per gene was less than or equal to 3.

We identified the differentially expressed circRNAs between two groups (normal vs cancer). The R package DESeq was used to find the differentially expressed circRNAs. We collected 10 samples from five patients: for each patients two samples were collected one from cancer tissue and other from the adjacent normal tissue. We also performed DE analysis between normal and cancer groups for each patient. The volcano plot was used to find the target circRNAs with the thresholds p-value<0.05 and |log2FoldChange| >1 (Figure 3). The number of DE circRNAs for all samples were given in Table 2. Out of 45,783 total circRNAs, 478 were differentially expressed between normal vs cancer groups. From Table 2, we observed that about half of the DE circRNAs were upregulated, whereas the other half were downregulated. The expression patterns of these circRNAs were shown in the heatmap (Figure 4). The dendrograms showed that the cancer and normal samples were clearly distinguishable. Firstly, we identified DE circRNAs for all patients (5 normal samples vs 5 cancer samples) and performed DE analysis for each individual patient and isolated those circRNAs which were DE in four or more patients. 4 upregulated and 7 downregulated circRNAs were found to be DE in more than four patients. We selected 8 circRNAs (4 upregulated and 4 downregulated) from the DE analysis of individual patient. These 8 circRNAs were not only DE patient-wise but also DE while performed DE analysis combining all patients (5 normal samples vs 5 cancer samples). Finally, these 8 circRNAs were selected for circRNA-miRNA interaction analysis. The details of the 8 circRNAs were given in Table 3.

We got a total of 478 circRNAs as DE between normal vs cancer groups. Among these circRNAs, we investigated those which were DE in more number of patients. We chose those circRNAs for interaction analysis, which were DE in four or more patients (out of five patients). We selected the 4 upregulated and 4 downregulated (top 4 out of 7 downregulated) circRNAs for miRNA interaction analysis. We got a total of 1990 interactions with 8 circRNAs and 1,278 miRNAs. The miRNA sequences were downloaded from the mirbase database. From the online search we found 64 GC related miRNAs out of 1,278 miRNAs in the network. Next, we extracted the interactions of these 64 miRNAs and obtained a total of 81 interactions. Then, a circRNA-miRNA-gene network was constructed using Cytoscape. Here, the genes were the host genes of the circRNAs. Finally, we constructed a subnetwork with top 5 hub circRNAs (Figure 5). From the network, we observed that two genes were functionally related to GC. Two circRNAs circCEACAM5 and circCOL1A1 were interacted with GC related miRNAs and their host genes were also related to GC. This indicated that these two circRNAs might be potential biomarker for GC. The detail list of top 5 hub circRNAs, their interacted GC related miRNAs and their host genes were given in Table 4.

We performed the differential expression analysis of the miRNAs in gastric cancer tissue and the adjacent normal tissue using the GEO dataset (GSE158315). R package limma (Ritchie et al., 2015) was used to find the DE miRNAs. The expression patterns of the DE miRNAs were shown in Figure 6A. The dendrograms showed that the cancer and normal samples were clearly distinguishable. The volcano plot (Figure 6B) was used to find the DE miRNAs with the threshold p-value<0.05 and |log2FoldChange|>1. We found a total of 176 DE miRNAs of which 60 were upregulated and 116 were downregulated. The number of targeted miRNAs for circRNAs circCEACAM5 and circCOL1A1 were 9 and 5 respectively. The unique targeted miRNAs for these two circRNAs were 13 (Figure 6C). MiRNAs hsa-miR-634 and hsa-miR-429 were common between the 9 targeted miRNAs of circCEACAM5 and the 176 DE miRNAs. Again, miRNAs hsa-miR-634, hsa-miR-21-3p and hsa-miR-145-5p were common between the 5 targeted miRNAs of circCOL1A1 and the 176 DE miRNAs. MiRNAs hsa-miR-429, hsa-miR-21-3p and hsa-miR-145-5p were found to be significantly upregulated in GC (Figures 6D–F). Therefore, miRNA hsa-miR-429 might be the potential target of the circRNA circCEACAM5, and miRNAs hsa-miR-21-3p and hsa-miR-145-5p might be the potential targets of the circRNA circCOL1A1.

We selected top 8 circRNAs (4 upregulated and 4 downregulated) for the downstream analysis. These 8 circRNAs were originated from 6 host genes. The expression of these 6 host genes were checked by the TCGA data of 408 gastric tumor samples and 36 normal samples obtained from Gepia (version 2) (Tang et al., 2019). We found that 4 genes were upregulated in the GC tissue compared to the normal tissue Figures 7A–D. In addition 2 genes were downregulated in GC compared to the normal samples Figures 7E,F. The 4 upregulated host genes might be the potential targets of their respective circRNAs.

We also performed the receiver operating characteristic (ROC) curve analysis for our two candidate circRNA biomarkers circCEACAM5 and circCOL1A1 (Figures 8A,B). The area under the ROC curves (AUC) for these two circRNAs were significant (p-value<0.05). The AUC values 0.96 and 0.94 represented respectively the accuracy of circCEACAM5 and circCOL1A1 in distinguishing GC and normal patients.

We constructed the PPI network from STRING database for the host genes of the DE circRNAs and selected the top 50 hub genes using cytoscape. To explore the functions of circRNAs, we performed the GO and KEGG pathway analyses for the top 50 hub genes of the DE circRNAs from the PPI network. We used the cut-off p < 0.05 for finding the significant pathways and GO terms. The GO term biological process analyses showed that some of the genes were enriched cell-cell adhesion, extracellular matrix organization, DNA repair, collagen catabolic process, collagen fibril organization, response to virus etc. The GO term molecular function analyses showed that some host genes were able to bind many kinds of molecules including ATP, protein, platelet-derived growth factor, DNA, microtubule, chromatin DNA, double-stranded RNA etc. Some other genes were enriched in cadherin binding involved in cell-cell adhesion, extracellular matrix structural constituent, microtubule motor activity etc. The GO term cellular component analyses showed that some genes were enriched in nucleus, nucleoplasm, cytosol, cytoplasm, mitotic spindle, basement membrane, mitotic spindle etc. The KEGG pathway analysis indicated that some host genes were enriched in several significant pathways including ECM-receptor interaction, PI3K-Akt signaling pathway, Focal adhesion, Cell cycle etc. The results of the significant GO term (top 30) and KEGG pathway analyses were provided in Tables 5, 6 respectively.

We selected 8 circRNAs (4 upregulated and 4 downregulated) for miRNA interaction analysis. These 8 circRNAs were DE in four (out of five) or more patients. Among these 8 circRNAs, we selected 5 circRNAs (3 upregulated and 2 downregulated) for experimental validation by qRT-PCR. The expression of these 5 circRNAs were shown in Figure 9. From Figure 9, we observed that the expression of circCEACAM5, circMISP, circCOL1A1 was upregulated in GC, which is consistent with the sequencing result. The expression of circLIPF and circPGC was decreased in GC in line with the sequencing result. Hence, the computational results for these 5 circRNAs were validated by the qRT_PCR experiments.