Genome-Wide Analysis of the Expression of Circular RNA Full-Length Transcripts and Construction of the circRNA-miRNA-mRNA Network in Cervical Cancer

Increasing evidence suggests that circular RNA (circRNA) plays an important role in tumorigenesis by regulating gene expression at the transcriptional and post-transcriptional levels. Alternative splicing events permit multiple transcript isoforms of circRNA to be produced; however, changes in the expression of circRNA full-length transcripts in cervical cancer remain unclear. Here, we systematically explored the dysregulation circRNA full-length transcripts and constructed an improved circRNA-miRNA-mRNA regulatory network to provide potential biomarkers and possible treatment targets in cervical cancer. We identified 9359 circular full-length transcripts from RNase R-treated RNA-seq data in cervical cancer, of which 353 circular full-length transcripts were significantly differentially expressed (DE) between the tumor and normal group. A total of 881 DE mRNA transcript isoforms were also identified from total RNA-seq data in cervical cancer, of which 421 (47.8%) transcript isoforms were up-regulated, and 460 (52.2%) transcript isoforms were down-regulated in tumor samples. Two circRNA-miRNA-mRNA competitively regulated networks, including 33 circRNA transcripts, 2 miRNAs, and 189 mRNA transcripts were constructed. Three genes (COPE, RAB3B, and TFPI) in the network were significantly associated with overall survival (P < 0.05), which indicated that these genes could act as prognostic biomarkers for patients with cervical cancer. Our study revealed genome-wide differential expression of full-length circRNA transcripts and constructed a more accurate circRNA-miRNA-mRNA network at the full-length transcript expression level in cervical cancer. CircRNA may thus be involved in the development of cervical cancer by regulating the expression of COPE, RAB3B, and TFPI. However, the specific regulatory mechanism in cervical cancer requires further study.


INTRODUCTION
Cervical cancer is the fourth most common type of cancer diagnosed in women (Bray et al., 2018) and the fourth main cause of cancer mortality in women; cervical cancer thus poses a serious threat to the health of women worldwide. Statistics from the Global Cancer Annual Report from 2018 revealed that there were approximately 570,000 new cases of cervical cancer and that approximately 311,000 cancer patients died of cervical cancer (Bray et al., 2018). Cervical cancer accounts for 80% of all cancers attributable to human papillomavirus (HPV) infection. According to statistics released by the China Cancer Center in 2015, the incidence of cervical cancer in China showed consistent annual increases from 2000 to 2011, and the age of onset shifted to younger ages . Cervical cancer is the sixth most common cause of female malignant tumors, accounting for 6.25% of cases. In China, there are approximately 100,000 new cases of cervical cancer and 30,000 deaths from cervical cancer annually. Current treatment methods for cervical cancer include surgery, radiotherapy, chemotherapy, and comprehensive treatment. Cervical adenocarcinoma is prone to lymphatic metastasis in the early stage, and the prognosis is relatively poor. Patients with advanced cervical cancer have lost the opportunity for surgical resection because of tumor metastasis and cannot bear the side effects of radiotherapy and chemotherapy (Cohen et al., 2019). There is thus a need to explore the internal mechanism of cervical cancer and identify new therapeutic targets.
With the development of next-generation sequencing technology and the integrated analysis and application of multiple omics data in recent years, non-coding RNA (ncRNA), such as microRNA (miRNA), long non-coding RNA (lncRNA), piRNA, and snoRNA, has been shown to be involved in the biological processes of various cancers. Researchers have also recently reported a new type of non-coding RNA with a special structure-circular RNA (circRNA) (Memczak et al., 2013;Salzman et al., 2013;Jeck and Sharpless, 2014;Li et al., 2015;Xu et al., 2017). CircRNAs are a large class of covalently closed transcripts that do not have a 5'-end cap structure and a 3'-end polyadenylated (PolyA) tail; thus, they are not easily degraded by exonucleases and are structurally stable and tissue specific. Previous studies have shown that circRNA has the potential to compete with endogenous RNA. CircRNA can act as a "sponge" of miRNAs by naturally isolating and competitively inhibiting the activity of miRNA, thereby indirectly regulating miRNA target genes (Hansen et al., 2013a;Liu et al., 2016;Xie et al., 2016;Yu et al., 2016;Han et al., 2017;Zhang et al., 2017). Like lncRNA and miRNA, circRNA also plays an important role in cervical cancer (Tornesello et al., 2020). For example, Tang et al. found that hsa_circ_0000515 can promote the development of cervical cancer (Tang et al., 2019). Specifically, hsa_circ_0000515 acted as a competitive endogenous RNA that specifically bound to miR-326 and inhibited its expression, thereby increasing the expression level of the target gene ETS transcription factor ELK1 (ELK1). The enhancement of ELK1 expression led to the enhancement of cervical cancer cell proliferation and invasion ability but repressed cell apoptosis and autophagy. Similarly, Song et al. found that reducing the expression of hsa_circRNA_101996 can significantly suppress the proliferation, migration, and invasion of cervical cancer cells (Song et al., 2019). Specifically, hsa_circRNA_101996 can act as a sponge of miR-8075 and down-regulate its expression level, thereby regulating miR-8075's target gene TPX2 and affecting the occurrence and development of cervical cancer. Zhao et al. also reported that the virus-encoded circular RNA-circE7 in HPV is preferentially located in the cytoplasm and that circE7 can be modified by N6-methyladenosine (m6A) and translated to produce E7 oncoprotein. The increased expression of E7 oncoprotein contributes to certain transforming activities of HPV cells . In sum, both circRNA in the human body and circRNA produced by viruses play an important role in the biological processes of cervical cancer.
A large number of circular RNAs have been identified in multiple species to date, and a small number of circRNAs have been shown to be involved in gene regulation , growth and development (Kristensen et al., 2018), and the innate immune response , and these circRNAs play an important role in various types of diseases, including cancer (Hirsch et al., 2017;Li P. et al., 2017;Piwecka et al., 2017;Yu et al., 2017;Chen et al., 2019;Vo et al., 2019). Existing circRNA-related research primarily focuses on using the total read counts aligned to back-splice junctions to estimate circRNA expression. However, the internal structure of alternative splicing circRNA transcripts and the expression of full-length circRNA transcripts in cancer remain unclear. Previous studies have suggested that alternative splicing events can cause circRNA to produce multiple transcript isoforms (Gao et al., 2016;Zhang X.-O. et al., 2016;Wu et al., 2019). Our understanding of the functions of specific circRNAs, including the evolutionary characteristics of circular RNAs between species and the differences in circular RNA transcripts between individuals and groups, remains limited. Moreover, the binding sites of circRNAs and protein-coding potentials also require information on the entire lengths of circRNAs. The diagram of this mechanism is shown in Figure 1A. For example, there are circRNAs derived from exon regions that contain five exons. Through different alternative splicing methods, circRNA can produce 4 full-length circRNA transcript types with the same back splicing site but with different numbers of exons ( Figure 1A a-d).
Here, we used a series of bioinformatics tools to systematically explore the differentially expressed (DE) full-length circRNA transcripts and more precisely predict the regulation target sites of circRNA in cervical cancer ( Figure 1B). Generally, our aim was to aid the search for potential biomarkers and possible treatments for cervical cancer. expression of full-length circRNA transcripts in cervical cancer. The public dataset included the RNA-seq data of 5 cervical cancer patients with primary tumors and matched adjacent normal samples, which were processed by ribonuclease R (RNase R) to remove linear RNAs and ribosomal RNA (rRNA).
To further verify and explore the function of differential full-length circRNA transcripts, we collected 4 pairs of samples from the primary tumors and adjacent normal tissue of cervical cancer patients from The First Affiliated Hospital with Nanjing Medical University to characterize differences in the abundances of mRNA transcripts. Sample collections were approved by the Human Research Ethics Committee of The First Affiliated Hospital with Nanjing Medical University (No.2020-SR-387), and permission from the cervical cancer patients was also obtained.

Total RNA Extraction and Sequencing
First, we cut the clinical samples into small pieces and added TRizol reagent (Life Technologies) to homogenize samples with a homogenizer. The total RNA extraction kit (Tiangen, China) was used to extract the total RNA from each sample per the manufacturer's instructions. Whole transcriptome libraries were constructed using the TruSeq Ribo Profile Library Prep Kit (Illumina, United States) per the manufacturer's protocol. Briefly, 10 mg of total RNA was depleted of rRNA with an Illumina Ribo-Zero Gold kit and purified for end repair and 5'-adaptor ligations. Next, reverse transcription was performed with random primers containing 3' adaptor sequences and randomized hexamers. Finally, cDNA fragments of 250-400 bp were purified and amplified by 15-20 cycles of PCR. The libraries were subjected to 150 nt paired-end sequencing with an Illumina Hiseq4000 system (BGI, China). At least 45 million read pairs were generated for each library, and adapters were removed with Cutadapt (Martin, 2011) to obtain clean reads.
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive  in National Genomics Data Center, Beijing Institute of Genomics (China National Center for Bioinformation), Chinese Academy of Sciences, under accession number HRA000394 that are publicly accessible at http://bigd.big.ac.cn/gsa-human.

Identification of Dysregulated circRNA Full-Length Transcripts and Linear mRNA Transcripts
Ten RNase R-treated RNA-seq samples (GSE107472) and 8 total RNA-seq samples were used to identify dysregulated circRNA full-length transcripts and linear mRNA transcripts, respectively. First, the raw sequence reads were cleaned by trimming the low-quality bases (Q < 20). The cleaned RNAseq reads were then mapped to the human reference genome (GRCh37/hg19, UCSC Genome Browser) by HISAT2 (Kim et al., 2015) (version 2.1.0) and BWA-MEM (Heng, 2013) software. These two tools are capable of detecting canonical splicing events. Next, CIRI2 (Gao et al., 2018) (version 2.0.6) was used to identify circRNAs in each sample by recognized the backsplicing reads (≥ 2). CIRI2 is an algorithm based on a multipleseed-matching strategy. Specifically, it extracts seed sequences in descending order of length for genomic regions with lower comparison quality and quickly matches them with the front and rear flanking genomic regions. At the same time, a maximum likelihood estimation model was established to determine the actual source of the seed sequence and eliminate interference from linear transcripts or splicing byproducts, thereby greatly improving the accuracy of circular RNA molecular recognition. Finally, StringTie (Pertea et al., 2015) (version 1.3.4d) and CircAST  were used to assemble linear mRNA transcripts and full-length circRNA transcripts, respectively. CircAST is a tool for assembling and quantifying full-length circRNA transcripts in RNase R-treated RNA-seq datasets based on graph theory and the minimum path coverage algorithm. To evaluate the relative expression of mRNA transcripts between tumor samples versus adjacent normal samples, DESeq2 (Love et al., 2014), an R package for differential expression analysis based on negative binomial generalized linear models, was used with the cutoff values [adjusted p-value < 0.05, | log2(fold change)| > 1]. Read counts mapped to the sequence of mRNA transcripts were used. Student's t-tests were used to determine DE full-length circRNA transcripts between tumor and adjacent normal samples [P < 0.05, | log2(fold change)| > 1]. The relative expression level of full-length circRNA transcripts in fragments per kilobase per million mapped reads (FPKM) was also used.

Construction of the circRNA-miRNA and ceRNA Network
To predict the relationship between the DE full-length circRNA transcripts and dysregulated miRNAs and the relationship between the dysregulated miRNAs and its target mRNA transcripts, miRanda (Betel et al., 2010) (August 2010 Release) and TargetScan (Agarwal et al., 2015) (Release 7.1) pipelines were used to predict the circRNA-miRNA interaction and circRNA-miRNA-mRNA competitively regulated network. The sequences and the annotation of miRNAs were obtained from the miRbase (Kozomara et al., 2019) database. In the miRanda pipeline, we set parameters with match scores higher than 150 and minimum free energy less than -20 kcal/mol to improve the reliability of our prediction. The visualization software Cytoscape (Shannon et al., 2003) (version 3.7.1) was used to display the above networks.

Integrated Functional Enrichment Analysis
The corresponding parental genes of full-length circRNA transcripts and the significantly dysregulated mRNA transcripts in the ceRNA network were both used for the functional enrichment analysis. Gene Ontology (GO) enrichment analysis (including Biological Processes, Cellular Components, and Molecular Function) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were performed using the WEB-based Gene Set Analysis Toolkit (WebGestalt) (Liao et al., 2019) with default parameters.

Evaluation of Overall Survival for Genes
We utilized the Gene Expression Profiling Interactive Analysis (GEPIA) database (Tang et al., 2017) to explore the survival curves. GEPIA captures the relationship between patient survival information and gene expression based on TCGA datasets. We selected cervical squamous cell carcinoma (CESC) and input the DE genes corresponding to DE mRNA transcripts in the network to explore the association between DE genes and overall survival (OS). The genes with P < 0.05 were considered as critical genes.

Identification of Full-Length circRNA Transcripts in Cervical Cancer
We first used CIRI2 to identify circRNAs in the tumors of 5 cervical cancer patients and adjacent normal samples. A total of 14,125 circRNAs with different back splicing sites were identified. According to the chromosomal location of circRNAs, we concluded that 82.8% of circRNAs were derived from exonic regions and 10.7% of circRNAs were derived from intronic regions; only 6.5% of circRNAs were derived from intergenic regions. Next, we used CircAST software to assemble and quantify the full-length transcripts of the circRNA derived from exons in cervical cancer. We successfully constructed a total of 9359 different full-length circRNA transcripts in primary tumor and adjacent normal tissues in cervical cancer. We found that most circRNAs (81%) had only one circular transcript isoform and that 14% and 2% of circRNAs had two and three circular transcript isoforms, respectively (Figure 2A). Approximately 2% of circRNAs contained more than 3 circular transcript isoforms. Although most circRNA have only one circRNA transcript isoform in cervical cancer, alternative splicing can occur inside circRNA, resulting in the formation of multiple different fulllength circRNA transcripts.
Next, we analyzed the internal composition of these fulllength circRNA transcripts in cervical cancer. Approximately 71% of the full-length circRNA transcripts had less than 4 internal exons, and nearly half of the circular transcripts were primarily composed of two exons (23.27%) or three exons (24.96%) ( Figure 2B). We also found that there were 32 (0.34%) fulllength circRNA transcript isoforms with more than 14 exons, of which circSYNE2 (chr14:64529469-64612941) had the most exons (composed of 35 exons). The circRNA transcripts derived from exons in cervical cancer were not particularly long, as most of the full-length circRNA transcript isoforms (5857/9359, 62.58%) were between 250-550 bp ( Figure 2C). Transcripts 300-350 bp in length accounted for the highest proportion (12.07%, 1130/9359), and only 0.48% of the transcripts were longer than 2500 bp. We statistically analyzed the number of circRNA transcripts and the number of genes in different chromosomes ( Figure 2D) and also compared the number of circular transcripts with the number of linear transcripts and their distribution in chromosomes (Supplementary Table S1). We found that the number of full-length circRNA transcripts was lower than the number of genes on each chromosome, but the lowest variance in gene and transcript numbers was observed on chromosome 18 (total of 252 circRNA transcripts and 289 genes), and the highest variance in gene and transcript numbers was on chromosome 19 (total of 205 circular transcripts and 1485 genes) (Figure 2D). Although the number of full-length circRNA transcripts on the chromosome was much lower than the number of corresponding linear transcripts of parental genes, the ratio of circular transcripts to linear transcripts was relatively stable and was mostly concentrated around 5% (Supplementary Table  S1). Chromosome 19 had the lowest proportion (1.61%), whereas chromosome 13 had the highest proportion (8.44%). Figure 2E shows the number of full-length circular transcripts produced by the parental genes of circRNAs. We found that most of the parental genes (77.8%) produced less than 4 full-length circular transcripts in cervical cancer, of which 44.1% genes produced only one full-length circular transcript, 21.7% genes produced two full-length circular transcripts, and 12.0% genes produced three full-length circular transcripts; 3.1% genes could produce more than 10 full-length circular transcripts. Figure 2F shows the length distribution of circ-exons contained in full-length circular transcripts composed of different numbers of exons. The length of circRNA transcript isoforms composed of one exon had a median value of 404.5 bp, which was significantly higher than that of circular transcripts composed of other numbers of exons (median value was approximately 248.1 bp, p < 0.05).

Analysis of DE Full-Length circRNA Transcripts in Cervical Cancer
We used CIRI2 software to identify the circRNAs in the tumors of 5 cervical cancer patients and adjacent normal samples. A total of 6183, 5140, 2532, 3346, and 3692 circRNAs were detected in each sample of the normal group, and 3961, 4102, 3483, 2927, and 2813 circRNAs were detected in the corresponding tumor group samples (Figure 3A). These circRNAs were supported by at least two reads that spanned the back-splicing site. We then used the CircAST software to assemble the full-length transcripts of circRNAs identified in each sample. A total of 4148, 3391, 1445, 1949, and 2308 circRNA transcripts were successfully assembled in each sample of the normal group; in the corresponding tumor group, 2456, 2535, 2223, 1661, and 1768 circRNA transcripts were successfully assembled ( Figure 3B). Thus, with the exception of patient No. 3, the number of circRNA and the number of fulllength circRNA transcripts in the tumor sample were lower than that of adjacent normal samples in cervical cancer. Individual differences between patients may explain the higher number of circRNA and full-length circular transcripts in the tumor sample of patient No. 3 relative to the adjacent normal sample.
Next, we compared the full-length transcripts of circRNA in each sample of the cancer group and the normal group. We found that 6991 and 5655 full-length circRNA transcripts were observed in the tumor group and adjacent normal group, respectively. Among them, 3704 (53.0%) specific full-length circular transcripts were identified in the cancer group, and 2368 (41.9%) specific full-length circular transcripts were successfully assembled in the normal group; there were 3287 full-length circular transcripts expressed in both groups ( Figure 3C). Next, we analyzed the parental genes corresponding to the full-length circRNA transcripts between the two groups ( Figure 3D). A total of 2718 and 2503 genes were found to correspond to the fulllength circRNA transcripts in the tumor and normal group, respectively, of which 842 (31.0%) parental genes were only expressed in the tumor group, 627 (25.0%) parental genes were only recognized in the adjacent normal group, and 1876 genes were expressed in both groups ( Figure 3D). To analyze differences in the expression level of the same full-length circRNA transcripts between the two groups, we normalized the expression level of each circular transcript obtained by the CircAST software according to formula 1 and used the relative expression level as estimated by FPKM for further analysis. Formula 1: Relative expression = log 2 expression + 1 Most of the data points were distributed on both sides of the blue line (Figure 3E), which means that the expression levels of most full-length circRNA transcripts were not significantly different between the cancer group and the adjacent normal group. The blue line represents the same transcript expression level between the two groups, and the red line represents a straight line fitted according to the distribution of the relative expression of all of the circular transcripts in the two groups. The red line was below the blue line, indicating that the differential full-length circular transcripts tended to be highly expressed in the adjacent normal group but weakly expressed in the cancer group. Thus, the number of down-regulated expressed full-length circRNA transcripts was higher than the number of up-regulated expressed full-length circRNA transcripts in tumors.
To further understand the role and potential functions of full-length circRNA transcripts in cervical cancer, we first defined full-length circRNA transcripts expressed in at least 4 samples as high-confidence full-length circRNA transcripts. After filtering, we obtained a total of 2133 high-confidence full-length circRNA transcript isoforms; next, Student's t-tests were used to determine the DE full-length circRNA transcripts in tumor and adjacent normal samples. Statistical analysis revealed that 353 full-length circular transcripts were significantly DE between the tumor and normal group (| log2 (fold change)| > 1, adjusted p-value < 0.05). Among them, 107 circular transcripts (30.3%) were significantly up-regulated, and 246 circular transcripts (69.7%) were significantly down-regulated in tumor samples ( Figure 3F). We also conducted a hierarchical cluster analysis to reveal the differential expression of full-length circRNA transcripts across tumor and normal samples. The tumor and the adjacent normal samples could be accurately classified into two different branches (Figure 4).

Construction of the circRNA-miRNA Regulatory Network in Cervical Cancer
Previous studies have shown that circRNA can compete with miRNA and regulate miRNA-mediated downstream targeted gene expression in a variety of diseases. However, previous studies have failed to reveal the internal composition of circRNA; consequently, our understanding of the interaction between miRNA and circRNA containing multiple exons remains incomplete.
Regulatory networks are an important tool for studying cancer, including cervical cancer (Song et al., 2018;Yi et al., 2019). We constructed a circRNA-miRNA co-expression regulatory network based on the DE full-length circRNA transcripts to further explore the potential biological functions of circRNA in cervical cancer. First, we employed a dataset from the HMDD (Human microRNA Disease Database) , which contained the experimentally DE miRNAs in cervical cancer. A total of 145 miRNAs were collected. Next, miRanda and TargetScan pipelines were used to predict the interaction between DE full-length circRNA transcripts and DE miRNAs in cervical cancer (see details in Materials and Methods). Finally, a total of 33 interaction pairs were obtained, and the circRNA-miRNA regulatory network was constructed (Figure 6). Our analysis demonstrated that miR-98-5p can interact with 22 circRNA transcripts, and miR-98-3p can interact with 11 circRNA transcripts. This finding suggested that the biological function of has-miR-98 during the development of cervical cancer is to act as a regulatory mediator of DE circRNA transcripts.

Construction of the circRNA-miRNA-mRNA Regulatory Network in Cervical Cancer
To further explore the downstream target genes regulated by miR-98-5p and miR-98-3p in cervical cancer, we collected 4 cervical cancer patient tumors and adjacent normal tissue samples for total RNA sequencing. HISAT2 and StringTie pipelines were used to determine mRNA transcript expression in cervical cancer (see details in "Materials and Methods"). A total of 169,739 linear mRNA transcripts were identified in all samples, of which 149,838 transcripts were annotated in gene transfer format file (GTF, GRCh37/hg19, downloaded from Ensembl Genome Browser) and 19,901 transcripts were newly discovered. Based on read counts, we used DESeq2 to identify the DE mRNA transcripts. We made volcano plots for mRNA transcripts in cervical cancer-paired samples (| log2 (fold change) | > 1, adjusted p-value < 0.05) (Supplementary Figure S1A). In total, 881 mRNA transcript isoforms were found to be DE in cervical cancer, of which 421 (47.8%) transcript isoforms were up-regulated and 460 (52.2%) transcript isoforms were FIGURE 6 | CircRNA-miRNA co-expression regulatory network in cervical cancer. The green and purple nodes represent miRNAs and circRNA full-length transcripts, respectively. Edge denotes their relationship.
down-regulated in tumor samples. Cluster analysis based on the differential expression of mRNA transcripts across tumor and normal samples showed that the tumor and the adjacent normal samples could be accurately classified into two different branches (Supplementary Figure S1B).
We then predicted the interaction between DE miRNAs and DE mRNA transcript isoforms in cervical cancer, and a total of 226 miRNA-mRNA interaction pairs were obtained. By combining the interaction pairs in the circRNA-miRNA regulatory network and positive relationship with circRNA and mRNA, two circRNA-miRNA-mRNA competitive regulatory networks with different regulated trends were then constructed (Figures 7A,B). One network included 11 up-regulated full-length circRNA transcripts, 2 miRNAs, and 82 up-regulated mRNA transcripts ( Figure 7A). Another network included 22 down-regulated full-length circRNA transcripts, 2 miRNAs, and 107 down-regulated mRNA transcripts ( Figure 7B).
To explore the function of aberrantly expressed mRNA transcripts in the ceRNA network, the annotated genes corresponding to the mRNA transcripts were analyzed by WebGestalt. We analyzed the functions of a total of 105 DE genes ( Figure 7C). GO Biological Process term enrichment analysis showed that the target genes in the ceRNA network were enriched in muscle structure development, neural tube development, regulation of translation, protein-DNA complex subunit organization, and muscle contraction (P < 0.05). The GO Molecular Function term analysis revealed that the DE genes were enriched in ATPase activity, histone binding, pyrophosphatase activity, and hydrolase activity (P < 0.05).

Overall Survival of Target Genes in the Network
To explore the relationship between our identified target genes and clinical features, we utilized the Gene Expression Profiling Interactive Analysis (GEPIA) database. A total of 3 out of 105 DE genes from the ceRNA network were significantly correlated with CESC OS (P < 0.05) (Figure 8). The expression level of COPE was negatively correlated to OS (i.e., survival was longer for patients with higher expression levels in CESC). The expression level of RAB3B and TFPI were positively correlated to OS (i.e., survival was shorter for patients with higher expression levels in CESC).

DISCUSSION
Cervical cancer is the fourth most common cancer and the fourth main cause of cancer mortality in women worldwide (Bray et al., 2018). Although previous studies have reported that HPV plays a critical role in the development of cervical cancer (Walboomers et al., 1999), the specific molecular mechanism following viral infection remains unclear. Previous studies have shown that non-coding RNA, such as miRNAs (Liu et al., 2018;Pardini et al., 2018) and lncRNAs (Bhan et al., 2017;Peng et al., 2017;Luan and Wang, 2018), participate in various biological processes of cervical cancer. Recent work shows that circular  RNA (circRNA), a new type of non-coding RNA with a special structure, also plays an important role in cancer. The functions of most circRNAs are thought to include acting as miRNA sponges (Hansen et al., 2013a,b;Memczak et al., 2013), interacting with RNA-binding proteins (Ashwal-Fluss et al., 2014;Dudekula et al., 2016), and translating peptides (Legnini et al., 2017;Pamudurti et al., 2017;Zhang et al., 2018). Previous studies have shown that circRNA has the potential to compete with endogenous RNA in cervical cancer. However, as in RNAs formed by linear splicing, alternative splicing events also occur in circRNAs and can lead to the production of multiple transcripts (Gao et al., 2016;Zhang X.-O. et al., 2016;Wu et al., 2019). Few studies have examined the expression of full-length circRNA transcripts in cervical cancer. To address this knowledge gap, we systematically explored the DE full-length circRNA transcripts and more precisely predicted the regulation target sites of circRNA in cervical cancer.
We characterized a total of 14,125 circRNAs in cervical cancer and successfully assembled 9,359 full-length transcripts of circRNA derived from exons. Furthermore, we found that 81% of circRNAs had only one circular transcript isoform, indicating that most exon-derived circRNAs were not derived from alternative splicing events in cervical cancer. Nevertheless, 2% of circRNAs contained more than 3 circular transcript isoforms. The length of full-length circRNA transcripts formed by one exon (median value was 404.5 bp) was significantly higher than that of circular transcripts composed of other numbers of exons (median value was 248.1 bp, p < 0.05). Furthermore, 353 full-length circRNA transcripts and 881 lineal mRNA transcripts were identified to be DE in tumors and adjacent normal tissues. The circRNA-miRNA-mRNA regulatory network was also established in our study. The network included 33 DE circRNA transcripts, 2 DE miRNAs, and 189 DE mRNA transcripts in cervical cancer. Analysis of a TCGA dataset revealed that a total of 3 target genes (COPE, RAB3B, and TFPI) in the regulatory network were correlated to survival curves with clinically significant outcomes in cervical cancer.
Previous studies have reported the overexpression of RAB3B in various types of cancer, such as laryngeal squamous cell carcinoma  and lung squamous cell carcinoma (Zhang C. et al., 2016). In addition, overexpression of RAB3B has been associated with the prognosis, clinicopathological features, and the proliferation ability of tumor cells Ye et al., 2014). All of these results thus indicate that RAB3B is an oncogene. By further analysis, there are 33 full-length circRNA transcripts associated with RAB3B. We suspect that circRNA may be involved in the development of cervical cancer by regulating the expression of RAB3B; however, additional experiments of regulation pairs are needed to confirm this speculation.
In sum, our work provided a detailed list of changes in the expression of full-length circRNA transcripts in cervical cancer. This list of circRNAs may potentially serve as a novel biomarker or therapeutic target. The ceRNA network also provides potential targets for future functional studies of circRNAs in cervical cancer.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of The First Affiliated Hospital with Nanjing Medical University (No. 2020-SR-387). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
TX conceived and designed the study, collected the public data, performed the analysis, and wrote the manuscript. XS assisted with manuscript review and revision. YW provided the technical support and performed the analysis. SF and PH collected the clinical samples and assisted with manuscript review. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
This paper is recommended by the 5th Computational Bioinformatics Conference.