{'fraction_non_alphanumeric': 0.0293817832943575, 'fraction_numerical': 0.04089464531990166, 'mean_word_length': 5.001439366678661, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantifying full-length circular RNAs in cancer

Circular RNAs (circRNAs) are abundantly expressed in cancer. Their resistance to exonucleases enables them to have potentially stable interactions with different types of biomolecules. Alternative splicing can create different circRNA isoforms that have different sequences and unequal interaction potentials. The study of circRNA function thus requires knowledge of complete circRNA sequences. Here we describe psirc, a method that can identify full-length circRNA isoforms and quantify their expression levels from RNA sequencing data. We confirm the effectiveness and computational efficiency of psirc using both simulated and actual experimental data. Applying psirc on transcriptome profiles from nasopharyngeal carcinoma and normal nasopharynx samples, we discover and validate circRNA isoforms differentially expressed between the two groups. Compared with the assumed circular isoforms derived from linear transcript annotations, some of the alternatively spliced circular isoforms have 100 times higher expression and contain substantially fewer microRNA response elements, showing the importance of quantifying full-length circRNA isoforms.

Two types of biases need to be corrected, both by adjusting the effective lengths of 895 transcripts. The first type of bias is caused by a drop of sequencing depth at the ends of 896 linear transcripts, since fragments produced near the end of a transcript could be too short 897 to be sequenced by the experimental procedure. This type of bias does not apply to circular 898 transcripts, assuming that circular transcripts can be cut anywhere with equal chance. To 899 correct for this bias, the effective length of each linear transcript is defined as its actual 900 length minus the average fragment length (i.e., half the fragment length from each end).

901
The second type of bias is caused by the GC content of transcripts and the corresponding 902 non-uniform sequence amplification. To correct for this bias, the occurrence counts of all 903 6-mers on the sequencing reads are compared with their expected occurrences based on the 904 estimated abundance of each transcript isoform, as proposed previously (Bray et al., 2016). 905 To ensure only relevant sequencing reads are used to quantify the expression level of a 906 transcript, we additionally i) check the orientations and insert sizes of the pseudo-alignments 907 of read pairs, which for instance allow us to handle the reverse-overlap (RO) cases (Zheng 908 et al., 2019), ii) verify that reads that span BSJs are properly aligned on both sides, rather 909 than assuming it based on the positional information of pseudo-alignments alone, and iii) 910 consider the number of repeated sequences within a transcript that may lead to non-unique 911 pseudo-alignments. 912 Generation of simulated RNA-seq data 913 Sequencing reads of linear and circular transcripts were generated based on the RefSeq gene 914 annotation of the human genome reference GRCh38.p12. We first extracted from the an-915 notation all 23,851 genes with at least one transcript isoform containing at least one CDS 916 feature. We generated two sets of data based on these genes, namely an easier set with only 917 one linear transcript isoform and one circular transcript isoform per gene, and a second set 918 S4 with more complex gene structures. In both cases, transcript copy numbers were drawn in-919 dependently for each transcript isoform from the uniform distribution U (0, 1000). 9,000,000 920 fragments were generated for each group of genes (see below) according to these transcript 921 copy numbers with a normal length distribution N (250, 25), and these fragments were then 922 used to generate paired-end sequencing reads of 100bp long using Polyester-LC, which was 923 an extension of Polyester (Frazee et al., 2015)  transcript isoforms in total. Based on these linear transcript isoforms, we then generated a 935 linear transcript isoform and a circular transcript isoform of the same length with a certain 936 overlap ratio between them, where the overlap ratio ranged from 0% to 100%, forming 11 937 groups of data. Specifically, the synthetic linear transcript isoform occupies the beginning 938 of the real linear transcript isoform while the synthetic circular transcript isoform occupies 939 the end, with their lengths determined according to the target overlap ratio ( Figure S6). For 940 example, when the overlap ratio was 0%, the synthetic linear transcript isoform was equal to 941 the first half of the real linear transcript isoform and the synthetic circular transcript isoform 942 was equal to the second half. As a result, this data set contained 11 groups each with 1,000 943 synthetic transcript isoforms including 500 linear ones and 500 circular ones. Sequencing 944 reads were then generated from these 11 groups of transcript isoforms as explained above.

945
For the second set of data, we collected the first 10 genes on the TOR-sorted gene list 946 S5 each with at least 30 linear transcript isoforms. We then randomly sampled a certain number 947 of transcript isoforms from each gene, where the number ranged from 1 to 10, forming 10 948 groups. For each linear transcript isoform sampled, we also synthesized a corresponding 949 circular transcript isoform with the same sequence as the linear counterpart. As a result, we 950 had 20-200 transcript isoforms per group, and sequencing reads were generated from them 951 according to the procedure described above. We repeated the whole process 10 times with 952 10 sets of random transcripts sampled from the 10 genes.

953
Information contained in the output files of the NPC data analysis

954
The results of the NPC data analysis are summarized in four files, namely a file of all detected 955 BSJs, an annotation file of all candidate full-length linear and circular transcript isoforms, a 956 file of the corresponding sequences of these full-length isoforms, and a file that provides the 957 quantification results of these isoforms.

958
In the first file, provided in a tab-delimited format, each line lists a BSJ indexed by the 959 two genomic positions right next to the junction in the two respective exons. All candidate 960 full-length transcript isoforms that contain the BSJ are also listed, which serve as a mapping 961 between this file and the second file. The number of reads that support the junction in each 962 sample is listed next. The mean supporting read count (normalized by the total number 963 of reads in each sample) in each of the two sample groups is also shown, together with the 964 log 2 fold change of these mean values and the corresponding differential count p-value based 965 on Wilcoxon rank-sum test and its q-value. To focus on the most interesting cases, only 966 junctions with a non-zero read count in at least 70% NPC or NP samples are included.

967
In the second file, provided in the GTF format, each line lists an exon and the corre- w and x indicate the 3' and 5' exons that define the BSJ, and each pair of y i and z i defines 971 a junction that involve non-adjacent exons, and there can be zero, one or more of them. For 972 example, "HP1BP3 E9B1 2A4" means that this circular isoform involves a BSJ from exon 973 S6 9 of the HP1BP3 gene back to exon 1, and exon 2 is connected to exon 4 directly skipping 974 exon 3, followed by exons 4 to 9.

975
In the third file, provided in the FASTA format, the actual sequences of the full-length 976 transcript isoforms defined in the second file are listed.

977
In the fourth file, provided in a tab-delimited format, each line lists a candidate full-length 978 transcript isoform and its expression level in each of the samples in TPM values. The mean 979 expression in each of the two groups, log 2 fold change, differential expression p-value and  In the comparison between the quantification results of CircAST, isoCirc and psirc, we 988 used TPM values of CircAST and psirc, and normalized isoCirc's read counts by transcript 989 length before computing TPM values. Since both CircAST and psirc were designed for exonic 990 circRNAs, we considered only circRNAs formed purely by annotated exons. Since CircAST 991 used a strategy to filter out isoforms of low abundance, we also considered only circRNA 992 isoforms reported by psirc with a TPM of at least 2 −10 . 993 We also tried to run FUCHS (Metge et al., 2017) (https://github.com/dieterich-lab/ 994 FUCHS, git commit 82ce9602) but found that its outputs were inconsistent with the input gene 995 annotations in terms of exon lengths and exon boundaries, and therefore omitted it from the 996 comparisons.

S7
Full details of effective BSJ detection with low time and memory requirements 999 We verified the ability of psirc in identifying BSJs using three data sets. The first data 1000 set contained rRNA-depleted RNA-seq data from 11 human fetal samples (Table S1) consistency suggests that at least some of the BSJs identified by them but not CIRCexplorer2 1014 were legitimate, and thus they achieved higher sensitivity than CIRCexplorer2 on this data 1015 set.

1016
To evaluate the precision of the five methods, we tested them on two additional data sets, 1017 which contained matched rRNA-depleted and RNase R-treated RNA-seq data from human 1018 HeLa and Hs68 cells, respectively (Table S1). We applied each method to identify BSJs 1019 from the rRNA-depleted RNA-seq data, and then used the corresponding RNase R treated 1020 RNA-seq data to classify these junctions into Enriched, Unaffected, Depleted, and Abolished 1021 cases, defined as junctions with increased, similar, decreased but non-zero, and zero support 1022 in the RNase R data, respectively (Materials and Methods). we also ran it using only one thread, and found that the resulting elapsed time was still 1035 comparable to CircMarker and faster than CIRCExplorer2 even these methods were run 1036 with 25 threads, while the CPU time of psirc was further substantially reduced by 50-60% 1037 as compared to running it with 25 threads ( Figure S3).

1038
Together, these results show that in terms of identifying BSJs, psirc requires much less 1039 computational resources but still achieves comparable sensitivity and precision as the best 1040 of the other four methods.

Effects of read length 1042
To explore how read length may affect circRNA calling results, we have made use of in-house 1043 RNA-seq data of an NPC sample (NPC43). The data were produced using a protocol with 1044 rRNA depletion and the sequencing reads were 150bp long. Based on the original reads, 1045 we produced two additional data sets with read lengths of 100bp and 90bp, respectively, by 1046 artificially trimming the last bases of each read. The 100bp setting is for mimicking shorter 1047 reads that are still commonly produced nowadays, while the 90bp setting is for mimicking 1048 100bp reads with some bases trimmed due to low sequencing quality. 1049 We applied psirc to each of the three data sets to identify BSJs and full-length circRNAs.
From the results, read length does have some effects on the identification of BSJs and full-read lengths, there were also more than 2,000 that were only called when the reads were 1053 150bp long, suggesting that pseudo-alignments of some reads were different when they were 1054 trimmed to shorter lengths. There were also around 1,500 called only when the read length was 90bp or 100bp, which could be false positives due to incorrect pseudo-alignments. The 1056 same observations are also made for the full-length isoforms.

1057
In general, when sequencing reads are shorter, there is a higher chance of reads aligning to 1058 locations that they do not originate from due to non-unique sequences in the genome/transcriptome.

1059
If such an alignment appears to support a back-splicing junction, which would not have been 1060 supported with longer reads, a false positive back-splicing junction will be called.  Table S1: Short-read RNA-seq data sets used for comparing the performance of different methods. In the sample IDs of data sets 2 and 3, "R" refers to RNase R-treated data and "C" ("control", following terminology of previous studies) refers to untreated data, while "1" and "2" refer to two replicates.    Table S4: Full-length circRNA isoforms with a differential expression p-value ≤ 0.05 between the NPC and NP samples with at least 10 MREs of a single miRNA family. In the miRNA family name column, two miRNAs without a family ID in the TargetScan database are shown with the prefix "hsa-" or "ebv-", while all the others are shown without such a prefix. S13 BSJ Primers (F:forward; R:reverse) circAKT3 E9B1 F:5'-CAGATGCAGCCACCATGAAGA R:5'-CCCCCAACTTGGAGAAATGG circAKT3 E7B3 F:5'-CCTTGAAATATTCCTTCCAGACAAAAGACC R:5'-CCTCTGGAGTATCTACATGAAATGTTCTCTCT circRAPGEF5 E7B2 F:5'-GTGCAGAGAGAGCTAGCAGCT R:5'-CTGGGAAGGCGTATGTTCAAGG circRAPGEF5 E5B2 F:5'-GTACTGTGAATTTGAAAGAGAAGAAGAATGGC R:5'-CTGGGAAGGCGTATGTTCAAGG circROBO1 E2B2 F:5'-CATGGACGGAAAAGTAGACCTGATGAAG R:5'-GCTTTGCAGTTCAAAGTTGCAGGT circROBO1 E6B4 F:5'-GAAAATTAGGAAGGTGACAGCTGGTGAC R:5'-CTTCACTCTCACGTTCCCCAACC circNETO2 E7B2 F:5'-GATGTAATGCTTAAAACAGGAATTGGAGTGAT R:5'-GTTTTTTGGGCCACGGCAATCC circVEGFC E4B2 F:5'-CCAATCACACTTCCTGCCGATG R:5'-CACTTGTACATTTTCCAATATTCTGGGTAGAGTACA circITCH E20B19 F:5'-TGAGTCCAACTGGTTTGTTCAAGATACG R:5'-ACTTCTCCGTTGACAAAGAAATTCTAGGTG circNTRK2 E12B10 F:5'-GTCTATGCTGTGGTGGTGATTGCG R:5'-CTTCATAAATTACATCAGGATAATTTGGGTTTGC circNTRK2 E12B3 F:5'-GTCTATGCTGTGGTGGTGATTGCG R:5'-AGATTTCTCAGTCCCACATAAGCTTCAAC circNTRK2 E6B3 F:5'-TCCAGTCCAGACACTCAGGATTTGTAC R:5'-AGATTTCTCAGTCCCACATAAGCTTCAAC circRPMS1 E6B4 F:5'-CTAGTGCTGCATGGGCTCCT R:5'-GTCATACGCCCGTATTCACACATTC  Transcript isoform Primers (F:forward; R:reverse; P:probe) circNTRK2 E12B10 F:5'-CGGTCGGGAACATCTCTCG R:5'-GTCATTCGCTGCAGTTCCATAAT P:5'-TTTGGCATGAAAGGTGCAAACCCAAA circRPMS1 E6B4 F:5'-CTAGTGCTGCATGGGCTCCT R:5'-GTCATACGCCCGTATTCACA P:5'-GTCGACGGGCAAGGTCCGGCGTGTCC Table S7: Primers used in the RT-qPCR experiments for verifying the differentially expressed full-length circular transcript isoforms between the NPC and NP groups.   CT value Figure S5: Comparison of different pipelines based on their deduced normalized read counts that support each BSJ from the human fetal samples. Explanations are the same as in Figure 2. For BSJs not detected by the BSJ-detection method of a pipeline, their read supports are shown as zero (plus a small constant added to avoid taking log of zeroes). CircMiner failed to complete its run on GSM1567915 in five days and so BSJs from this sample are excluded from its plots.  Figure S6: A schematic figure illustrating the procedure for generating the first set of simulated data. The black line on the top indicates a reference sequence used as the template for defining linear and circular transcript isoforms. Eleven pairs of linear and circular transcript isoforms were then defined based on the reference sequence, with their sequence overlap ranging from 0% to 100%, which correspond to cases of increasing difficulty. In total, 500 reference sequences were used to define these simulated linear and circular transcript isoforms, and sequencing reads were generated based on these isoform definitions.    Figure S9: The workflow for quantifying the expression levels of linear and circular transcript isoforms. Starting from a set of linear and circular transcript isoforms (A), the indexing step (B) constructs a T-DBG based on the order of k-mers in the isoforms, with nodes in each non-branching path condensed into a single node, and each such node given a set of colors that represent the corresponding isoforms supporting it, called its k-compatibility class. The pseudo-alignment step (C) then pseudo-aligns each sequence read to the isoforms from which it may be originated using the T-DBG. Finally, the quantification step (D) uses the pseudoalignment results to quantify the expression level of each linear and circular isoform.  Figure S10: Validation of the BSJs in the RT-PCR products of the selected BSJ and fulllength circRNA candidates by Sanger sequencing.

Method
File S1: BSJs identified from the NPC and normal NP samples. Each row corresponds 1065 to a junction. The different columns provide information about the potential full-length 1066 isoforms that involve the current junction, the number of reads that support the junction 1067 in each sample, the average number of reads that support the junction in each of the two 1068 groups, and the corresponding log 2 fold change, p-value and FDR q-value for the differential 1069 supporting read counts.