Is this really the end? Template-switching artifacts resemble alternative polyadenylation

Background: Alternative polyadenylation is commonly examined using cDNA sequencing, which is known to be affected by template-switching artifacts. However, the effects of such template-switching artifacts on alternative polyadenylation are generally disregarded, while alternative polyadenylation artifacts are attributed to internal priming. Results: Here, we analyzed both long-read cDNA sequencing and direct RNA sequencing data of two organisms, generated by different sequencing platforms. We developed a filtering algorithm which takes into consideration that template-switching can be a source of artifactual polyadenylation when filtering out spurious polyadenylation sites. The algorithm outperformed the conventional internal priming filters based on comparison to direct RNA sequencing data. We also showed that the polyadenylation artifacts arise in cDNA sequencing at consecutive stretches of as few as three adenines. There was no substantial difference between the lengths of poly(A) tails at the artifactual and the true transcriptional end sites even though it is expected that internal priming artifacts have shorter poly(A) tails than genuine polyadenylated reads. Conclusions: Our findings suggest that template switching plays an important role in the generation of spurious polyadenylation and support the need for more rigorous filtering of artifactual polyadenylation sites in cDNA data, or that alternative polyadenylation should be annotated using native RNA sequencing.

filtering algorithm which takes into consideration that template-switching can be a source of artifactual polyadenylation when filtering out spurious polyadenylation sites. The algorithm outperformed the conventional internal priming filters based on comparison to direct RNA sequencing data. We also showed that the polyadenylation artifacts arise in cDNA sequencing at consecutive stretches of as few as three adenines. There was no substantial difference between the lengths of poly(A) tails at the artifactual and the true transcriptional end sites even though it is expected that internal priming artifacts have shorter poly(A) tails than genuine polyadenylated reads. Conclusions: Our findings suggest that template switching plays an important role in the generation of spurious polyadenylation and support the need for more rigorous filtering of artifactual polyadenylation sites in cDNA data, or that alternative polyadenylation should be annotated using native RNA sequencing.

Background
The majority of human genes utilize alternative polyadenylation (APA) sites [1,2]. APA is known to substantially influence gene expression [3,4] and plays a role in disease development [5]. cDNA sequencing greatly facilitates the analysis of APA [6]; however, it is influenced by internal-priming artifacts. In the case of internal priming, the oligod(T) primer attaches to an adenine A-rich region of the transcript and initiates transcription from this region rather than the poly(A) tail [7] (Figure 1a). RNA ligation can be applied to enable specific amplification of the 3'-ends of transcripts and to negate the effects of internal priming [8,9]. Regular poly(A)-seq data generated using oligod(T) primers are usually filtered so that poly(A) (pA) sites in A-rich genomic regions are discarded. A-rich regions are often defined as stretches of 6 or more consecutive As or 20-nt-long windows comprising more than 60% adenines [10][11][12][13][14]. In a recent long-read cDNA sequencing study of the human cytomegalovirus transcriptome, we described potentially artifactual pA sites arising from homopolymer stretches-sometimes as short as only three As [15]. Based on this finding, we propose that such artifacts are produced by template switching (TS).
TS refers to the ability of DNA polymerase to discontinue elongation while still binding the newly synthesized strand and to reinitiate synthesis at a homologous locus of another nucleic acid strand. [16] (Figure 1b). Both DNA-dependent and RNA-dependent DNA polymerases reportedly participate in TS [17,18]. This phenomenon has been shown to occur more frequently if the concentration of the templates is high, the homologous sequences are long, or the RT temperature is low [19,20]. Polymerase pausing may also facilitate TS [20]. Another study found that direct repeats of three to six nucleotides can trigger TS; however, longer homologous sequences (i.e., 12-24 nt) resulted substantially in more artifacts [21].

Results
We hypothesized that artifactual polyadenylation events at shorter stretches of adenines are caused by TS. To characterize artifactual and genuine transcriptional end sites, we analyzed publicly available datasets in which both direct (d)RNA and cDNA sequencing data were available for the same cell lines. Potential pA sites were determined based on the cDNA sequencing data and then compared with direct (d)RNA sequencing data to 4 identify artifacts in the cDNA sequencing results. In total, 87,980 and 403 potential pA sites were identified in the human cDNA dataset [22] and human cytomegalovirus (HCMV) cDNA dataset, respectively [23]. Figure 2a (see alsoSupplementary Figure 1a)shows that the more As located upstream of a pA site, the less likely it was to be confirmed by dRNA sequencing. A decrease in the ratio of confirmed pA sites was already observed at relatively low numbers of As. It should be noted that not all potential pA sites missing from the dRNA sequencing data are artifacts; it is also possible that a pA site was not supported by dRNA reads simply because dRNA sequencing had lower coverage in that region. However, the adenine content of a region is not expected to reduce dRNA sequencing coverage. Therefore, decreased dRNA support for potential poly(A) sites in Arich regions points to an increase in the number of artifacts. The list of potential human poly(A) sites was also compared with PolyA_DB, a database of poly(A) sites validated by the 3'READS+ method, which uses RNase H digestion and RNA ligation to prevent internal priming [24].
Based on the hypothesis that TS-not internal priming-produces artifacts at shorter stretches of As, we devised an algorithm that differentiates between artifacts and genuine transcriptional end sites (TES). The algorithm considered the number of As in the region immediately upstream of a pA site, the number of reads in the proximity of the pA site that fell outside of A-rich regions, and the ratio of polyadenylated reads to the coverage of the region (Supplementary Figure 2). The HCMV dataset comprised of four different cDNA sequencing experiments. The potential pA sites called differently in different experiments (e.g. called as artifact according to the Sequel data but called as TES in the MinION sequencing data) were regarded as TESs. Of the 859 calls, 24 were discordant out of which 18 were confirmed by dRNA sequencing to be genuine TESs, which suggests that the algorithm is more likely to generate false negative than false positive results. The 5 algorithm proved to be slightly more specific and substantially more sensitive than the conventional filtering method (≥6 consecutive or ≥12 As in a 20-nt region) (Supplementary Figure 1b). Specificity was increased by excluding potential pA sites with  Figure 1c). Nevertheless, many putative artifacts were detected in regions with as few as 3-5 As, and the likelihood of these artifactual pA sites to be detected by dRNA sequencing did not decrease when these sites contained more As in a 20-nt window (Supplementary Table 1). The sites confirmed by PolyA_DB data were the most likely to have been confirmed by dRNA sequencing data, although many sites not in PolyA_DB were also detected by dRNA sequencing.
Interestingly, potential pA sites in PolyA_DB were less likely to have been confirmed by dRNA sequencing if they were in A-rich regions, although this phenomenon was not as prominent as that for other pA sites (Figure 2d and Supplementary Figure 1d). All filtering methods performed worse at regions containing more than ten adenines.
The putative TESs identified by the algorithm differed greatly from the putative TS artifacts ( Figure 3). The nucleotide composition surrounding TESs showed specific motifs commonly observed around cleavage sites ( Figure 3a). Putative TESs were often preceded by common polyadenylation signals (PAS), whereas putative TS artifacts generally lacked such signals (Figure 3b). PAS usage in HCMV, like in other herpesviruses [25], is very similar to its host. Accordingly, the PAS usage of HCMV TESs was very similar to that of human TESs, but different from putative artifacts (Figure 3b). In cases where putative artifacts were preceded by PASs, the signal was often not at the expected distance of 25 nt, as observed at putative TESs ( Figure 3c). Polyadenylation at a given pA site does not always occur at the same nucleotide; rather, it may occur at any of several nucleotides around the most frequently cleaved nucleotide [15,26,27]. This phenomenon was observed at putative TESs in both the human and HCMV datasets but absent at artifactual pA sites (Figure 3d). The accumulation of many artifactual reads at certain positions is due to an erroneous alignment to homopolymer As, whereas the genuine cleavage sites are more spread around a given position. Figure

Discussion
We analyzed poly(A) + cDNA sequencing data of two species (human and HCMV), stemming from three different long-read sequencing platforms (RSII, Sequel and MinION), generated 7 by three different library preparation methods (Iso-Seq, Cap and poly(A)-selection, as well as only poly(A)-selection), and then compared them to dRNA sequencing data obtained by the MinION platform. Our analyses confirmed that artifacts arising in A-rich regions complicate the study of alternative polyadenylation. This phenomenon is generally accredited to internal priming [7]. Given our findings, we argue that TS is more likely to be responsible for these artifacts as many artifacts were detected in regions with rather few As (sometimes three to five), which make oligod(T) primer binding unlikely. Further, it would be expected that reads ending at artifactual sites produced by internal priming would not contain poly(A) tails substantially longer than the oligod(T) primer. However, we found that poly(A) tails at artifactual sites were longer than the primer and not shorter than poly(A) tails at bona fide TESs. We thus developed a filtering algorithm to differentiate TS artifacts from genuine TESs. Based on comparison with dRNA sequencing data, the filtering algorithm performed better than conventional internal priming filters but was less specific than using a curated list of pA sites. We suggest that, although internal priming is likely to contribute to the number of artifacts in very A-rich regions, artifacts in regions with lower adenine content are generated by TS. Even though it was not part of the filtering criteria, sites that the algorithm classified as TESs were likely to contain consensus polyadenylation motifs. This result indicates that machine learning algorithms can distinguish even more specifically between TESs and artifacts using more sequence information. However, a large training dataset would be required for machine learning to be efficient, and such datasets are not available at the time.
Our findings were obtained using long-read sequencing datasets. While it may seem sensible to extend our conclusions to short-read sequencing data and other results obtained by cDNA sequencing, it must be noted that some aspects of long-read sequencing promote the production of template-switching artifacts. Firstly, long-read sequencing usually necessitates reverse transcription of the whole transcript, not only its most 3' fragment, which is an option for short-read sequencing. Reverse-transcribing more genomic regions provides more potential templates for TS. Secondly, SMART technology, which is widely used in long-read sequencing studies to produce full-length transcripts, requires ideal conditions for TS [28,29]. Whereas the SMART protocol allows reverse transcription to be carried out at 50°C, the second strand synthesis in the same reaction mixture must occur at 42°C to allow strand switching. The characteristics of long-read sequencing library preparation increase the impact of TS; nevertheless, similar artifacts could influence other reverse-transcription-based methods as well.

Conclusions
TS is known to produce cDNA artifacts, however its effects on the analysis alternative polyadenylation have never been discussed until now. Considering that the poly(A) tail is likely the most frequent template in most transcriptomic libraries, polyadenylation artifacts may be even more prevalent than the more reviewed splicing artifacts. The effects of TS on short-read sequencing can be mitigated by higher reverse-transcription temperatures [30] or by employing high read-count thresholds that are easier to implement due to the higher throughput of these sequencing methods. Long-read cDNA sequencing approaches are currently more prone to TS artifacts, but these artifacts can be ruled out by dRNA sequencing or curated pA-site databases when available for the study organism. If such datasets are unavailable or inappropriate, we advise strict filtering that also considers the effects of TS. The filtering method presented here can be applied to data from any long-read sequencing platform and performs better than the conventional filtering method. An important advantage of this filtering method is its higher sensitivity that allows utilization of more data, which is crucial for long-read RNA sequencing as it has a lower throughput than short-read sequencing methods [31]. 9

Data acquisition
Two long-read cDNA and dRNA sequencing datasets were downloaded and analyzed during the study. (1) The human cDNA and dRNA sequencing FASTQ reads of the Nanopore WGS Consortium (https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md) were generated by extracting RNA from the GM12878 human cell line (Ceph/Utah pedigree) and sequenced on MinION flow cells (FLO-MIN106) using R9.4 chemistry (SQK-RNA001 and SQK-LSK108 kits) [22]. This dataset will be referred to as the human dataset.

Mapping and read processing
The computational pipeline of the study is summarized in Supplementary Figure 2. The processing steps for the human and HCMV data were the same. The reads were mapped using minimap2 [33] to the human genome (hg19) and to the HCMV strain Towne varS genome (LT907985). Reads from the HCMV dataset were only mapped to the viral genome, reads from the viral infected host were not used. The mapper settings were "-ax splice -Y -C5" for the cDNA and "-ax splice -uf -k14" for the dRNA sequencing reads. Coverage and dRNA read endings were determined using bedtools [34]. As dRNA sequencing does not accurately sequence the terminal poly(A) tail of the reads, every dRNA read ending was counted. A genomic locus was confirmed as a poly(A) site confirmed by dRNA sequencing, if at least 0.5% of the overlapping dRNA reads ended in the 21-nt window (10nt upstream + the locus + 10nt downstream = 21nt) around the locus.

Identifying potential poly(A) sites in the cDNA sequencing data
The LoRTIA toolkit (https://github.com/zsolt-balazs/LoRTIA) was used to identify potential poly(A) sites in the cDNA sequencing data. A genomic locus was considered a potential poly(A) site when at least two reads and at least 0.1% of the overlapping reads ended at a given nucleotide. In a 21-nt window, the genomic position with the highest number of poly(A)+ reads was selected as the potential poly(A) site. The separate experiments of the HCMV dataset were analyzed separately and the results were joined to create the list of potential HCMV poly(A) sites. The LoRTIA toolkit was also used to mark reads which ended in A-rich genomic regions (three or more consecutive As as potentially artefactual reads).
When characterizing high confidence calls only the sites where more than ten reads ended in a 21-nt window around the locus were analyzed.

Defining A-rich regions
We have deployed a slightly different definition of A-rich regions than it is commonly used in the literature. A common approach is to count the number of consecutive As in a region surrounding the poly(A) site. Another method is to count the number As in a given, often 20-nt-long window (because 20nt is the primer length). Instead, we iterated the 20nt upstream of a potential poly(A) site and incremented a counter each time an A was iterated, all the other nucleotides were counted as -1. If the counter reached -1, the iteration was halted and the highest count was regarded as the A-count of the region (Supplementary Figure 2).

Filtering out template-switching (TS) artefacts
The potential poly(A) sites which were not at A-rich loci, were accepted by our script as transcriptional end sites (TES). The potential poly(A) sites at A-rich loci were accepted as TES if the number of reads in a 21-nt window around that loci contained either more reads which ended in a non-A-rich region than reads which ended in an A-rich region or a proportion of overlapping reads greater than 0.81+2-100(120-n -0.08), where n is the number of As in the A-rich region. The potential poly(A) sites which did not meet these requirements were classified as TS artefacts. The genomic loci which contained at least 100 times more 5'-ends from the opposite strand than 3'-ends were also classified as artefactual poly(A) sites. The sites which were classified as TS artefacts, but were no further than 50nt away from sites that were classified as TES, were excluded from comparisons and sequence analyses in order to prevent the characteristics of the TES to interfere with the characteristics of TS artefacts (e.g. a polyadenylation signal that is upstream of a TES could also be upstream of a TS artefact, similarly, a direct RNA read ending in the proximity of a TES could also fall into the 21-nt window of an artefactual site).

Detecting polyadenylation signals (PAS)
The 12 most commonly utilized human PAS [26] were searched for in the region 40nt upstream of each potential poly(A) site. Only exact matches were accepted. If the sequence upstream of a potential poly(A) site contained multiple PAS sequences, the most common PAS was assigned to the site.

Measuring poly(A) tail length
The length of the sequenced poly(A) tail of the cDNA reads was measured for each read that contained a poly(A) tail recognized by the LoRTIA toolkit. The sequence containing the 3'-terminal 30nt of the mapped part of the read and the first 150nt of the 3'-terminal soft clipped (unmapped) part of the read was aligned to a stretch of 180 As using the Smith-Waterman algorithm (match score: +2, mismatch score: -3, gap opening score: -3,

Availability of supporting data and materials
The human dataset used in this study was downloaded from GitHub site of the WGS consortium (https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md).

Supplementary Information
Supplementary Figure 1 One criterion was that the number of reads ending in not A-rich regions (less than three As) in a 21-nt window, had to be higher than the number of reads ending in A-rich regions (three or more As). The second criterion was that the number of supporting reads had to be higher than the percentage of overlapping reads calculated by a logistic function that takes into consideration the adenine content of the region. A potential poly(A) site was classified as a TES if any of the two criteria was fulfilled.