The Landscape of Circular RNA in Cancer

Abstract

Circular RNAs (circRNAs) are an intriguing class of RNA due to their covalently closed structure, high stability, and implicated roles in gene regulation. Here, we used an exome capture RNA sequencing protocol to detect and characterize circRNAs across >2,000 cancer samples. When compared against Ribo-Zero and RNase R, capture sequencing significantly enhanced the enrichment of circRNAs and preserved accurate circular-to-linear ratios. Using capture sequencing, we built the most comprehensive catalog of circRNA species to date: MiOncoCirc, the first database to be composed primarily of circRNAs directly detected in tumor tissues. Using MiOncoCirc, we identified candidate circRNAs to serve as biomarkers for prostate cancer and were able to detect circRNAs in urine. We further detected a novel class of circular transcripts, termed read-through circRNAs, that involved exons originating from different genes. MiOncoCirc will serve as a valuable resource for the development of circRNAs as diagnostic or therapeutic targets across cancer types.

Keywords: biomarkers; cancer; circRNA; circRNA database; exome capture sequencing; non-coding RNA; read-through transcripts.

Figure 1.. Validation of the exome capture RNA-seq method for circRNA detection.

A) When compared to matched Ribo-Zero samples, capture transcriptome sequencing consistently detected more circular RNA (circRNA) in six paired libraries from clinical samples and cell lines. Capt: exome capture RNA-seq; Ribo: Ribo-Zero. Details about sequencing depths and the number of circRNAs stratified by detected number of backspliced reads can be found in Table S1. To validate that the higher numbers of detected circRNAs in capture sequencing libraries were not due to the differences in sequencing depth, some libraries were down-sampled to make sure that the sequencing depth of any capture library was no more than the depth of its matched Ribo-Zero library. B) In a VCaP cell line library, the overlap of detected circRNA species in exome capture and Ribo-Zero sequencing platforms was significant (Fisher exact test, P-value < 1×10⁻¹⁶). A threshold of two backspliced reads was applied. C) In VCaP, capture sequencing retained the relative abundance of circRNAs to their linear counterparts, comparable to Ribo-Zero library (Spearman Rank Correlation ρ = 0.65) (see STAR Methods). A threshold of two backspliced reads was applied. D) Backspliced events called from our MiOncoCirc pipeline were elevated (Mann-Whitney-Wilcoxon < 2.2×10⁻¹⁶) post RNase R treatment in VCaP and 22RV1 cell lines, further confirming that they were true circRNAs and not ligation artifacts. See also Table S1.

Figure 2.. Construction and overview of the MiOncoCirc compendium.

A) 868 high-depth, paired-end RNA-seq samples from previously published data sets as well as cell line panels and normal tissues were included. Additional details and abbreviations can be found in Table S1. B) Exome capture RNA-seq protocol and the bioinformatics pipeline for creation of MiOncoCirc. The unmapped reads from chimeric aligner (STAR) were annotated against the exon junctions. CIRCexplorer was used to call circRNA transcripts, and CODAC was used to annotate circRNAs involving two genes. FeatureCounts was used to quantify gene expression. C) MiOncoCirc is an online database that enables querying and downloading of circRNAs abundance across different tissues. Additional genomic data can be retrieved from previous studies (STAR Methods). See also Table S1.

Figure 3.. Features of circRNAs in MiOncoCirc and properties associated with expression.

A) The overlap of circRNA species in MiOncoCirc and CircBase was significant (Fisher Exact Test P-Value < 1×10⁻¹⁶). Only high confidence circRNAs which appeared in five or more samples were included in this comparison. B) Genes can form multiple circRNA transcripts. The number of circular transcripts increased proportionally with the number of exons per gene (binned to 10). C) Average expression of circRNA abundance (in normalized backspliced reads) vs. average expression of parental expression (in FPKM). Parent gene expression was grouped into bins of 50. Overall, there was no different in the mean of the bins (ANOVA P-Value = 0.12). This result agreed with Figure S2A in that the correlation was weak (Spearman’s ρ = 0.12). D) The distribution of Spearman’s rank correlation between circRNA abundance and their cognate parent gene expression, across all samples (gray), in prostate adenocarcinoma (PRAD, blue), and breast cancer (BRCA, red). Overall, the correlations were low (all medians < 0.28). E) Circular RNA abundance (in normalized backspliced reads) vs. sample fraction (%). There was a small portion of circRNAs (<2% of all circRNAs, generated from 589 genes, marked as “high”) that were detected in more than 90% of all samples. They also had higher expression compared to the median of all circRNAs (marked as “high” in the density plot of Figure S2A). F) These “high” circRNAs were flanked by significantly longer introns (Mann-Whitney-Wilcoxon P < 2.2×10⁻¹⁶) than the remaining 98% of circRNAs. Genome-wide introns were included as the control. See also Table S2 and Figures S1–S2.

Figure 4.. Identification of novel read-through (rt) circRNA species.

A) Schematic showing that genomic tandem duplications and circRNAs involving two genes can appear similar in paired-end RNA-seq. Specifically, when mates of a paired-end read were aligned in divergent orientation to exons of two adjacent genes, the result could be interpreted as either a duplication of a group of exons from two genes (Scenario 1), or a circularization from the downstream gene back to the upstream gene (Scenario 2). B) Schematic depicting the circular read-through event that can be generated from two adjacent genes, and their genomic features. C) The introns flanking rt-circRNAs were longer than genome-wide introns (Mann-Whitney-Wilcoxon P < 2.2×10⁻¹⁶). D) The introns flanking rt-circRNAs harbored more repetitive elements than genome-wide introns (Mann-Whitney-Wilcoxon P < 3×10⁻⁹). E) The frequency and distribution of the top 30 most abundant backspliced events involved neighboring genes in our compendium. F) The circular read-through event involving exon 3 of TTTY15 and exon 3 of USP9Y was chosen for validation in LNCaP cells. Post RNase R treatment, only the RT-qPCR product of outward facing primers was resistant to exoribonuclease degradation (see STAR Methods). See also Tables S3 and S7 and Figures S3–S5.

Figure 5.. Expression patterns and characteristics of circRNAs in cancer.

A) Tissue-specific heatmap of genes that can generate circRNAs, as demonstrated in 17 cancer cohorts from the MiOncoCirc compendium. A gene was considered to be consistently detected if it generated at least one high-confidence circRNA in more than 30% of samples of any given lineage (see STAR Methods). B) Volcano plot of circRNA abundances in 25 matched pairs of normal/localized prostate adenocarcinoma. Horizontal dash-line corresponded to FDR = 0.05. Vertical dash-line corresponded to fold-change > 1.5× (up-regulation) and fold-change < −1.5× (down-regulation). C) The correlation of log fold-change (FC) of circular RNA vs. log FC of linear expression. Again, circRNA abundances were downregulated overall (mean circular logFC = −0.9). We further stratified genes into groups based on the relationship between the linear and circular fold change. Group 1 (red) were circRNAs that were upregulated in cancer because their parent genes were upregulated. Group 2 (purple) were those circRNAs that were downregulated in cancer because their parent genes were also downregulated. However, there was a subset of circRNAs (Group 3, blue) downregulated in cancer with no corresponding change in parent gene expression. D) Total circRNA correlated with prostate cancer mRNA (“m”) proliferation markers calculated in FPKM (MCM10, TOP2A, MKI67, PCNA, KIAA0101, and NUSAP1). The size and the color scale of the dots indicate the values of pair-wise Spearman Rank Correlation. GAPDH mRNA expression was included as a negative control. Circular FBXO7, a highly abundant circRNA also showed remarkable negative correlation with proliferation index, even though its parental gene expression did not correlate. See also Tables S3–S4 and Figures S6–S7.

Figure 6.. Differential circRNAs in neuroendocrine prostate.

A) NEPC cases from our cohort, classified by cell morphology from pathology assessments, were all characterized by the upregulation of neuroendocrine markers (SYP, CHGA, CHGB, NCAM1, ENO2, ASCL1, MYCN, and AURKA) and downregulation of genes in the AR signaling pathway (AR, AMACR, KLK3, KLK2, FKBP5, and PSCA) compared to castration-resistant prostate cancer (CRPC) cases. B) The heatmap of 34 upregulated and 48 downregulated circRNAs with statistical significance (P < 0.01) in NEPC compared to CRPC cases. C) Comparing NEPC to CRPC, the most significantly upregulated circRNA was circ-AURKA (Mann-Whitney U test P = 3.17×10⁻⁹); the most significantly downregulated circRNA was circ-AMACR (Mann-Whitney U test P = 0.002). D) RT-qPCR of outward-facing primers of AURKA (backspliced from exon 6 to exon 3) in RNase R treated NCI-H660, a NEPC cell line, confirmed the circular structure of this molecule. P < 0.0001 calculated from one-way ANOVA. E) RT-qPCR of circular and linear AURKA in prostate cancer cell lines. Both circ-AURKA and linear-AURKA were expressed higher in NCI-H660 than in two non-NEPC cell lines, LNCaP and VCaP. See also Tables S5 and S7.

Figure 7.. Circular RNAs are more stable than cognate linear transcripts and can be detected in urine samples from prostate cancer patients.

A) Compared to their linear counterparts, circRNAs were resistant to RNase R degradation. Linear transcripts were detected by inward-facing RT-qPCR primers, while circular transcripts were detected by outward-facing RT-qPCR primers (**P < 0.0001, calculated from Student’s t test). B) After transcription inhibition by actinomycin D in LNCaP cells, linear transcripts (Linear) degraded faster than their corresponding circular transcripts (Circular). Samples were harvested at 0, 2, 4, 8, and 24 hours post-treatment. GAPDH was used as the control. The fold changes were calculated relative to the starting time point. circHIPK2 was selected to represent “high” class circRNAs. circLUZP2 represented “low” class circRNAs but with elevated expression in prostate cancer compared to normal. C) After incubating VCaP RNAs in plasma, the circular-to-linear ratio of circRNAs increased over time. Samples were harvested at 0, 15, 30, 45, 60, and 75 minutes. D) Circular RNAs were detected by exome capture RNA-seq of 3 urine samples from prostate cancer patients. These circRNAs greatly overlapped with circRNAs identified in prostate cancer tissues from the MiOncoCirc cohorts. See also Tables S6–S7 and Figure S7E.