Comment on "A comprehensive overview and evaluation of circular RNA detection tools"

A recent paper published in PLOS Computational Biology [1] provided a comprehensive evaluation of various circular RNA (circRNA)-detection tools. The authors compared 11 different circRNA-detection tools using four different datasets, including three simulated datasets (positive, background, and mixed datasets) and one real dataset. Since the advent of highthroughput next-generation sequencing technology, dozens of computational tools have been developed and used to successfully detect thousands of circRNAs in a diverse range of species. However, there are great discrepancies in the results obtained using different tools [2–7], and systematic evaluations of their performance have not been available. Indeed, the cited work has provided a useful guideline for researchers engaged in circRNA studies. Nevertheless, it seems inappropriate to use all CircBase-deposited circRNA candidates (14,689 events) identified in silico from RNA-seq data of HeLa cells [8] as the positive dataset. The qualification of the 14,689 candidates requires further evaluation. We suggest that three main confounding factors, which may affect the fairness of the evaluation of circRNA-detection tools, should be considered. First, it has been shown that non-co-linear (NCL) junctions (including circRNA and transspliced RNA junctions) that do not match annotated exon boundaries tend to be unreliable and are more likely to stem from mis-splicing [9–12], although we cannot eliminate the possibility that a few true backspliced junctions indeed originate from unannotated gene loci. Since circRNA candidates are regarded to be less or more reliable if their normalized read counts are depleted or enriched after RNase R treatment, respectively [13], we reexamined the circRNA candidates detected on the HeLa RNase R-treated and untreated samples (the circRNA candidates and the corresponding read counts were downloaded from the cited study). Of the circRNA candidates with unannotated exon boundaries, we can find that 50%~100% of them were “completely” depleted (not detected) after RNase R treatment, whereas only <8% of them were “significantly” enriched (i.e., 5-fold increase in normalized read count) after RNase R treatment (Fig 1). This result revealed that the candidates with unannotated exon boundaries are more likely to be false calls. Thus, we suggest that the CircBase circRNA candidates with unannotated exon boundaries (1,046 events; Table 1) should be excluded from the positive dataset. At least, since circRNA junctions were observed to be predominantly located at canonical splice sites [14–16], the candidates with junctions that have not canonical splice site sequences (GT-AG, GC-AG, or AT-AC) should be removed (778 events; Table 1). Second, ambiguous alignments originating from repetitive sequences or paralogous genes often result in false positive circRNA detection. In CircBase, most circRNA candidates were identified by find_circ [8]. It has been reported that some of find_circ-identified candidates were mis-predicted from paralogous genes [17]. Therefore, the factor of alignment ambiguity should be considered when using CircBase circRNAs as true positives. To this end, we concatenated the exonic sequence flanking the circRNA junction (within -100 nucleotides

to +100 nucleotides of each CircBase-identified junction; see Fig 2A). We then aligned the 200 bp concatenated sequence (the concatenated sequence may be shorter than 200 bp if the exonic circRNA sequence is shorter than 200 bp) against the reference genome and NCBI RefSeq-/GENCODE-identified mRNAs using BLAT [18]. It is worth noting that 2,316 concatenated sequences exhibited ambiguous alignments (Table 1), of which 514 events contained alternative co-linear explanations (an example was illustrated in Fig 2B) and 1,802 events mapped to multiple positions with similar BLAT mapping scores (difference of mapping scores < 5) (an example was illustrated in Fig 2C). These 2,316 CircBase circRNAs are not appropriate for inclusion in the positive dataset. We further examined the 2,316 potentially Only the circRNAs with �2 supporting NCL junction reads without RNase R treatment were considered. Of note, no circRNAs with unannotated exon boundaries were detected by UROBORUS, KNIFE, and NCLscan on the HeLa samples without RNase R treatment. false calls and found that they contributed to 4%~14% of the identifications among the 11 tools ( Fig 2D). Although these false calls can be considerably eliminated by merging two cir-cRNA predicting tools, the filtering performance is dependent on tools used ( Fig 2E). These results also suggested that these false calls were not the find_circ specific issue in CircBase.
Third, of the CircBase circRNAs, 3,580 events were identified by all 11 tools. We observed that the number of identified backspliced junction reads varied remarkably between tools depending on the level of strictness of the filtering steps used (Fig 3A and a similar result shown in Fig. 4C of the cited study [1]). With relatively strict filtering criteria, UROBORUS, NCLscan, circRNA_finder, and CIRCexplorer provided a higher percentage of identified cir-cRNAs with one supporting junction read than did the other tools ( Fig 3B). These four tools and DCC provided a relatively low percentage (<10%) of identified circRNAs with >10 supporting junction read, whereas such percentages were more than 60% for CIRI and MapSplice ( Fig 3B). For example, for the CircBase circRNA candidate of COQ4 (Fig 3C, top), CIRI identified five supporting junction reads, whereas UROBORUS and NCLscan identified one supporting junction read only (the information of supporting junction reads were extracted from the github website generated by the cited study at https://github.com/linatbeishan/circRNA_ detection [1]). However, four of the five CIRI-provided junction reads were unqualified reads with ambiguous alignments of an alternative co-linear explanation (Fig 3C, middle and bottom), if realigned against the reference genome and annotated mRNAs. Such unqualified reads were not included in the UROBORUS and NCLscan results. To evaluate the impact of the unqualified reads on circRNA predictions, we took the result of CIRI identification as an example (because only CIRI provided the IDs of supporting backspliced junction reads). We found that as high as 50% (57,288 reads) of the CIRI-provided backspliced junction reads were unqualified (Fig 3D, left). After removing these unqualified reads, the number of the detected circRNA candidates (Fig 3D, left) and the median number of supporting backspliced junction reads per candidate (Fig 3D, middle) remarkably decreased; meanwhile, the percentage of the CIRI-identified circRNAs with one supporting junction read significantly increased (from 0% to 5%; Fig 3D, right). These results revealed that the unqualified reads would considerably affect the accuracy of circRNA predictions and should be carefully controlled. In addition, as also noted in the cited study [1], different tools use different counting methods to deal with a read pair spanning the same backspliced junction. Thus, it would be unfair to evaluate sensitivity among tools while directly considering candidates with � 2 supporting junction reads Table 1. Uncertain circRNA candidates in CircBase (see also S1 Dataset).

Confounding factors Number of circRNAs
Junctions not matching annotated exon boundaries (A) 1,046 Non-canonical splice site 1 778 Ambiguous alignment (B) 2 2,316 Alternative co-linear explanation 514 Multiple hits 1,802 Total uncertain circRNA candidates((A)[(B)) 3,110 1 The backspliced junctions have not canonical splice site sequences (GT-AG, GC-AG, or AT-AC). 2 The extracted concatenated sequences were aligned against the reference/annotated transcripts using BLAT with two parameter sets: (-titleSize = 9, -stepSize = 9, -repMatch = 32768; default parameters of the Windows version) and (-titleSize = 11, -stepSize = 11, -repMatch = 1024; default parameters of the stand-alone version). A concatenated sequence is regarded as an "alternative co-linear explanation" if its BLAT-alignment results contain at least one colinear explanation.  (i.e., Table 1 of the cited study [1]). By the same reason, it is improper to assess sensitivity among tools at the read level without accounting for read qualification (i.e., Fig. 4 of the cited study [1]). Identification of NCL RNAs is often hampered by false positives arising from sequencing/ alignment errors (particularly alignment ambiguity as stated above) and experimental artifacts (particularly in vitro artifacts arising from template switching during reverse transcription (RT)) [11,12,[19][20][21][22]. Selection of positive/negative datasets (positive datasets particularly) for evaluating accuracy of circRNA-detection tools should be careful. Generally, positive datasets are generated or selected from: (1) simulation, (2) real RNA-seq datasets from samples with different treatments, and (3) wet-laboratory validated datasets. For simulated datasets, the colinear (false positive or background) and circRNA (true positive) reads are generated from well-annotated co-linear mRNAs. Alternatively, in the cited study, the true positive reads were generated from in silico identified circRNA candidates (i.e., CircBase circRNAs) [1]. Several read simulators such as Mason [23], flux [24], ART [25], and CIRI-simulator [17] (used in this cited study) have been developed and successfully applied to generating reads. Since the authenticity of the synthetic datasets is totally unknown, the ambiguity/uncertainty of the simulated positive datasets should be minimized. We suggest that the simulated circRNA reads should not be generated from gene families, pseudogenes, mitochondrial/ribosomal genes, or unannotated exon boundaries (or junctions not having canonical splice site sequences). In addition, it would be better to evaluate accuracy of tools based on simulated data with a variety of data conditions of different read depths and of different read lengths, because accuracy of circRNA detection is often susceptible to these two factors [2,7]. The main flaw of these simulation-based approaches is that they are harmful to assess the effect of experimental artifacts on tools. For real RNA-seq datasets, although comparisons of NCL events detected with and without RNase R-treatment [13] (or oligo-dT selection [26]) have been used to demarcate cir-cRNAs versus non-circRNA NCL events, such approaches cannot distinguish between in vivo circRNAs and RT-based artifacts (e.g., template switching events). To control for experimental artifacts, an approach based on Drosophila hybrid mRNAs (D. melanogaster females vs. D. sechellia males) and a mixed mRNA-negative control sample was developed and successfully detected a considerable number of experimental artifacts [22]. However, it would be impossible to apply this approach to human studies. Recently, on the basis of the concept that RTasedependent RNA products are likely to be RT artifacts [2,11,12,27,28] and comparisons of different RTase products could act as effectively as RTase-free validation when detecting RTbased artifacts [11,12], we recently described an alternative pipeline to systematically extract potentially true positive circRNAs with controlling for experimental artifacts [29]. The RNAseq reads were generated from Avian Myeloblastosis Virus (AMV)-and Moloney Murine Leukemia Virus (MMLV)-derived samples, respectively. Only the non-poly(A)-selected RNA products, in which their NCL junctions were supported by both AMV-and MMLV-based reads, were regarded as true positive circRNAs [29]. Finally, wet-laboratory validated cir-cRNAs often serve as true positives and are employed for evaluating the sensitivity of tools. However, the collected circRNAs tested in various tissues or cell lines. It has been shown that most circRNAs were expressed in only a few tissues/cell types [17,[29][30][31][32][33]. As stated in the cited study [1], many of the collected circRNAs may not be expressed in the examined samples. In addition, we especially emphasize that the collected circRNAs should not include circRNAs in the results of the 11 tools (see also S2 Dataset). (E) Percentage of the ambiguous CircBase circRNA candidates of a tool when considering another circRNA predicting tool. The numbers in parentheses were ambiguous CircBase circRNA candidates detected by each tool.
https://doi.org/10.1371/journal.pcbi.1006158.g002 Example of unqualified RNAseq reads (i.e., ambiguous alignments with an alternative co-linear explanation or multiple hits). In this case (CircBase circRNA of COQ4, top), CIRI identified five supporting circRNA junction reads, four of which were unqualified reads with an alternative co-linear explanation. The BLAT-alignment result with an alternative co-linear explanation for Read validated by RT-PCR experiments with only one type of RTase because they may be RTdependent RNA products and derived from in vitro artifacts [12]. To minimize potential RTartifacts, the selected circRNAs should be confirmed by both RT-and non-RT-based experiments (e.g., Northern blot or RNase protection assay [34]) or at least, by different two types of RTase-based experiments [11,12,29].
With more and more circRNAs are detected, the reliability and function of most identified circRNAs remains an open question worthy of further investigation. To reduce the cost of further validation and functional analysis, it is necessary to carefully evaluate the reliability of the used circRNA-detecting tools. The abovementioned factors may considerably affect the evaluation results and partially explain the discrepancy observed in the sensitivity of the examined tools using synthetic (52-92%; Table 1 of the cited study [1]) and real (46-63% on HeLa cells; Fig. 5 of the cited study [1]) datasets. Therefore, we suggest that the accuracy of the 11 tools should be reevaluated (e.g., Table 1 and Fig. 4 of the cited study [1]), taking into account the factors discussed above when using the positive datasets.

Data access
The publicly available toolkit for downstream filtering of circRNA predictions with ambiguous alignment is downloadable at https://github.com/TreesLab/NCLtk. Supporting information S1 Dataset. The 3,150 uncertain CircBase circRNAs (summarized in Table 1) and the 3,580 CircBase circRNAs identified by all the 11 tools examined (used in Fig 3A).   Writing -original draft: Trees-Juen Chuang.
(simulate:514081) was shown (bottom). E, exon. (D) Comparison of CIRI-identified circRNA candidates before and after removing the unqualified reads. For the left panel, the numbers of CIRI-identified circRNA candidates were provided in parentheses. For the middle and right panels, the P values were evaluated using the two-tailed Wilcoxon rank-sum and Chisquare (for given probabilities) tests, respectively.