Circular RNA profile in gliomas revealed by identification tool UROBORUS

Abstract

Recent evidence suggests that many endogenous circular RNAs (circRNAs) may play roles in biological processes. However, the expression patterns and functions of circRNAs in human diseases are not well understood. Computationally identifying circRNAs from total RNA-seq data is a primary step in studying their expression pattern and biological roles. In this work, we have developed a computational pipeline named UROBORUS to detect circRNAs in total RNA-seq data. By applying UROBORUS to RNA-seq data from 46 gliomas and normal brain samples, we detected thousands of circRNAs supported by at least two read counts, followed by successful experimental validation on 24 circRNAs from the randomly selected 27 circRNAs. UROBORUS is an efficient tool that can detect circRNAs with low expression levels in total RNA-seq without RNase R treatment. The circRNAs expression profiling revealed more than 476 circular RNAs differentially expressed in control brain tissues and gliomas. Together with parental gene expression, we found that circRNA and its parental gene have diversified expression patterns in gliomas and control brain tissues. This study establishes an efficient and sensitive approach for predicting circRNAs using total RNA-seq data. The UROBORUS pipeline can be accessed freely for non-commercial purposes at http://uroborus.openbioinformatics.org/.

Figure 1.

The UROBORUS pipeline for identifying circRNA based on total RNA-seq. The artificial paired-end seed (20 bp) was first extracted from two ends of reads in an unmapped.sam file, and then aligned to the reference genome. The results would have two cases of reads spanning the spliced site: balanced mapped junction (BMJ) reads and unbalanced mapped junction (UMJ) reads. The UROBORUS pipeline designed algorithm to deal with BMJ and UMJ reads, and detects more circRNA supported reads.

Figure 2.

CircRNAs supported by BMJ and UMJ reads. UMJ reads are neglected by other approaches, leading to underestimate of the circRNA expression level or complete missing of those only UMJ supported circRNAs. (A) UROBORUS can identify both BMJ and UMJ reads supported circular ERC1, while find_circ and CIRCexplorer can only detect fraction of BMJ reads of circular ERC1, leading to expression level underestimation of circular ERC1. Sanger sequencing confirmed ERC1 circRNA. (B) Unlike find_circ and CIRCexplorer, one of the most powerful features in UROBORUS is its ability to identify those only UMJ reads supported circular RNA such as CORC1C, which were further confirmed by Sanger sequencing.

Figure 3.

Comparison of three pipelines. (A) CircRNAs (supported by 2 reads) identified by the three methods were compared using a Venn Diagram; (B) UROBORUS and CIRCexplorer identify more common circRNAs as shown by the factions when the expression level of circRNAs increases; (C) The comparison of the three methods in language used, depending tools, memory and running time.

Figure 4.

Primer design, PCR and Sanger sequencing validation for circ_ERC1. Two pairs of primers (convergent or divergent primers) were designed for each candidate circRNA (primer information in Supplementary Table S2). The junction region of circular RNA sequenced by Sanger sequencing validated the circular ERC1 formation; Gel electrophoresis validated the existence of circular ERC1.

Figure 5.

Genomic features of circRNA in gliomas and normal tissue. (A) There is no significant difference in the exon number of highly expressed circRNAs between normal and glioma tissues; (B) The genomic distance of the back-splicing site in most circRNAs is within 50 kb both in normal and glioma tissue, with only a few circRNAs spanning 100–300 kb; (C) The number of alternative circularization numbers in normal tissue is significantly higher than in glioma tissue; (D) Exon length of circRNA in all samples; (E) RIMS circRNA has more alternative circularization in normal tissue than in glioma tissue; (F) VCAN circRNA isoforms are similar in normal and glioma tissue.

Figure 6.

CircRNA distribution in normal tissue and gliomas. The upper panel shows the circRNA distribution in different chromosomes in normal tissue and gliomas. The circRNA number in GBM was significantly lower compared to normal tissues or oligodendroblastoma (Wilcoxon rank-sum test, P-value = 1.944e-09, 0.0002117 separately); no significant difference in the circRNA number between normal tissue and oligodendroblastoma (Wilcoxon rank-sum test, P-value = 0.1516). The lower panel shows the circRNA number in each chromosomes in normal tissue and gliomas.

Figure 7.

Differential circRNA expression and their parental mRNA expression in normal tissue and gliomas. (A) The heat map of expressed circRNAs in normal, oligodendroblastoma, and GBM (RPM > 0.1); (B) The heat map of parental gene related to circRNAs in normal, oligodendroblastoma and GBM.

Figure 8.

Pearson correlation coefficient between circRNAs and their parental mRNA. (A) The genes were arranged from left to right according to the RSME value from lower to higher of PCC among normal, GBM, and oligodendroblastoma. Thus, the left 185 circRNAs and their parental mRNAs have the similar correlation pattern in normal or gliomas (RMSE < 0.15). The correlations of the remaining about 387 circRNAs with their parental mRNAs have different patterns in normal and gliomas; (B) The PCC Frequency distribution in normal and gliomas.