
==== Front
RNA Biol
RNA Biol
RNA Biology
1547-6286
1555-8584
Taylor & Francis

39257052
10.1080/15476286.2024.2395718
2395718
Version of Record
Research Article
Research Paper
A systematic analysis of circRNAs in subnuclear compartments
A. BREZSKI ET AL.
RNA BIOLOGY
Brezski Andre a
Murtagh Justin b
Schulz Marcel H. b c d
Zarnack Kathi a *
a Buchmann Institute for Molecular Life Sciences (BMLS) & Institute of Molecular Biosciences, Goethe University Frankfurt , Frankfurt am Main, Hesse, Germany
b Department of Medicine, Institute for Computational Genomic Medicine and Institute of Cardiovascular Regeneration, Goethe University Frankfurt , Frankfurt am Main, Hesse, Germany
c Cardio-Pulmonary Institute, Goethe University Frankfurt , Frankfurt am Main, Hesse, Germany
d German Center for Cardiovascular Research, Partner site Rhein-Main , Frankfurt am Main, Hesse, Germany
CONTACT Kathi Zarnack kathi.zarnack@uni-wuerzburg.de Buchmann Institute for Molecular Life Sciences (BMLS) & Institute of Molecular Biosciences, Goethe University Frankfurt, Frankfurt am Main, Hesse 60438, Germany
* Current affiliation for Kathi Zarnack is Department of Bioinformatics, Theodor Boveri Institute, Julius Maximilians University Würzburg, Würzburg, Germany.

10 9 2024
2024
10 9 2024
21 1 116
Integra06 9 2024
Integra06 9 2024
21 6 2024
13 8 2024
© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.
2024
The Author(s)
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

ABSTRACT

CircRNAs are an important class of RNAs with diverse cellular functions in human physiology and disease. A thorough knowledge of circRNAs including their biogenesis and subcellular distribution is important to understand their roles in a wide variety of processes. However, the analysis of circRNAs from total RNA sequencing data remains challenging. Therefore, we developed Calcifer, a versatile workflow for circRNA annotation. Using Calcifer, we analysed APEX-Seq data to compare circRNA occurrence between whole cells, nucleus and subnuclear compartments. We generally find that circRNAs show higher abundance in whole cells compared to nuclear samples, consistent with their accumulation in the cytoplasm. The notable exception is the single-exon circRNA circCANX(9), which is unexpectedly enriched in the nucleus. In addition, we observe that circFIRRE prevails over the linear lncRNA FIRRE in both the cytoplasm and the nucleus. Zooming in on the subnuclear compartments, we show that circRNAs are strongly depleted from nuclear speckles, indicating that excess splicing factors in this compartment counteract back-splicing. Our results thereby provide valuable insights into the subnuclear distribution of circRNAs. Regarding circRNA function, we surprisingly find that the majority of all detected circRNAs possess complete open reading frames with potential for cap-independent translation. Overall, we show that Calcifer is an easy-to-use, versatile and sustainable workflow for the annotation of circRNAs which expands the repertoire of circRNA tools and allows to gain new insights into circRNA distribution and function.

KEYWORDS

CircRNA
subnuclear compartments
workflow
nuclear speckles
circRNA translation
German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) 403584255 This work was supported by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) via TRR 267 [Project ID 403584255] Project A01 to K.Z., and Project Z03 to M.H.S.
==== Body
pmcIntroduction

CircRNAs are a type of RNA which is characterized by a circular structure. The functions of circRNAs are multifaceted, ranging from important roles in human development to the specific regulation of parental genes and also disease markers [1,2]. The expression of circRNAs is often specific to a certain cell type and developmental stage [3,4]. Some circRNAs are highly expressed in human diseases including cancer and cardiovascular diseases such as hypertension, making them promising candidates as biomarkers or even therapeutic targets [5,6].

CircRNAs are generated by an unconventional pre-mRNA splicing process called back-splicing, in which a 5’ splice site is linked to an upstream 3’ splice site [7–11]. The biogenesis of circRNAs partially competes with linear splicing (i.e. the linkage to a downstream 3’ splice site) and can be modulated by RNA-binding proteins (RBPs) or complementary sequences in the flanking introns, among others [12–14]. Analogous to linear alternative splicing, one or more circRNAs can be generated from a single parental pre-mRNA by using different 3’ and/or 5’ splice sites [15]. Although circRNAs mostly consist of the same exons as their linear counterparts, they usually have much longer half-lives. This is mainly due to their resistance to exonuclease-mediated decay, which allows circRNAs to accumulate over time [10,15,16]. However, their slow biogenesis rates may also lead to a dilution of circRNAs in fast proliferating cells [17,18].

Within the cells, many circRNAs are exported to the cytoplasm [19], where they have been described to exert various functions including the sponging of microRNAs (miRNAs) and RBPs [20–23]. Moreover, recent studies showed that translation of circRNAs is also possible [24–26]. Nonetheless, a certain fraction of circRNAs remains in the nucleus. Here, they can, for instance, recruit proteins to chromatin or increase the retention of proteins in the nucleus [26–28]. Despite these nuclear functions, the spatial distribution of circRNAs within the nucleus remains largely unknown. Of note, an emerging experimental protocol termed APEX-Seq (named after the utilized peroxidase enzyme APEX2) allows to study the spatial distribution of RNAs using high-throughput sequencing. In this proximity labelling-based method, all RNAs in the vicinity of an APEX2-tagged marker protein specific for a given compartment are biotinylated and then extracted [29]. Yet, such data has not been used to assess the spatial distribution of circRNAs.

Computationally, circRNAs can be detected in high-throughput sequencing data of total RNA that includes all RNA species irrespective of the presence of a polyA tail (referred to as total RNA-seq). The distinction between linear transcripts and circRNAs is exclusively possible via sequencing reads that span across the back-splice junction (BSJ), indicated by a non-linear alignment of the two read fragments to the genome (referred to as chimeric alignment). Thus, only a very limited number of sequencing reads is available for circRNA analyses. Accordingly, the detection and quantification of circRNAs is challenging, especially for datasets with a short read length or a low sequencing depth. While most tools use these chimerically mapped BSJ reads, some of the newer tools also work with machine learning [30,31]. Although there are many different tools for the detection and functional characterization of circRNAs, each of them usually covers only a part of the analyses [15,32,33].

Here, we present Calcifer, an efficient workflow for the automated analysis of circRNAs in total RNA-seq data. It effectively combines circRNA detection by established tools with an in-depth analysis of their characteristics and possible functions, including miRNA binding, RBP binding, and possibly translation. In addition, Calcifer provides a quantitative assessment of circRNA abundance, also in relation to the corresponding linear transcripts, which allows for comparisons between multiple conditions. The integrity of the results is ensured by several quality control and filtering steps that check for the properties of the detected circRNAs. At the end of the workflow, the entire results are summarized in a tab-separated text file, which can be further analysed. To showcase the performance of Calcifer, we analysed published APEX-Seq data to investigate the distribution of circRNAs in different subnuclear compartments [34]. We found that circRNAs are depleted from nuclear speckles, indicating that the excess of splicing factors in this compartment may promote canonical splicing and thereby prevent back-splicing.

Results

Calcifer – an automated workflow to detect and analyze circRNAs

To allow for an automated analysis of circRNAs, we developed Calcifer, a workflow that efficiently identifies circRNAs from datasets with one or multiple conditions and offers comprehensive downstream analyses related to circRNA biogenesis and function (Figure 1). For the initial circRNA detection, we utilize two distinct tools based on different genomic alignment algorithms: CIRCexplorer2 [15] uses splice-aware alignments from STAR [35] which directly include chimeric alignments, whereas CIRI2 [32] takes as input the continuous alignments from BWA [36] and follows a seed matching approach. These tools were selected because they showed good results in direct comparison with other circRNA detection tools [37]. In addition, both tools take different approaches for circRNA detection, which is beneficial to the yield. We chose this dual approach since using a single tool generally captures fewer circRNAs compared to employing multiple tools [31,38]. Consistently, the combination of tools increases the total number of detected circRNAs, underlining that the two tools complement each other (Figure S1A). The resulting circRNAs from both tools are then merged and filtered for a high-confidence circRNA set as shown before [39] (Figure 1A): We perform a secondary search for BSJ reads of all detected circRNAs based on the chimeric alignments of STAR and thus double-check the findings of the two tools. We re-count the supporting reads and use the chimeric alignment as ground truth. The remaining circRNAs are then filtered further. Therefore, only circRNAs which are detected by at least one tool and have sufficient BSJ support in the chimeric alignments are considered for further analysis. In brief, we expect circRNAs to possess canonical splice site sequences, no overly long distances between splice sites (<100 kilobases), and no encompassing junctions (i.e. BSJ between the read mates for paired-end reads). We also filter for support by ≥2 unique BSJ reads, which reduces the number of false positives from only one isolated BSJ read. If available, the 2 required unique BSJ reads can also originate from different replicates and do not have to appear in the same sample. For a more stringent filtering, the minimum number of unique BSJ supporting reads can be freely set by the user. Uniqueness of alignments is identified using the CIGAR strings of the chimeric alignments. Altogether, the combination of tools and filtering steps in Calcifer allows to obtain a high-confidence set of circRNAs for the major part of the workflow, the circRNA annotation. Figure 1. The Calcifer workflow. (A) Overview of the Calcifer workflow to detect, quantify and analyse circRNAs. The input consists of RNA-seq data as well as genome indexes, a reference and databases for the downstream analysis. CircRNA detection is based on the output from CIRCexplorer2 and CIRI2, which is further processed and filtered. The identified circRNAs are then subjected to downstream analyses, including the prediction of miRNA target sites, RBP binding sites as well as putative open reading frames (ORFs) and resulting peptides. The predicted peptides are further filtered for m6A RNA modification sites or internal ribosomal entry site (IRES)-like sequences upstream of the respective ORF and (sub)peptides that are unique within the whole proteome. The output of the workflow contains results for each analysis step, ordered in folders, and a general output file with an overview of all results for each circRNA. (B) Detailed overview of different downstream analyses of the Calcifer workflow. Circular-to-linear ratio is computed based on back-splice junction (BSJ) read support and respective linear junction read support. The sequence-based downstream analysis includes miRNA target site prediction, RBP binding site prediction and an extensive ORF analysis combined with detection of putative unique (sub-)peptides. (C) Cumulative runtime of the different steps of the workflow. Runtimes are based on a 2 × 100 nt paired-end dataset with 4 million reads. (D) Runtimes for the workflow depending on the number of replicates. Each replicate consists of 4 million 2 × 100 nt paired-end reads. (E) overlap of all circRNAs detected in the whole-cell samples (67.78%), all circRNAs detected in the nuclear samples (71.06%) and circRNA detected exclusively in the nuclear samples (30.64%) with the circRNA database circBase [40,41].

Once the circRNAs are identified, Calcifer performs several downstream analyses (Figure 1B): Among these, it calculates a relative ‘percent circularized’ value from the BSJ reads versus all linearly spliced junction reads that originate from the same splice sites. This value can be considered when comparing different conditions where inverse changes between the linear and circular isoform expression might occur. Calcifer also extracts the exonic sequence for each circRNA, considering all intervening exons between the back-splice sites, which is needed for the subsequent sequence analyses. However, it should be mentioned that the exact sequence of the circRNAs cannot be determined by short-read sequencing, as the only evidence for the circRNAs is given by the reads that lie across the BSJ. Regarding the biogenesis of circRNAs, Calcifer additionally checks for RBP binding sites in a fixed 250 nucleotide (nt) window around the back-splice sites of each circRNA. Related to possible functions of the circRNAs, the workflow provides in silico analyses of RBP and miRNA binding on the circRNA, as well as open reading frame (ORF) predictions including analyses of the possible resulting peptides (Table 1). For the RBP, miRNA and ORF analysis, a pseudo-circular sequence is also generated for all circRNAs: an overhang of 25 nt over the BSJ is appended to the respective opposite end, thus simulating the sequence over the BSJ. In order to consider ORFs that are located several times around the circRNA, a multi-cycle sequence is also defined, which consists of four repetitions of the linear sequence. Furthermore, a count matrix is created for all linear and BSJ reads that can be used for differential gene expression analysis.Table 1. Calcifer output.

Column	Description	Example	
CircID	Back-splice junction (BSJ) position as circID	chrX:131749305 -131,794,466:-	
Parental_gene	Name of the parental gene	FIRRE	
Type	CircRNA type as detected by CE2/CIRI2	exon	
Unique_bsr	Sum of unique back-splice junction reads	80	
Con_1_clr	Circular-to-linear ration (CLR) of circRNA per condition	0.74744	
Mirna_binding_site_density	MiRNA binding sites per nt	0.09624	
Most_mirna	Most single miRNA binding site counts and proportion	hsa-miR-877-5p:5:0.03049	
Rbp_circ_binding	RBP binding sites on the circRNA	PABPC4:1	
Rbp_bsj_binding	RBP binding sites around the BSJ	HuR:1	
Linear_seq_orf	Unique ORFs on the linear sequence and maximum length	12:61	
Pseudo_circular_seq_orf	Unique ORFs on the circular sequence and maximum length	0:0	
Multi_cycle_seq_orf	Unique ORFs on the multi-cycle sequence and maximum length	0:0	
Unique_peptide_regions	Partially or complete unique putative circRNA peptides	ORF1:0-21:ireslike_0:m6a_1	
Description and example for each column in the summary output of Calcifer. The output is tab-separated and contains information about each circRNA and all downstream analyses. Further results can be found in the respective output folder.

In order to test the runtime of the entire workflow, we analysed a dataset with 4,000,000 paired-end reads with a combined length of 200 nt (Figure 1C). For this analysis, we utilized 10 CPUs and 157.35 GB RAM. The runtime for the complete analysis was 45 minutes, with 27 minutes for the circRNA detection and 18 minutes for the downstream analysis. It can be seen that STAR, BWA, CIRI2, the count matrix generation and the miRNA target site prediction had most impact on the runtime (Figure 1C). We also repeated this analysis with 2–5 replicates, all containing the same number of raw paired-end reads. Here, it can be seen that the runtime increases at most linearly per replicate (Figure 1D).

In addition to a fully automated run, the modules of Calcifer can also be executed individually, such as the analysis with CIRCexplorer2 and CIRI2, the merging of results as well as all downstream analyses. This can be useful if partial results are already available. There is also an alternative mode to perform the downstream analyses on any list of circRNAs, making it possible to combine Calcifer with other circRNA detection tools or even apply it to complete circRNA databases. In summary, Calcifer is a versatile tool for circRNA detection and analysis. This is made possible by its modular structure and the various entry points into the workflow that can run on raw sequencing data as well as lists of already detected circRNAs.

Calcifer can be used to efficiently process a large-scale RNA-seq dataset

To test the performance of Calcifer and its different modes on large-scale data, we processed an APEX-Seq dataset which provides transcriptome data for eight different subnuclear compartments (Figure S1B) [34]. The complete dataset consists of 14 conditions with 2–4 replicates per condition, including whole-cell samples with total RNA (i.e. the pre-pull-down [pre-PD] samples) as well as nuclear samples enriched for RNAs from eight subnuclear compartments (i.e. the pull-down [PD] samples), cumulating in 96 fastq files with a total of 2.3 billion reads (75.8 GB). The mean sample size was 84.0 million reads for the pre-PD and 80.4 million reads for the PD samples (Figure S1C). In the genomic alignment step, STAR retrieved a mean of 365,202 chimeric reads for the pre-PD and 346,103 for the PD samples (Figure S1D). Processing the complete dataset with Calcifer took approximately 10 days on a Linux system with 10 CPUs and 157.35 GB RAM.

Using the union of all samples, CIRCexplorer2 and CIRI2 detected 15,773 and 9,754 different circRNAs, respectively. After merging and applying the general filters, this yielded an initial set of 17,168 circRNAs. By further filtering for ≥2 unique BSJ reads, Calcifer identified 7,208 high-confidence circRNAs in the complete APEX-Seq dataset, which were used for further analysis. An overlap of the detected circRNAs with the circRNA database circBase [40,41] showed that 67.78% of all circRNAs detected in the whole-cell samples as well as 71.06% of all circRNAs in the nuclear samples had been reported previously, supporting the validity of our results (Figure 1E, left and middle). The overlap was considerably smaller for circRNAs that were exclusively detected in the nuclear samples (30.64%, Figure 1E, right), indicating that the targeted APEX-Seq approach increases the sensitivity for circRNA detection in the nucleus.

Automated downstream analyses allow for a detailed characterization of the detected circRNAs

After detecting the circRNAs from the raw sequencing data, we performed the downstream analyses implemented in Calcifer. These include circRNA quantification, miRNA target site detection and ORF prediction. We used the Calcifer output to globally characterize the high-confidence circRNAs from the complete APEX-Seq dataset. Overall, the 7,208 circRNAs originated from 3,633 genes, with mostly one circRNA isoform per gene (n = 2,125), while a few genes hosted more than 6 isoforms (n = 135) (Figure 2A). From the different circRNA host genes, 3,443 were protein-coding genes, while a smaller fraction of circRNAs originated from 148 long noncoding RNA (lncRNA) genes and 23 pseudogenes (Figure 2B). The mean length of all detected circRNAs was 1,235 nt, while the median length was 681 nt (considering all intervening exons; Figure 2C). Figure 2. Characterization of the circRNA repertoire in the APEX-Seq data. (A) Amount of different circRNA isoforms per gene. Most genes (n = 2,126) have a single circRNA isoform. (B) Pie chart of the types of genes which host at least one circRNA. Protein-coding genes are the majority (n = 3,443). (C) Log10-scaled circRNA length distribution (mean = 1,235 nt; median = 681 nt). (D) Quantitative analysis of miRNA target sites on circRNAs. The ratio of the single most abundant miRNA target site to total binding sites is shown against the target site count for this specific miRNA. Marked are ciRS-7 (37 target sites for hsa-miR-7-5p) and circATP9B(5, 6, 7) (32 target sites for hsa-miR-6794-3p) as circRNAs with possible miRNA sponging function. The remaining circRNAs either have a low ratio of single miRNA target sites or low target site counts. (E) Pie chart of circRNAs with predicted RBP binding sites within the circRNA sequence, around the back-splice sites or on both sequences. There are 6,433 circRNAs which have at least one RBP binding site around the back-splice sites and on the circRNA sequence. 30 of the detected circRNAs possess no RBP binding sites. (F) Top 10 RBPs in regard of amount of circRNAs with RBP binding sites on the circRNA sequences. SRSF1 (n = 3,152) and HuR (n = 2,781) are the two RBPs with the most circRNAs possessing at least one binding site on the circRNA sequence. (G) Top 10 RBPs in regard of amount of circRNAs with RBP binding sites around the BSJ. HuR (n = 5,049) and TIA1 (n = 3,373) are the two RBPs with the most circRNAs possessing at least one binding site around the BSJ. (H) Top 10 RBPs in regard of amount of circRNAs with RBP binding sites around the BSJ on both sides. HuR (n = 1,958) is the RBP with the most circRNAs possessing a binding site on both sides of the BSJ. (I) Percentage of total circRNAs with different types of ORFs. 6,464 (89.68%) circRNAs possess at least one ORF on the linear sequence, whereas 3,447 (47.82%) and 3,951 (54.81%) circRNAs have at least one ORF on the circular and multi-cycle sequence, respectively. (J) Comparison of putative unique peptides originating from non-linear ORFs on circRNAs, which possess a DRACH motif and/or IRES-like sequence upstream of the start codon. In total, 2,036 circRNAs have at least one non-linear ORF with a putative unique peptide, which is not found elsewhere in the proteome.

As part of the downstream analysis, Calcifer uses the algorithm miRanda [42] to predict putative miRNA target sites in the circRNA sequences. As a result, it returns the number of target sites for the most abundant miRNA as well as its relative contribution to the total number of miRNA target sites predicted on this circRNA (referred to as density of miRNA target sites). We assume that a circRNA, which functions as a miRNA sponge, should have similar properties like the well-studied circRNA ciRS-7 (also known as CD1as), which harbours a high density of target sites for one prominent miRNA, in this case miR-7 (hsa-miR-7-5p) [22]. Consistently, ciRS-7 has been shown to regulate the activity of miR-7 via miRNA sponging [22]. In the present dataset, we observed that most circRNAs either possessed a generally high density of target sites for many different miRNAs or, if one miRNA target site was predominant, harboured only a small overall number of target sites, but rarely a combination of both features (Figure 2C). The notable exceptions were ciRS-7 and circATP9B(5, 6, 7) (naming of circRNA isoforms follows the nomenclature proposed by [43], listing all potentially included exons), which we newly detected as a putative miRNA sponge with 32 target sites for the miRNA hsa-miR-6794-3p (Figure 2D).

RBP binding in the vicinity of the back-splice sites can modulate circRNA biogenesis, while RBPs on the mature circRNAs may influence their downstream functions and export. To address both possibilities, Calcifer uses the motif search algorithm FIMO [44] to predict RBP binding sites in a 250 nt window around the BSJs into the intron and 25 nt into the circRNA sequence (window 1) as well as on the complete mature circRNA sequences (window 2) using known binding motifs for 85 RBPs [45]. We note that the window size around the back-splice sites (550 nt, window 1) is close to the median length of all circRNAs in the dataset (681 nt, window 2). For the circRNAs from the APEX-Seq dataset, Calcifer predicted a least one RBP binding site around the back-splice sites for 6,746 circRNAs (Figure 2E). In addition, we found 6,866 circRNAs with putative RBP binding sites within the circRNA sequence, including 6,433 circRNAs with hits in both regions. 30 circRNAs did not yield any results. Taking a closer look at the RBPs predicted to bind within the circRNA sequences, the top 3 RBPs occurring on most circRNAs were SRSF1, HuR (ELAVL1) and RBM5 (Figure 2F). For the sequences around the back-splice sites, HuR was also most often predicted to bind, followed by TIA1 and HNRNPC (Figure 2G). In terms of RBP binding on both sites of the BSJ, HuR, TIA1 and HNRNPC were the top 3 RBPs (Figure 2H). In total, we predicted 81 different RBPs to bind around the back-splice sites, 66 RBPs to bind on both sides of the back-splice sites, and 81 RBPs to bind on the circRNA sequences.

Several recent reports suggested cap-independent translation into polypeptides as a putative function of circRNAs [24,26]. Calcifer, therefore, screens all circRNAs for complete ORFs using a custom approach. Since ORFs may span across the BSJ and possibly run around the circRNA multiple times, Calcifer considers three sequence variants, including (1) the linear sequence consisting of all included exons, (2) a pseudo-circular sequence which additionally extends 25 nt across the BSJ at both ends, and (3) a multi-cycle sequence obtained from concatenating four times the linear sequence (Figure S1G, H). Remarkably, we found that 6,952 out of 7,208 circRNAs (96.45%) harboured at least one complete ORF (Figure 2I). These included 6,464 circRNAs with at least one complete ORF on the linear sequence, and 3,447 circRNAs with at least one ORF on the pseudo-circular sequence, i.e. starting and/or terminating within 25 nt across the BSJ. Moreover, 3,951 circRNAs showed at least one ORF on the multi-cycle sequence, i.e. running at least once around the circRNA. Multiple occurrences of an ORF were counted only for the first appearance to prevent redundant duplicated ORF detections, especially on the multi-cycle sequence. Therefore, all ORFs determined on the circular or multi-cycle sequence are unique and only segments of these ORFs might also appear on the linear sequence (Figure S1E).

Previous reports suggested different scenarios to facilitate the cap-independent translation of circRNAs, involving either N6-methyladenosine (m6A) RNA modifications or short internal ribosomal entry site (IRES)-like sequences to initiate translation [24,26]. Thus, Calcifer includes a sequence search for the consensus motif DRACH for m6A sites (with D = A, G or T, R = A or G, and H = A, C or T) as well as for IRES-like motifs (Table S1) taken from Fan et al. [24]. Both motif types are searched for in the 10 nt upstream of the start codon of each ORF. Finally, to facilitate experimental validation, e.g. by using mass spectrometry, the putative ORFs are in silico translated into the encoded peptide sequences. These are further evaluated for unique sequence stretches that specifically arise from the circRNA ORF but not from the canonical ORF in the cognate linear transcript. The identification for these unique circRNA peptides is achieved via comparing all putative circRNA peptides against the whole proteome and scanning for previously unknown peptide sequences and subsequences of at least 10 amino acids in length. For all circRNAs with a non-linear ORF and an upstream IRES-like or DRACH motif, we detected 9,315 unique circRNA peptides originating from translation of 2,036 different circRNAs (Figure 2J).

To assess whether the predicted ORFs on circRNAs are translatable, we compared our data with the results of CircCode, a Ribo-seq analysis of circRNAs with ORFs detected in MCF7 cells [46]. In total, 149 circRNAs detected in the APEX-seq data showed evidence of translation in the CircCode analysis (Figure S1F). Of these, 144 possessed at least one predicted ORF, and 122 also harboured at least one m6A or IRES-like site upstream of the ORF and the putative translation product contained unique (sub-)peptides. The generally low number of overlapping circRNAs may be partly due to the insufficient depth of the Ribo-seq data, the low abundance of circRNAs in cells and the fact that different cell lines were used.

Overall, these results highlight the versatility of Calcifer to detect, quantify and characterize circRNAs. The reported metrics provide useful insights into their putative interaction partners and shed light on potential functions, offering starting points into further in-depth analyses.

Most circRNAs are more abundant in the cytoplasm

The APEX-Seq dataset consists of a series of whole-cell samples with total RNA (i.e. including cytoplasmic RNAs) as well as nuclear samples in which RNAs associated with different subnuclear compartments have been selectively enriched [34]. In a first step, we compared the whole-cell versus nuclear samples to investigate the distribution of circRNAs between the nucleus and the cytoplasm. We found that the number of BSJ reads was significantly reduced in the nuclear samples, with on average 2,017 normalized BSJ reads per sample in the whole-cell samples compared to only 398 in the nuclear samples (Figure 3A). When looking at the overlap of circRNAs between the sample types, the vast majority of nuclear circRNAs were also detected in the whole-cell samples (722 out of 895, 80.7%). Overall, the whole-cell samples harboured a total of 7,035 circRNAs, being over 7-fold more than in the nucleus. Only 173 circRNAs were exclusively detected in the nucleus (Figure 3B). To visualize the distribution of circRNAs between the nucleus and the cytoplasm, we computed z-scores of normalized BSJ read counts across all samples. For comparison, we looked at total read counts for linear transcripts, which split into distinct subsets of predominantly nuclear or cytoplasmic transcripts (Figure 3C). In contrast, the z-score distribution of BSJ read counts supported a universal depletion of circRNAs from the nucleus (Figure 3D). Together, these observations indicated that circRNAs are mostly exported to the cytoplasm and hence underrepresented in the nucleus. Figure 3. CircRNA comparison between whole cell and subnuclear compartments. (A) Boxplot of the normalized sum of BSJ read counts per sample divided by whole cell and the subnuclear compartments. There is a significant reduction in BSJ reads between whole cell and subnuclear compartments. (B) Overlap of detected circRNAs from whole cell and nucleus. 173 circRNAs were only detected in the nucleus. (C) Heatmap of z-scores to compare read counts for all expressed genes between whole cell and subnuclear compartment samples. (D) Heatmap of z-scores to compare BSJ reads for each circRNA between whole cell and subnuclear compartment samples. Most circRNAs show lower abundance in the subnuclear compartments. (E) Comparison of proportion of circRNAs containing at least one ORF on the linear, circular or multi-cycle sequence, as well as a putative resulting unique peptide. CircRNAs unique to the nucleus show a significantly higher proportion of a possible ORFs and putative unique peptides compared to circRNAs unique to the whole cell.

Following up on the different abundance of circRNAs in nucleus and cytoplasm, we compared the properties of circRNAs that were exclusively detected in the whole-cell (n = 6,313) or nuclear samples (n = 173). For the RBP binding sites, we found a significantly higher amount of binding sites on circRNAs unique to the whole-cell samples (Figure S2A). The comparison of the top 5 RBPs across all circRNAs showed no significant differences for the circRNA sequences (Figure S2B). However, the amount of circRNAs with binding sites for HNRNPC and HNRNPCL1 around the BSJs was significantly lower for nuclear circRNAs (Figure S2C). For miRNA binding sites, we found no differences between unique circRNAs from nuclear and whole-cell samples (Figure S2D). Interestingly, however, we found that all types of predicted ORFs as well as unique peptides were significantly decreased for the nuclear circRNAs (Figure 3E). The increased translatability of cytoplasmic circRNAs may be taken as evidence that circRNAs in the cytoplasm are more often translated than previously anticipated. Consistently, the nuclear circRNAs are depleted for this functionality.

CircRNAs are globally depleted from nuclear speckles

Although the circRNAs were generally more abundant in the cytoplasm, a substantial number was still expressed in the nucleus. We, therefore, proceeded to analyse the distribution of circRNAs within the nucleus. The APEX-Seq dataset differentiates between eight different subnuclear compartments, namely Cajal bodies, histone locus bodies, nuclear speckles, the nucleolus, PML bodies, SAM68 bodies, and the nuclear periphery, with each compartment between represented by one or more marker proteins (Figure S1B) [34]. To investigate to which extent a given RNA is present in its circularized form, we compared the percent circularized values for all circRNAs in the different compartments (Figure 4A). As previously observed [39], most circRNAs had low values, indicating that the circRNAs contributed only a minor fraction compared to the associated linear transcripts. The minimum average percent circularized value was 10.4% for circRNAs in SAM68 bodies, while the maximum average percent circularized value was 21.7% for circRNAs in the nuclear periphery. Surprisingly, we found that circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) (NE = novel exon) and ciRS-7 were not only present in most compartments, but also showed a very high percent circularized value. Comparing the sum of normalized BSJ reads among subnuclear compartments revealed similar values (Figure 4B). Interestingly, the exception were nuclear speckles which showed a significantly lower amount of BSJ reads. This observation fits well with the function of nuclear speckles as splicing factories, which store excess splicing factors, such as SR proteins, which were previously reported to promote linear splicing fidelity and suppress back-splicing [47]. To directly test this assumption, we looked at the behaviour of circRNA expression upon SRSF6 overexpression in HeLa cells [48] (Figure 4C). Indeed, we found that circRNAs were exclusively down-regulated upon SRSF6 overexpression, whereas linear transcripts responded equally in both directions. This observation supported the notion that an excess of SR proteins, as found in nuclear speckles, generally prevents back-splicing. Figure 4. Comparison of circRNA properties between subnuclear compartments. (A) Violin plots of percent circularized values for circRNAs of each subnuclear compartment. The percent circularized shows a slight increase for the histone locus body, the nuclear periphery and the nucleolus. CiRS-7, a circRNA originating from CDR1AS, is nearly completely circularized and is only absent in the histone locus body. CircFIRRE has an average percent circularized value of 46.6% and was detected in all subnuclear compartments. (B) Boxplots of the normalized BSJ read counts for each subnuclear compartment. The distribution of the normalized BSJ reads is relatively even between all subnuclear compartments, besides nuclear speckles which show a significant reduction compared to control. This plot shows the same data as in Figure 3A, with distinct highlighting. (C) Violin plot of log2-transformed fold changes of significantly regulated circRNAs (n = 1,826) and linear up- (n = 8,309) and down-regulated (n = 12,931) RNAs upon SRSF6 overexpression in HeLa cells (adjusted P value < 0.05). (D) Heatmap of z-scores to compare circRNA read counts between different subnuclear compartments. Nuclear speckles show a general depletion of circRNAs. (E) Heatmap of z-scores to compare circRNA read counts between nucleoplasm control (NLS pull-down) and SAM68 body and nuclear speckles, respectively. For the SAM68 body, there is an enrichment of multiple circRNAs like circCANX(9). We note that the enrichment for ciRS-7 is due to a batch effect and not significant. In comparison, there are only depleted circRNAs for nuclear speckles compared to control.

Taking a closer look at the circRNAs in the different compartments, we found that only eight circRNAs were shared between all compartments, while most detected circRNAs were unique for just one compartment (Figure S3A). As before, we used z-scores of BSJ read counts to compare circRNA expression between the different subnuclear compartments (Figure 4D, E). This confirmed that the nuclear speckles showed a general depletion of circRNAs, even for circRNAs that were generally found across compartments like circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8). In contrast, most other compartments like SAM68 bodies or the nucleolus showed an enrichment in a few circRNAs.

The distinct nuclear distribution of circRNAs was exemplified by circCANX(9) (Figure 5A, B). CircCANX has been previously described in the context of osteosarcoma and as a miRNA sponge [49], although our analyses did not support the latter. In total, we detected eight different isoforms of circCANX (Figure S3C), with the major isoform consisting of only exon 9 (Figure 5A). In contrast to most other circRNAs, circCANX(9) was more frequently circularized in the nucleus than in the cytoplasm (Figure 5C). Within the nucleus, circCANX(9) showed an enrichment in the SAM68, PML and Cajal body samples compared to the nucleoplasm (i.e. the NLS pull-down control; Figure 5D). Figure 5. CircRNA analyses for circCANX and circFIRRE. (A) IGV visualization of the major isoform of circCANX, spanning over one exon. (B) Boxplot of BSJ read counts for circCANX(9) between whole cell and nucleus. There is a non-significant increase in comparison of the nucleus (mean = 20.71 reads) to the whole cell (mean = 15.50 reads). (C) Boxplot of percent circularized values for circCANX(9) between the whole cell (mean = 0.40%) and the nucleus (mean = 0.85%). Circularization is significantly increased in the nucleus compared to the whole cell. (D) Boxplot of BSJ read counts for circCANX(9) between all subnuclear compartments. The SAM68 body shows the highest amount of normalized BSJ reads (mean = 16.12 reads), whereas the nuclear periphery has the lowest amount (mean = 1.99). (E) Boxplot of BSJ read counts for circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) between whole cell and nucleus. There is a significant decrease in the nucleus (mean = 22.36 BSJ reads) compared to the whole cell (mean = 70.57 BSJ reads). (F) Boxplot of normalized BSJ read counts for circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) between all subnuclear compartments. The nucleolus shows the highest amount of normalized BSJ reads (mean = 15.97 reads), whereas nuclear speckles have the lowest amount (mean = 3.95 reads). (G) Boxplot of percent circularized values for circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) between the whole cell (mean = 63%) and nucleus (mean = 46%). Circularization is significantly decreased in the nucleus compared to the whole cell. (H) IGV visualization of the major isoform of circFIRRE and further selected isoforms. The majority of the read coverage is seen in the range of the major circRNA isoform.

circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) is the predominant isoform of the lncRNA FIRRE

CircFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) is a circRNA that arises from the lncRNA FIRRE [50]. We found that this circRNA had a high BSJ read support in all whole-cell and nuclear samples (Figure 5E). Consistently, circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) was detected in nearly all subnuclear compartments (Figure 5F). Interestingly, it showed high percent circularized values across all samples, reaching up to 59.6% in the nuclear periphery (LMNA pull-down), indicating that the majority of FIRRE transcripts are present in the circularized form (Figure 5G). This was further supported by looking at total read coverage which was substantially elevated in exons 3–8 forming the circRNA. In total, we detected 13 different isoforms of circFIRRE, with one major isoform across all sample types (Figure S3D). We were able to determine the composition of the major circFIRRE isoform due to its high circularization and low read coverage outside the circRNA segments. The major isoform circFIRRE (3, NE, NE, 4, NE, 5, 6, 7, 8) contained several annotated exons (FIRRE exons 3–8) as well as 3 non-annotated exons (Figure 5H).

Thus, we conclude that circRNAs are generally more abundant in the cytoplasm but can be enriched in certain subnuclear compartments. Intriguingly, we observe a strong depletion of circRNAs from nuclear speckles, indicating that this compartment disfavours back-splicing or circRNA accumulation in this compartment.

Discussion

Calcifer – an automated workflow to detect and analyze circRNAs

Calcifer is an easy to instal and versatile workflow for an in-depth look into circRNAs from all kinds of datasets. The different run modes allow the workflow to be applied in a targeted way and to be utilized for already existing circRNA sets. The focus of Calcifer is to link a strongly filtered high-confidence set of circRNAs with a broad downstream analysis. Although the amount of circRNAs in the fully processed Calcifer results was lower than in the initial output of CIRCexplorer2 and CIRI2, the final circRNAs are expected to be more reliable due to the usage of more than one detection tool and a stringent filtering. Calcifer combines the detection of circRNAs with a detailed downstream analysis and annotation. Therefore, this enables an automatized approach without manual usage of additional tools. With the analysis of the APEX-Seq dataset, we demonstrated the functionality of Calcifer which allows evaluation of diverse datasets and comparison of high-confidence circRNAs between different conditions. In addition, there is a broad downstream analysis for each circRNA, which includes a comprehensive prediction of ORFs and the potentially resulting peptides. Together, these features make Calcifer a versatile tool for the combined detection and analysis of circRNAs.

Most circRNAs lack evidence for miRNA sponging but contain putative ORFs for translation

Using the detailed downstream analysis results of Calcifer, we investigated possible functions of the detected circRNAs. One of the most frequently reported purported roles for circRNAs is miRNA sponging, i.e. restricting the availability of a certain miRNA by sequestration to recurrent miRNA target sites in the circRNA. However, our results support the notion that this function is in fact limited to very few circRNAs including CDR1as (also known as CiRS-7) [22]. In the APEX-Seq dataset, we detected only a few circRNAs, which harboured a high number of recurrent target sites for a single miRNA, in line with what has previously been observed [51]. Even then, the relative proportion among all found miRNA target sites was generally low, unlike for example for ciRS-7, where 27.5% of all miRNA target sites belong to just a single miRNA, namely hsa-miR-7. Therefore, contrary to what has been described in many other publications, our analyses suggest few circRNAs to be effective miRNA binders. This is in line with more critical studies that consider miRNA sponging as a general function of circRNAs controversially [52]. One exception was CDR1as which we clearly found as a miRNA-binding circRNA [22]. Moreover, we newly identify circATP9B(5, 6, 7) as a putative miRNA sponge for miR-6794-3p. Of note, the database miRTarBase [53] lists more than 110 putative targets of miR-6794-3p. Although all listed targets are categorized as less strong evidence, they include for instance ZNF518B and FOXI2, which are associated with transcriptional regulation [54,55] and could be responsive to miR-6794-3p sponging by circATP9B(5, 6, 7). Further analyses, such as miRNA expression profiling or integration with AGO2 binding information, would be required to substantiate whether circATP9B(5, 6, 7) indeed acts as miRNA sponge.

For the RBP binding-site prediction, we saw equal numbers of circRNAs with at least one predicted RBP binding site around the BSJ and on the circRNA sequence. This might be related to RBP involvement in circRNA biogenesis [13,56–58] as well as regulation of circRNA export from the nucleus [59]. In line with a previous report [39], the main RBP detected around the BSJ was HuR. Note that since FIMO was chosen as a motif-based approach, only RBP binding sites present in the dataset [45] can be detected. However, additional motifs can be added to Calcifer if available.

The prediction of ORFs on circRNAs showed that many complete ORFs already exist in the linear sequences, but a considerable number of ORFs is present only in the circularized sequence. Together with the ORFs, which become accessible through multiple cycles of translation of the same circRNA, there is a large number of ORFs that uniquely exist on circRNAs alone. We investigated these ORFs for opportunities for cap-independent translation. For circRNAs, it has been shown that m6A RNA modifications and IRES-like sequences can enable translation and we found these possible motifs in front of nearly 50% of all non-linear ORFs [24,26]. In addition, we could show that non-linear ORFs could also encode unique peptides. Of particular note is the large amount of ORFs and associated unique peptides, especially in the context that circRNAs with complete ORFs are depleted in the nucleus compared to the entire cell. This further supports the possibility that these unique peptides might be translated from circRNAs and might have particular functions. It should be mentioned that a broad translation of circRNAs is still viewed critically, in particular with regard to scientifically proven evidence. This includes the facts that circRNA identification can be noisy, and that ribosome profiling and proteomics analyses have in some cases failed to find any evidence of circRNA translation products [60,61]. Nonetheless, there are several circRNAs for which protein products have been well established, as well as in silico experiments showing direct translation of circRNAs [24,37,62]. In any case, further research is needed to investigate the possible role of circRNA translation in a scientifically critical way and, if present, to detect and functionally analyse the resulting proteins.

CircRNAs are less abundant in the nucleus and particularly depleted from nuclear speckles

The first result of our analyses was that the vast majority of circRNAs are less abundant in the subnuclear compartments compared to the whole cell. This is in line with previous studies showing that circRNAs are also present in the nucleus, but are primarily exported to the cytoplasm [7,9,63]. Recent studies showed that the export of a subset of exonic circRNAs might be linked to XPO4 [64] and that the export of circRNAs might be dependent on their length [19]. It was also shown that circRNAs bound by SRSF1, which is the most detected RBP binding site within our circRNAs, are more likely to remain in the nucleus rather than being exported into the cytoplasm [59]. The interplay of such mechanisms could potentially explain the different distribution of circRNAs between the nucleus and cytoplasm. Even though new mechanisms for circRNA export are constantly being described [65], a more comprehensive knowledge of the various mechanisms by which circRNAs are exported from the nucleus still remains an important question for future research.

Remarkably, we observed a significant depletion of circRNAs from nuclear speckles. These are membrane-less organelles that are strongly enriched for core splicing factors, including SR proteins. Previous studies showed that nuclear speckles might play a crucial role in mRNA splicing as well as the general regulation of gene expression [47,66–68]. Our data indicate that the excess of splicing factors in this compartment disfavours circRNA formation, possibly by promoting canonical linear splicing. In addition, other factors like the release timing of circRNAs might also add to the general depletion in nuclear speckles. Consistent with the idea of disfavoured circRNA formation by SRSF proteins, we could recently show that overexpression of the SR protein SRSF6 almost completely abolished circRNA formation in response to hypoxia human HeLa cells [48]. Although the dilution of circRNAs in SRSF6 overexpression might be favoured by the increased proliferation [17,18,69], we suggest that the absence of many abundant circRNAs may not be explained by this alone and be directly linked to SRSF6 promoting linear splicing.

The lncRNA FIRRE is mostly circularized in the cytoplasm and also to a lower degree in the nuclear compartments

Among the specific circRNAs, we found circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) as an interesting candidate that is abundantly expressed across nuclear compartments in the APEX-Seq data. Our data indicate that circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) is also efficiently exported to the cytoplasm where it accumulates. It is known that the FIRRE locus produces a lncRNA as well as circRNAs that are expressed during ES cell differentiation [50,70]. The lncRNA FIRRE plays a role in chromosomal localization in the nucleus [71,72], in adipogenesis in mouse [73], in development of different cancer types [74–76] and as trans-acting RNA in haematopoiesis [50]. We speculate that some of these functions could be modulated by the circular isoform circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8). Of note, this major circFIRRE isoform contains several unannotated exons that do not occur in the linear counterpart, as already observed in other circRNA isoforms of FIRRE detected in human brain tissue [77], and could therefore be used for specifically targeting circFIRRE(3, NE, NE, 4, NE, 5, 6, 7, 8) to investigate possible functions.

Conclusion

In this study, we showcased the workflow Calcifer for circRNA detection and analysis. We found a high amount of circRNAs in the whole cell but also in the subnuclear compartments. Surprisingly, many of the circRNAs were predicted to encode for complete ORFs, suggesting that a large number of peptides might be produced from these circRNAs. Experimental approaches such as mass spectrometry will be critical to evaluate these findings. Comparing between the subnuclear compartments, we observed a significant depletion of circRNAs from nuclear speckles, supporting the notion that core splicing factors such as SR proteins generally inhibit circRNA formation.

Methods

Initial quality control

We checked the overall quality of all APEX-Seq data using FastQC (version 0.11.9) [78].

Calcifer workflow

RNA-seq preprocessing

For the initial trimming of the sequencing reads, Flexbar (version 3.5.0) was utilized with default settings and –min-read-length 20 [79]. All remaining reads were then used for circRNA detection.

CircRNA detection and quantification

The Calcifer pipeline employs the two established tools CIRCexplorer2 (version 2.3.8) and CIRI2 (version 2.0.6) [32,80]. For CIRCexplorer2, RNA-seq reads were mapped against the reference genome (ENSEMBL release 108, GRCh38.108) using STAR (version 2.7.6a) with the following parameters: –outFilterMultimapNmax 1 –outFilterMismatchNmax 2 –alignSJDBoverhangMin 15 –alignSJoverhangMin 15 –chimSegmentMin 15 –chimScoreMin 15 –chimScoreSeparation 10 –chimJunctionOverhangMin 15 [35]. The chimeric junctions from STAR were used as input for CIRCexplorer2. For CIRI2, RNA-seq reads were aligned against the same reference genome using bwa (version 0.7.17) with the parameter -T 19 to pre-filter the resulting alignments [36]. CIRCexplorer2 and CIRI2 were both run with default parameters. For CIRCexplorer2, we performed a baseline filtering for circRNAs with at least 2 BSJ reads, which was also the default for CIRI2. The remaining circRNAs detected by these tools were then combined in one set and further filtered.

CircRNAs that lacked canonical splice sites or were longer than 100 kb (genomic distance between back-splice sites) were filtered out by the workflow as described before [39]. To check for canonical splice sites, we retrieved the 4 nt spanning over each back-splice site using samtools (version 1.17) faidx [81]. To merge the results between both tools and achieve comparable quantifications, we recounted the supporting BSJ reads for each circRNA from the chimeric alignments reported by STAR. Only circRNAs that had at least two unique mapped reads (i.e. with distinct CIGAR strings) over the BSJ were kept. If there were several replicas, a sum of 2 unique BSJ reads between two or more replicas was also sufficient. The same mapping of two reads in two different replicates counted as separate BSJ reads in the final results. This filtering step was done via comparing the CIGAR strings of all chimerically mapped reads for each BSJ. We used PyEnsembl (version 1.1.0) to retrieve the parental gene names for the Ensembl gene IDs reported by CIRCexplorer2 and CIRI2 (https://github.com/openvax/pyensembl).

Computation of the circular to linear ratio (CLR) for each circRNA

For each circRNA, we used the BSJ read count and divided it by the BSJ read count plus the linear splice junction read counts from the back-splice sites:CLR=circularcountlinearcount+circularcount

With a linear count of 0, this results in a CLR of 1, which indicates a complete circularization. A percent circularized value can be computed by multiplying the CLR with 100.

miRNA target site prediction on circRNAs

To detect potential miRNA target sites, the exonic sequence was determined for each circRNA. This was done based on gene annotation, taking into account all exons between the back-splice sites. To allow for target sites across the BSJ, we added 25 nt across the BSJ to both ends (‘pseudo-circular’). In order to predict miRNA target sites, miRanda (version 3.3a) was used [42]. The parameters for miRanda were chosen as suggested previously [39]. With -strict a perfect alignment in the seed region is required and with -sc 150 only hits with an alignment score ≥ 150 are considered. As miRNA input, all human miRNA sequences in fasta format from miRbase (v22) were utilized (accessed on 17 May 2021) [82]. For the general circRNA output we report the miRNA with the most target sites. In addition, we calculate the proportion of miRNA with the most target sites to the total detected target sites on each given circRNA:ProportionoftotalmiRNAtargetsites=maximumtargetsitessinglemiRNAtotalmiRNAtargetsites

To determine putative miRNA sponges that show a certain enrichment for target sites of one individual miRNA, we applied a cut-off of 15% for the proportion of total miRNA target sites. This cut-off is variable and can be changed by the user.

RBP binding site prediction on circRNAs and around BSJs

RBP binding site prediction was performed on the pseudo-circular circRNA sequence (i.e. adding 25 nt over the BSJ to both ends) as well as in a fixed window around the back-splice sites including 250 nt sequence into the intron until 25 nt into the exon. The RBP motifs (102 motifs, 7–8 nt width) were predicted on the circRNAs with FIMO (MEME version 5.4.1) and were taken from the MEME motif database [44,45]. FIMO was run with default parameters, and all predicted binding sites were used.

Open reading frame (ORF) prediction on circRNAs

Putative ORFs in circRNAs were predicted with Calcifer using a custom algorithm. The linear sequence containing all exons between the back-splice sites was used to determine putative ORFs. Fasta files with (a) the linear sequence containing all exons between the back-splice sites, as well as (b) the same sequence with added sequence (25 nt) over the BSJ (‘pseudo-circular’) and (c) a concatenation of up to four copies of the circRNA sequence were used as input. The latter is necessary to account for possible ORFs that run several times around the circRNA due to frameshifts or rolling circle translation as described in previous research [83]. The ORFs were detected via translation of the circRNA sequence in all different reading frames, searching for peptides with at least 10 amino acids length. Therefore, an ORF was defined by a start codon, the peptide sequence and a stop codon. Only the first appearance for each distinct ORF was considered in the further analyses.

Detection of m6A/IRES-like sites upstream of circRNA ORFs

Internal ribosomal entry sites (IRES) are the primary candidates for the translation of circRNAs. As previously described, m6A sites and also shorter so-called IRES-like sequences can initiate translation on circRNAs [24,26]. To find possible m6A and IRES-like sites, we examined 10 nt upstream of each complete ORF start codon on the circRNAs. For these sequences we generated a dictionary with all 5-mer (DRACH motifs) and 6-mer (IRES-like motifs, as described previously) [24] sequences and the circRNAs where they occur. For each of these relevant 5-mer and 6-mer sequences, we checked the dictionary for the occurrences and could thus determine all complete ORFs with these motifs on the circRNAs.

Determine uniqueness of putative circRNA-encoded peptides

For each complete ORF on a circRNA, with at least one putative m6A or IRES-like site, we determined the possible resulting peptide. For further analysis and a possible future detection in mass spectrometry data, it was necessary to screen these peptides for unique segments. For this, we analysed the overlaps of the peptides with the whole proteome of the organism. Peptides or peptide fragment with a minimum length of 10 amino acids that could not be detected in the proteome with MUMmer (version 3.23) were reported as unique [84].

DESeq2 analysis of differentially expressed genes and circRNAs

For DESeq2 (version 1.34.0), raw gene counts were determined via htseq-count (version 0.11.3) [85]. The raw circRNA counts were calculated by Calcifer. All raw counts were analysed together with DESeq2 to obtain the normalized counts for each circRNA and gene for each replicate [86]. DESeq2 was used without independent filtering for the differential expression analysis. Instead, only genes and circRNAs with a normalized mean count >1 were considered. The heatmaps were created with ComplexHeatmap (version 2.10.0) and the z-score was computed based on the normalized counts between the replicates [87].

Statistical analysis

The testing for statistical differences between two groups was performed with Welch Two Sample t-test. Statistical differences between two proportions from nuclear and whole-cell circRNAs were tested with Fisher exact test.

Highlights

Calcifer offers an all-round workflow for circRNA analysis including predictions of miRNA target sites, RBP binding, and translation.

CircRNAs differ between subnuclear compartments and are strongly depleted from nuclear speckles.

The circRNA circCANX(9) is uniquely enriched in the nucleus.

The lncRNA FIRRE predominantly occurs as a circRNA in the cytoplasm and the nucleus.

Supplementary Material

Supplemental Material

Acknowledgments

The authors would like to thank all members of the Zarnack group for valuable discussion as well as Antonella Di Liddo for initial support.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data for the initial runtime test is published [88] and online available in the SRA under SRR3420582, SRR3420583, SRR3420584, SRR3420586 and SRR3420587.

The APEX-Seq dataset is published [34] and online available under GSE176439 in the NCBI GEO database. The latest version of Calcifer is available on GitHub (https://github.com/ZarnackGroup/Calcifer).

The HeLa SRSF6 overexpression [48] and control [39] data is published and available in the NCBI GEO database under GSE198308 and GSE131379 respectively.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15476286.2024.2395718
==== Refs
References

[1] Lee ECS, Elhassan SAM, Lim GPL, et al., & others. The roles of circular RNAs in human development and diseases. Biomed & Pharmacother. 2019;111 :198–208. doi: 10.1016/j.biopha.2018.12.052
[2] Shao T, Pan Y, Xiong X. Circular RNA: an important player with multiple facets to regulate its parental gene expression. Mol Ther Nucleic Acids. 2021;23 :369–376. doi: 10.1016/j.omtn.2020.11.008 33425494
[3] Rybak-Wolf A, Stottmeister C, Glažar P, et al., & others. (Circular RNAs in the mammalian brain are highly abundant, conserved, and dynamically expressed. Mol Cell. 2015;58 (5 ):870–885. doi: 10.1016/j.molcel.2015.03.027 25921068
[4] Salzman J, Chen RE, Olsen MN, et al. Cell-type specific features of circular RNA expression. PLOS Genet. 2013;9 (9 ):e1003777. doi: 10.1371/journal.pgen.1003777 24039610
[5] Zhang F, Jiang J, Qian H, et al. Exosomal circRNA: emerging insights into cancer progression and clinical application potential. J Hematol Oncol. 2023;16 (1 ):67. doi: 10.1186/s13045-023-01452-2 37365670
[6] Zhao H, Tan Z, Zhou J, et al. The regulation of circRNA and lncRNA protein binding in cardiovascular diseases: emerging therapeutic targets. Biomed & Pharmacother. 2023;165 :115067. doi: 10.1016/j.biopha.2023.115067
[7] Capel B, Swain A, Nicolis S, et al. Circular transcripts of the testis-determining gene sry in adult mouse testis. Cell. 1993;73 (5 ):1019–1030. doi: 10.1016/0092-8674(93)90279-Y 7684656
[8] Hansen TB, Wiklund ED, Bramsen JB, et al. miRNA-dependent gene silencing involving Ago2-mediated cleavage of a circular antisense RNA. The EMBO Journal. 2011;30 (21 ):4414–4422. doi: 10.1038/emboj.2011.359 21964070
[9] Memczak S, Jens M, Elefsinioti A, et al., & others. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495 (7441 ):333–338. doi: 10.1038/nature11928 23446348
[10] Salzman J, Gawad C, Wang PL, et al. Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLOS ONE. 2012;7 (2 ):e30733. doi: 10.1371/journal.pone.0030733 22319583
[11] Zhang Y, Xue W, Li X, et al. The biogenesis of nascent circular RNAs. Cell Reports. 2016;15 (3 ):611–624. doi: 10.1016/j.celrep.2016.03.058 27068474
[12] Ashwal-Fluss R, Meyer M, Pamudurti NR, et al. circRNA biogenesis competes with pre-mRNA splicing. Molecular Cell. 2014;56 (1 ):55–66. doi: 10.1016/j.molcel.2014.08.019 25242144
[13] Conn SJ, Pillman KA, Toubia J, et al. The RNA binding protein quaking regulates formation of circRNAs. Cell. 2015;160 (6 ):1125–1134. doi: 10.1016/j.cell.2015.02.014 25768908
[14] Liang D, Wilusz JE. Short intronic repeat sequences facilitate circular RNA production. Genes & Devel. 2014;28 (20 ):2233–2247.25281217
[15] Zhang X-O, Dong R, Zhang Y, et al. Diverse alternative back-splicing and alternative splicing landscape of circular RNAs. Genome Res. 2016;26 (9 ):1277–1287. doi: 10.1101/gr.202895.115 27365365
[16] Starke S, Jost I, Rossbach O, et al. Exon circularization requires canonical splice signals. Cell Reports. 2015;10 (1 ):103–111. doi: 10.1016/j.celrep.2014.12.002 25543144
[17] Bachmayr-Heyda A, Reiner AT, Auer K, et al. Correlation of circular RNA abundance with proliferation – exemplified with colorectal and ovarian cancer, idiopathic lung fibrosis and normal human tissues. Sci Rep. 2015;5 (1 ):8057. doi: 10.1038/srep08057 25624062
[18] García-Rodríguez JL, Korsgaard U, Ahmadov U, et al. Spatial profiling of circular RNAs in cancer reveals high expression in muscle and stromal cells. Cancer Research. 2023;83 (20 ):3340–3353. doi: 10.1158/0008-5472.CAN-23-0748 37477923
[19] Huang C, Liang D, Tatomer DC, et al. A length-dependent evolutionarily conserved pathway controls nuclear export of circular RNAs. Genes Dev. 2018;32 (9–10 ):639–644. doi: 10.1101/gad.314856.118 29773557
[20] Du WW, Yang W, Liu E, et al. Foxo3 circular RNA retards cell cycle progression via forming ternary complexes with p21 and CDK2. Nucleic Acids Res. 2016;44 (6 ):2846–2858. doi: 10.1093/nar/gkw027 26861625
[21] Du WW, Zhang C, Yang W, et al. Identifying and characterizing circRNA-protein interaction. Theranostics. 2017;7 (17 ):4183. doi: 10.7150/thno.21299 29158818
[22] Hansen TB, Jensen TI, Clausen BH, et al. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495 (7441 ):384–388. doi: 10.1038/nature11993 23446346
[23] Lasda E, Parker R. Circular RNAs: diversity of form and function. RNA. 2014;20 (12 ):1829–1842. doi: 10.1261/rna.047126.114 25404635
[24] Fan X, Yang Y, Chen C, et al. Pervasive translation of circular RNAs driven by short ires-like elements. Nat Commun. 2022;13 (1 ):3751. doi: 10.1038/s41467-022-31327-y 35768398
[25] Shi Y, Jia X, Xu J. The new function of circRNA: translation. Clin Transl Oncol. 2020;22 (12 ):2162–2169. doi: 10.1007/s12094-020-02371-1 32449127
[26] Yang Q, Du WW, Wu N, et al., & others. A circular RNA promotes tumorigenesis by inducing c-myc nuclear translocation. Cell Death Differ. 2017;24 (9 ):1609–1620. doi: 10.1038/cdd.2017.86 28622299
[27] Li Z, Huang C, Bao C, et al. Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol. 2015;22 (3 ):256–264. doi: 10.1038/nsmb.2959 25664725
[28] Wang L, Long H, Zheng Q, et al. Circular RNA circRHOT1 promotes hepatocellular carcinoma progression by initiation of NR2F6 expression. Mol Cancer. 2019;18 (1 ):1–12. doi: 10.1186/s12943-019-1046-7 30609930
[29] Fazal FM, Han S, Parker KR, et al. Atlas of subcellular RNA localization revealed by APEX-Seq. Cell. 2019;178 (2 ):473–490. doi: 10.1016/j.cell.2019.05.027 31230715
[30] Chen L, Wang C, Sun H, et al. The bioinformatics toolbox for circRNA discovery and analysis. Briefings In Bioinformatics. 2021;22 (2 ):1706–1728. doi: 10.1093/bib/bbaa001 32103237
[31] Hansen TB. Improved circRNA identification by combining prediction algorithms. Front Cell Dev Biol. 2018;6 :20. doi: 10.3389/fcell.2018.00020 29556495
[32] Gao Y, Zhang J, Zhao F. Circular RNA identification based on multiple seed matching. Briefings In Bioinformatics. 2018;19 (5 ):803–810. doi: 10.1093/bib/bbx014 28334140
[33] Jakobi T, Uvarovskii A, Dieterich C. Circtools—a one-stop software solution for circular RNA research. Bioinformatics. 2019;35 (13 ):2326–2328. doi: 10.1093/bioinformatics/bty948 30462173
[34] Barutcu AR, Wu M, Braunschweig U, et al. (Systematic mapping of nuclear domain-associated transcripts reveals speckles and lamina as hubs of functionally distinct retained introns. Mol Cell. 2022;82 (5 ):1035–1052. doi: 10.1016/j.molcel.2021.12.010 35182477
[35] Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29 (1 ):15–21. doi: 10.1093/bioinformatics/bts635 23104886
[36] Li H, Durbin R. Fast and accurate short read alignment with burrows–Wheeler transform. Bioinformatics. 2009;25 (14 ):1754–1760. doi: 10.1093/bioinformatics/btp324 19451168
[37] Chen C-K, Cheng R, Demeter J, et al. Structured elements drive extensive circular RNA translation. Mol Cell. 2021;81 (20 ):4300–4318.e13. doi: 10.1016/j.molcel.2021.07.042 34437836
[38] Hansen TB, Venø MT, Damgaard CK, et al. Comparison of circular RNA prediction tools. Nucleic Acids Res. 2016;44 (6 ):e58–e58. doi: 10.1093/nar/gkv1458 26657634
[39] Di Liddo A, de Oliveira Freitas Machado C, Fischer S, et al. A combined computational pipeline to detect circular RNAs in human cancer cells under hypoxic stress. J Mol Cell Biol. 2019;11 (10 ):829–844. doi: 10.1093/jmcb/mjz094 31560396
[40] Glažar P, Papavasileiou P, Rajewsky N. circBase: a database for circular RNAs. RNA. 2014;20 (11 ):1666–1670. doi: 10.1261/rna.043687.113 25234927
[41] Maass PG, Glažar P, Memczak S, et al. A map of human circular RNAs in clinically relevant tissues. J Mol Med (Berl). 2017;95 (11 ):1179–1189. doi: 10.1007/s00109-017-1582-9 28842720
[42] Enright A, John B, Gaul U, et al. MicroRNA targets in Drosophila. Genome Biol. 2003;5 (1 ):1–27. doi: 10.1186/gb-2003-5-1-r1
[43] Chen L-L, Bindereif A, Bozzoni I, et al. A guide to naming eukaryotic circular RNAs. Nat Cell Biol. 2023;25 (1 ):1–5. doi: 10.1038/s41556-022-01066-9 36658223
[44] Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011;27 (7 ):1017–1018. doi: 10.1093/bioinformatics/btr064 21330290
[45] Ray D, Kazan H, Cook KB, et al. A compendium of RNA-binding motifs for decoding gene regulation. Nature. 2013;499 (7457 ):172–177. doi: 10.1038/nature12311 23846655
[46] Sun P, Li G. CircCode: a powerful tool for identifying circRNA coding ability. Front Genet. 2019;10 . doi: 10.3389/fgene.2019.00981
[47] Spector DL, Lamond AI. Nuclear speckles. Cold Spring Harbor Perspectives In Biology. 2011;3 (2 ):a000646. doi: 10.1101/cshperspect.a000646 20926517
[48] de Oliveira Freitas Machado C, Schafranek M, Brüggemann M, et al. Poison cassette exon splicing of SRSF6 regulates nuclear speckle dispersal and the response to hypoxia. Nucleic Acids Res. 2023;51 (2 ):870–890. doi: 10.1093/nar/gkac1225 36620874
[49] Song Y-Z, Li J-F. Circular RNA hsa_circ_0001564 regulates osteosarcoma proliferation and apoptosis by acting miRNA sponge. Biochem Biophys Res Commun. 2018;495 (3 ):2369–2375. doi: 10.1016/j.bbrc.2017.12.050 29229385
[50] Lewandowski JP, Lee JC, Hwang T, et al. The Firre locus produces a trans-acting RNA molecule that functions in hematopoiesis. Nat Commun. 2019;10 (1 ):5137. doi: 10.1038/s41467-019-12970-4 31723143
[51] Guo JU, Agarwal V, Guo H, et al. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014;15 (7 ):409. doi: 10.1186/s13059-014-0409-z 25070500
[52] Jarlstad Olesen MT, Kristensen LS. Circular RNAs as microRNA sponges: evidence and controversies. Essays Biochem. 2021;65 (4 ):685–696. doi: 10.1042/EBC20200060 34028529
[53] Huang H-Y, Lin Y-C-D, Cui S, et al. miRtarbase update 2022: an informative resource for experimentally validated miRNA–target interactions. Nucleic Acids Res. 2022;50 (D1 ):D222–D230. doi: 10.1093/nar/gkab1079 34850920
[54] Riffo-Campos ÁL, Castillo J, Vallet-Sanchez A, et al. In silico RNA-seq and experimental analyses reveal the differential expression and splicing of EPDR1 and ZNF518B genes in relation to KRAS mutations in colorectal cancer cells. Oncology Reports. 2016;36 (6 ):3627–3634. doi: 10.3892/or.2016.5210 27805251
[55] Wijchers PJEC, Hoekman MFM, Burbach JPH, et al. Cloning and analysis of the murine Foxi2 transcription factor. Biochim Et Biophys Acta (BBA)-Gene Struct And Expression. 2005;1731 (2 ):133–138. doi: 10.1016/j.bbaexp.2005.09.003
[56] Errichelli L, Dini Modigliani S, Laneve P, et al. FUS affects circular RNA expression in murine embryonic stem cell-derived motor neurons. Nat Commun. 2017;8 (1 ):14741. doi: 10.1038/ncomms14741 28358055
[57] Teplova M, Hafner M, Teplov D, et al. Structure–function studies of STAR family quaking proteins bound to their in vivo RNA target sites. Genes Dev. 2013;27 (8 ):928–940. doi: 10.1101/gad.216531.113 23630077
[58] Zhong Y, Yang Y, Wang X, et al. Systematic identification and characterization of exon-intron circRNAs. Genome Res. 2024;34 (3 ):376–393. doi: 10.1101/gr.278590.123 38609186
[59] Ron M, Ulitsky I. Context-specific effects of sequence elements on subcellular localization of linear and circular RNAs. Nat Commun. 2022;13 (1 ):2481. doi: 10.1038/s41467-022-30183-0 35513423
[60] Hansen TB. Signal and noise in circRNA translation. Methods. 2021;196 :68–73. doi: 10.1016/j.ymeth.2021.02.007 33588029
[61] Stagsted LV, Nielsen KM, Daugaard I, et al. Noncoding AUG circRNAs constitute an abundant and conserved subclass of circles. Life Sci Alliance. 2019;2 (3 ):e201900398. doi: 10.26508/lsa.201900398 31028096
[62] Wen S, Qadir J, Yang BB. Circular RNA translation: novel protein isoforms and clinical significance. Trends Mol Med. 2022;28 (5 ):405–420. doi: 10.1016/j.molmed.2022.03.003 35379558
[63] Werfel S, Nothjunge S, Schwarzmayr T, et al. Characterization of circular RNAs in human, mouse and rat hearts. J. Mol cell Cardiology. 2016;98 :103–107. doi: 10.1016/j.yjmcc.2016.07.007
[64] Chen L, Wang Y, Lin J, et al. Exportin 4 depletion leads to nuclear accumulation of a subset of circular RNAs. Nat Commun. 2022;13 (1 ):5769. doi: 10.1038/s41467-022-33356-z 36182935
[65] Ngo LH, Bert AG, Dredge BK, et al. Nuclear export of circular RNA. Nature. 2024;627 (8002 ):212–220. doi: 10.1038/s41586-024-07060-5 38355801
[66] Bhat P, Chow A, Emert B, et al. Genome organization around nuclear speckles drives mRNA splicing efficiency. Nature. 2024;629 (8014 ):1165–1173. doi: 10.1038/s41586-024-07429-6 38720076
[67] Chen Y, Belmont AS. Genome organization around nuclear speckles. Curr Opin In Genet & Devel. 2019;55 :91–99. doi: 10.1016/j.gde.2019.06.008 31394307
[68] Galganski L, Urbanek MO, Krzyzosiak WJ. Nuclear speckles: molecular organization, biological function and role in disease. Nucleic Acids Research. 2017;45 (18 ):10350–10368. doi: 10.1093/nar/gkx759 28977640
[69] She W, Shao J, Jia R. Targeting splicing factor SRSF6 for cancer therapy. Front Cell Dev Biol. 2021;9 . doi: 10.3389/fcell.2021.780023
[70] Izuogu OG, Alhasan AA, Mellough C, et al. (Analysis of human ES cell differentiation establishes that the dominant isoforms of the lncRNAs RMST and FIRRE are circular. BMC Genomics. 2018;19 (1 ):1–18. doi: 10.1186/s12864-018-4660-7 29291715
[71] Hacisuleyman E, Goff LA, Trapnell C, et al. Topological organization of multichromosomal regions by the long intergenic noncoding RNA firre. Nat Struct & Mol Biol. 2014;21 (2 ):198–206. doi: 10.1038/nsmb.2764 24463464
[72] Yang F, Deng X, Ma W, et al. (De Novo assembly of bacterial transcriptomes from RNA-seq data. Genome Biol. 2015;16 (1 ):1–17. doi: 10.1186/s13059-014-0572-2 25583448
[73] Sun L, Goff LA, Trapnell C, et al. (Long noncoding RNAs regulate adipogenesis. Proc Natl Acad Sci USA. 2013;110 (9 ):3387–3392. doi: 10.1073/pnas.1222643110 23401553
[74] Liu Q-H, Dai G-R, Wu Y, et al. LncRNA FIRRE stimulates PTBP1-induced Smurf2 decay, stabilizes B-cell receptor, and promotes the development of diffuse large B-cell lymphoma. Hematological Oncology. 2022;40 (4 ):554–566. doi: 10.1002/hon.3004 35416325
[75] Wang Y, Li Z, Xu S, et al. LncRNA FIRRE functions as a tumor promoter by interaction with PTBP1 to stabilize BECN1 mRNA and facilitate autophagy. Cell Death Dis. 2022;13 (2 ):98. doi: 10.1038/s41419-022-04509-1 35110535
[76] Zhou J, Liu T, Xu H, et al. LncRNA FIRRE promotes the proliferation and metastasis of hepatocellular carcinoma by regulating the expression of PXN through interacting with MBNL3. Biochem Biophys Res Commun. 2022;625 :188–195. doi: 10.1016/j.bbrc.2022.07.099 35988459
[77] Rahimi K, Venø MT, Dupont DM, et al. Nanopore sequencing of brain-derived full-length circRNAs reveals circRNA-specific exon usage, intron retention and microexons. Nat Commun. 2021;12 (1 ):4825. doi: 10.1038/s41467-021-24975-z 34376658
[78] Andrews S& others. Babraham bioinformatics-FastQC a quality control tool for high throughput sequence data. 2010. https://Www.Bioinformatics.Babraham.Ac.Uk/Projects/Fastqc
[79] Roehr JT, Dieterich C, Reinert K. Flexbar 3.0–SIMD and multicore parallelization. Bioinformatics. 2017;33 (18 ):2941–2942. doi: 10.1093/bioinformatics/btx330 28541403
[80] Zhang H, Jiang L, Sun D, et al. CircRNA: a novel type of biomarker for cancer. Breast Cancer. 2018;25 (1 ):1–7. doi: 10.1007/s12282-017-0793-9 28721656
[81] Danecek P, Bonfield JK, Liddle J, et al. (Twelve years of SAMtools and BCFtools. Gigascience. 2021;10 (2 ):giab008. doi: 10.1093/gigascience/giab008 33590861
[82] Griffiths-Jones S, Grocock RJ, Van Dongen S, et al. miRbase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34 (suppl_1 ):D140–D144. doi: 10.1093/nar/gkj112 16381832
[83] Abe N, Matsumoto K, Nishihara M, et al. Rolling circle translation of circular RNA in living human cells. Sci Rep. 2015;5 (1 ):16435. doi: 10.1038/srep16435 26553571
[84] Kurtz S, Phillippy A, Delcher AL, et al. Versatile and open software for comparing large genomes. Genome Biol. 2004;5 (2 ):1–9. doi: 10.1186/gb-2004-5-2-r12
[85] Anders S, Pyl PT, Huber W. Htseq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31 (2 ):166–169. doi: 10.1093/bioinformatics/btu638 25260700
[86] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15 (12 ):1–21. doi: 10.1186/s13059-014-0550-8
[87] Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32 (18 ):2847–2849. doi: 10.1093/bioinformatics/btw313 27207943
[88] Chujo T, Yamazaki T, Kawaguchi T, et al. Unusual semi‐extractability as a hallmark of nuclear body‐associated architectural noncoding RNAs. Embo J. 2017;36 (10 ):1447–1462. doi: 10.15252/embj.201695848 28404604
