BACKGROUND
recently, there has been an increased interest in a novel class of transcripts where the exon order differs from that found in the genome .
once considered cloning artefacts  <cit>  or products of aberrant splicing  <cit> , it is now established that the majority of these molecules represent circular rna species  , although some linear transcripts have been reported  <cit> .
thousands of these novel transcripts have now been identified in a variety of eukaryotic cells  <cit> , many are conserved across species  <cit> , suggesting functional relevance, and two  have been shown to harbour numerous mirna binding sites and act as mirna sponges to modulate gene expression  <cit> .
recent reports also implicate circrnas in synaptic development  <cit>  and some have expression patterns that correlate with diseases  and may act as biomarkers for ageing  <cit> .
despite these reports, the function of the vast majority of these transcripts has not been established.
the defining feature of these transcripts at the sequence level is the presence of a splice junction with exons in an order inconsistent with their position in the genome.
as this feature alone does not enable inference of structure or mechanistic origins, we use the term post-transcriptional exon shuffled  transcripts to collectively describe this population of rna molecules  <cit> .
recent reports have shown that the vast majority of these transcripts emanate from known genes  <cit> , utilise known splice junctions, and that their biogenesis competes with splicing of canonical transcripts  <cit> .
transcripts arising from ptes specifically exclude chimeric rnas without known splice junctions, and a class of circular rnas comprised of spliced introns, cirnas  <cit> .
many computational methods for identifying chimeric rna molecules from high-throughput rna sequence data have been described.
the majority of these are designed to detect heterotypic trans-splicing and fused genes, so are not suitable for ptes detection , or require post processing steps to analyse reads supporting different types of splicing events  <cit> .
recently, however, a number of programs for ptes discovery have been described and used primarily for circrna characterisation .
most analyse reads which fail to fully align to a reference sequence, and split these into two or more segments which are then independently aligned to define rearrangements.
some make use of gene annotation to guide discovery  <cit> , while others adopt an unbiased genome-wide approach to capture structures which do not utilise known splice junctions or are non-genic  <cit> .
in addition, the occurrence of ptes can be inferred when two paired end reads map to the transcriptome in a configuration that is not consistent with a linear transcript  <cit> .
the identification of ptes exon junctions within rnaseq data is, however, confounded by known artefacts.
false positives can arise from template switching during cdna synthesis  <cit> , from genes with duplicated exons  <cit> , from transcription read-through between genes in close proximity due to weak termination signals  <cit> , and from closely related genes within duplicons or tandem arrays  <cit> .
although experimental enrichment has been combined with informatic approaches to define bona-fide circrnas , many classes of false positive structures are not directly excluded by existing identification methods.
for example, reads defining  <dig> of the  <dig> most abundant human circrnas reported by memczak et al.
<cit>  map with high sequence identity to the reference sequence and include  <dig> which are indistinguishable from linear refseq entries .
furthermore, a recent experimental analysis of previously identified ptes trancripts concluded that many are template switching artefacts  <cit> , and template switching predominantly leads to rearragements where the breaks do not occur at splice junctions  <cit> .fig.
<dig> examples of intragenic false positives.
schematic diagrams showing inferred structure and key sequence relationships for  <dig> of the  <dig> most abundant circrnas reported in  <cit> .
in each case, the inferred structure shares 100 % identity to a linear transcript spanning the defining exon-exon junction.
within the top  <dig>  hsa_circ_ <dig>   <dig> and  <dig> show similar patterns of identity to multiple genomic locations.
blue – inferred donor exon, red – inferred acceptor exon, black – upstream or downstream refseq exon sharing 100 % identity to donor or acceptor exon at junction.
approximate chromosomal locations  are shown
here we present a method, ptesfinder, that identifies putative ptes structures by mapping rnaseq reads to sequence models generated using existing transcript annotation.
it then applies a series of mapping and alignment filters to systematically remove known classes of false positives.
it does not make use of paired end  mapping information as the lack of intervening sequence precludes such filtering and may affect specificity.
we first describe the implementation of this method, and then investigate the effects of different filtering criteria.
the program requires certain user adjustable parameters; we therefore also explore systematically the choice of these parameters.
finally, through analysis of real and simulated data, we compare ptesfinder to other publicly available methods  <cit>  which have been used to identify circrna transcripts in both cell lines and tissues.
implementation
pipeline for ptes discovery
ptesfinder requires as input files: rnaseq data in fastq format  <cit> , genomic reference in fasta format, and an annotated transcriptome reference in bed format  <cit> .
the pipeline is split into three phases : a discovery phase to identify putative ptes structures within rnaseq data and define ptes transcript models, an evaluation phase to assess these ptes models, and a filtering phase to exclude potential false positives.fig.
<dig> overview of ptes discovery pipeline.
the workflow includes three major phases: discovery phase, evaluation phase and filtering phase.
putative ptes structures discovered using 20 bp anchor reads are evaluated by aligning full fastq reads to the models.
the filtering phase includes stringent criteria designed to systematically exclude all known classes of false positive structures
discovery phase
short sequences from both ends of each fastq read  are aligned to the transcriptome reference using bowtie  <cit>  with tolerance for a single mismatch.
pairs of anchors from the same read that map to the same gene and in the same orientation, but which map in inverted order with respect to their order in the sequencing read, are then identified.
this excludes reads emanating from fused genes and sense-antisense template switching artefacts.
retained anchor pairs are then used to determine the exon junctions which define putative ptes events and create sequence models  of the inferred products.
these constructs are generated by concatenating the last 65 bp of the 5′ exon and the first 65 bp of the 3′ exon.
the segment size of  <dig> is used by default, with the full exon sequence used if an exon is smaller than 65 bp.
this parameter is adjustable to accommodate various rnaseq read lengths, and we recommend that it be at least 10 bp shorter than the read length to ensure that only reads mapping across ptes defining junctions are processed in subsequent filtering steps .
evaluation phase
all the original reads are then re-mapped to ptes constructs generated in the discovery phase.
this serves two purposes.
first, as rnaseq reads can be short, this enables reads with putative ptes exon junctions close to their termini to be accurately mapped.
second, it enables read mapping scores obtained using the ptes constructs to be directly compared to scores obtained from genomic and transcriptomic alignments for filtering purposes .
optionally, evaluation can also be ‘guided’ by supplying constructs of previously discovered ptes structures, effectively bypassing the discovery phase.
filtering phase
to eliminate potential false positives originating from the genome under investigation, all the original reads are mapped to both genomic and transcriptomic references.
the number of edits required for alignment , and the number of perfectly aligned base pairs, are used to remove reads which align as well or better to either of these reference sequences than to the ptes constructs.
to reduce template switching artefacts, which have heterogeneous junction points within short regions of often imperfect sequence homology  <cit> , reads which do not align perfectly to the exon junctions which define ptes are also removed using junctional filters.
first, a user adjustable minimum junction span  parameter is applied to ensure that there are no mismatches or indels within ‘n’ nucleotides either side of the junction position, where n is an even integer.
second, to eliminate reads with regions of low quality alignment, a user adjustable segment percent identity  parameter is also applied independently to the segments on either side of the ptes junction, such that for a read to be retained both must meet or exceed the specified pid when aligned to the ptes construct.
these user adjustable filters rely on alignment summaries provided by the nm field, md field and cigar in the sam files  <cit> .
the output includes the coordinates of the exon end involved in the junctions, a descriptor of the ptes  and the number of reads supporting the structure.
this is presented in bed format  <cit> .
a second file contains additional information, read counts of all canonical exon junctions from transcripts where a ptes structure has been identified, to facilitate comparison with ptes counts.
assessment of pipeline and comparisons to other methods
rnaseq data from jeck et al.
<cit>  were analysed at various combinations of jspan and pid .
all analyses were performed with and without genomic and transcriptomic filters applied to enable reads discarded by each filtering criterion to be identified.
the numbers of ptes structures identified and supporting reads were also recorded.
to assess sensitivity and specificity, simulated datasets were generated using all published ptes structures within circbase.org  <cit> .
for each dataset,  <dig> ptes junctions were randomly selected along with  <dig> canonical junctions, and constructs were generated for each junction by concatenating the full sequence of both exons involved in each case.
100 bp simulated reads with random start positions within each construct were then generated.
scripts published by memczak et al.
<cit>  , ciri v.  <dig>   <cit>  , circrna_finder  <cit>  , and mapsplice v.  <dig> .
<dig>  <cit>  used in  <cit>  , were compared to ptesfinder by analysing leukocytes cell line rnaseq data , fibroblasts rnaseq data , and simulated data.
for each simulation, transcripts correctly identified by each method were determined by comparing genomic coordinates of identified transcripts with the genomic coordinates of transcripts expected to be recovered from within each dataset.
the numbers of correctly identified ptes transcripts , incorrectly identified ptes transcripts , ptes transcripts incorrectly excluded , and canonical junctions correctly excluded , were used to estimate sensitivity: tp / , specificity: tn / , and false discovery rate : fp / .
RESULTS
ptesfinder uses established rnaseq tools  to identify putative ptes structures, and then systematically excludes known classes of false positive structures by applying genomic, transcriptomic and junctional  filters .
as an initial assessment of ptesfinder function, rnaseq data from human fibroblast total rna which has previously been mined for circrnas , were analysed both with and without the application of the genomic and transcriptomic alignment filters.
reads recovered during analysis, together with alignment edit distances of reads identified by each filter applied seperately, are shown in fig.
3a.fig.
<dig> summary of reads excluded by filters.
a edit distance distribution of reads filtered out by genomic, transcriptomic and junctional  filters.
inset: seven structures are supported by  <dig> reads with 100 % alignment but are excluded by the transcriptomic filter.
b 30 % of reads filtered out support a false positive structure from  <dig>  s rrna and are excluded in this plot.
c venn diagram showing number of reads excluded by filters.
majority of false positive reads are excluded by all three filters.
each filter also excludes a distinct population of false positive reads
filters target overlapping populations of reads but none is redundant
from a total of over  <dig> million reads analysed, approximately  <dig>  %  have shuffled coordinates with respect to exon position .
however, of these only  <dig>  map to ptes sequence constructs generated during the discovery phase, indicating that most of the reads with rearranged anchor pairs do not map to single genes and/or known exon junctions.
approximately 85 %  of the reads which map to ptes constructs are subsequently removed by the genomic, transcriptomic and junctional  filters, with the majority being identified by more than one filter.
for instance, over 98 % of reads excluded by the genomic filter are also excluded by the transcriptomic filter, and 60 %  of all filtered reads are identified by all three.
most of these have high edit distances  indicative of low quality alignment.
despite this, the genomic, transcriptomic and junctional filters  uniquely exclude ~ <dig>  % , ~ <dig>  %%  and  <dig>  %  of reads mapping to ptes models respectively , indicating that none is wholly redundant.
the subset of reads identified specifically by the junctional and transcriptomic filters are defined by low edit distances of between  <dig> and  <dig> , although a small number of reads excluded by the transcriptome filter  map perfectly to putative ptes constructs with nm =  <dig> .
fig.
3a also reveals a bimodal distribution of mapping qualities for reads excluded by all three filters with peaks at nm =  <dig> and nm =  <dig>  upon manual analysis, most of the excluded reads with nm =  <dig> were found to support a false positive structure from  <dig>  s rrna .
comparable rrna derived structures have been identified previously and filtered manually  <cit> .
in fig.
3b, reads supporting this structure have been removed to show the underlying distribution of mapping quality scores.
reads excluded by specific filters have different origins
to investigate the activity of specific filters further, the mapping co-ordinates of reads removed by the genomic filter were first compared to the co-ordinates of annotated pseudogenes and segmental duplications.
this established that ~74 % of reads excluded by the genomic filter had superior alignments to segmental duplications, and ~12 % had superior alignments to pseudogenes.
the  <dig> reads identified by the genomic filter but not by the transcriptomic filter were also found to be enriched for reads derived from segmental duplications and pseudogenes .
we next used blat  <cit>  to manually investigate the  <dig> reads excluded specifically by the transcriptome filter which mapped perfectly to putative ptes constructs .
these support  <dig> putative ptes structures from  <dig> genes .
however, blat analysis established that they all also mapped contiguously with ~100 % identity to the transcriptome due to high sequence identity between neighbouring exons.
for example,  <dig> reads which support a putative single exon ptes structure  map with ~100 % identity to exons  <dig> and  <dig> of the canonical hnrnph <dig> transcript  due to high sequence identity between these neighbouring exons.
as a result, these reads cannot be taken as supporting evidence for ptes.
it is noteworthy that such structures will pass any qualitative filter criterion requiring only unambiguous mapping to ptes constructs, illustrating the value of the transcriptome filter.
finally, manual analysis of a subset of the  <dig> reads identified only by the junctional filters established that these support structures with distinct patterns of suboptimal mapping, such as low alignment quality specific to only one of the two exons in the structure , and low sequence identity specifically at the junction , the latter being consistent with the expected pattern of alignment for template switching artefacts  <cit> .
as one further assessment of the filters, we analysed rnaseq data derived from fibroblast rna which had been pre-digested with rnase r. this selectively removes linear rnas, and enriches for circrnas  <cit> , and has been shown to significanty increase the recovery of ptes reads.
however, we would anticipate that this would also selectively remove false positives derived from pseudogenes and segmental duplications which mimic ptes structures, without necessarily reducing reverse transcription artefacts such as template switching.
only ~12 % of reads from the rnaser digested sample which map to ptes sequence constructs are identified by the genomic and transcriptomic filters , compared to 69 % in the undigested sample.
furthermore, only 17 % of these map to segmental duplications, compared to 74 % in the undigested sample.
in contrast, the proportion of reads excluded by the junctional filters is considerably higher after rnaser digestion, consistent with expectation.
pid has greater impact than jspan
to investigate the impact of varying the user defined jspan and pid parameters which comprise the junctional filter, the same data was re-analysed using  <dig> different combinations of these parameters, both with and without the genomic and transcriptome filters applied .
this established that varying the pid has a greater impact than varying the jspan, with  <dig> reads filtered with maximal pid  and lowest jspan  compared to only  <dig> reads filtered with the maximal jspan  at lowest pid .
furthermore, varying the pid between 60 % and 75 % has little impact at any jspan value, but above 75 % there is a linear relationship with the number of reads filtered.
as the default junctional filter parameters failed to identify some reads excluded by the other filters , this analysis was repeated using only these reads to establish the jspan and pid parameters required to identify them.
over 99 % of these reads are excluded with the most stringent junctional filter parameters .
furthermore, the vast majority are filtered with a pid of 85 %, suggesting this is a logical setting for this parameter.
the jspan setting only has a major impact at low pids .fig.
<dig> effect of varying junctional filter parameters.
a number of reads passing filter at different combinations of jspan and pid.
b percentage of reads only excluded by transcriptome and genomic filters at default settings, which are filtered at different combinations of jspan and pid
specificity, sensitivity and comparison with other methods
to assess the sensitivity and specificity of the pipeline and compare it to other methods, simulated reads were generated from previously identified ptes and associated canonical transcripts, and analysed at various read depths of coverage using default parameters.
in addition to assessing ptesfinder for de novo ptes discovery, the use of constructs of previously reported structures for guided discovery was also assessed , as were four publicly available methods which have previously been employed to identify circrna transcripts: mapsplice v  <dig> .
<dig>  <cit>  used in  <cit> , ciri v.  <dig>   <cit> , circrna_finder  <cit>  and the method used by memczak et al.
<cit> .
results from  <dig> simulated datasets are presented in fig.
5a-c, and illustrate that sensitivity varies considerably with coverage, and between methods.
at read coverage of  <dig>  the sensitivity of ptesfinder is below  <dig> .
this can be attributed to ptes junctions occurring within the terminal 20 bp of reads, as the low tolerance for mismatches during anchor mapping will result in their elimination.
however, sensitivity reaches >90 % at coverage of  <dig> or higher for both guided and unguided analyses, with guided ptesfinder being equally or more sensitive than all other methods at all read depths.
strikingly, the sensitivity of mapsplice is low, remaining below  <dig>  at all read depths.
in contrast, specificity is over  <dig>  for all methods at all read depths , although ptesfinder achieves the highest specificities averaged across all depths  for both de novo and guided ptes discovery, with all canonical junction reads being correctly identified as such within the simulated data.
only the memczak method has similar specificity when averaged across all read depths .fig.
<dig> sensitivity and specificity in comparisons to other methods.
a sensitivity and  specificity of ptesfinder and  <dig> other publicly available methods  analysed using simulated data .
c mean false discovery rate %  of all methods averaged across all read depths analysed
to compare performance using real data we first reanalysed data from jeck et al.
<cit>  using all  <dig> methods .
to allow direct comparison to ptesfinder, the number of putative circrna structures identified which utilise  <dig> refseq splice sites was recorded for all other methods , as the total numbers include structures from intergenic and intronic regions of the genome.
for all  <dig> samples analysed, ciri consumed >90gb of memory, resulting in incomplete analyses.
it was therefore not analysed further.
of the remaining  <dig> methods, ptesfinder identified on average 15 % more structures than the memczak method and ~70 % more than mapsplice.
the latter is consistent with our finding that mapsplice, which was used in their analysis  <cit> , has low sensitivity at all depths of sequence coverage .
however, circrna_finder reported the highest number of putative circrna transcripts from both exonic and non-exonic regions of the genome, reporting approximately 31–42 % more structures with refseq co-ordinates than ptesfinder .table  <dig>  number of ptes transcripts identified from human fibroblast samples using four methods
memczaka
mapsplicea
circrna_findera
acircrnas utilizing two refseq annotated splice sites in brackets
to investigate the origins of the refseq related structures identified exclusively by circrna_finder, reads defining these structures from  <dig> sample  were re-analysed using ptesfinder .
of  <dig> reads re-analysed, approximately 20 %  are defined as mutilocus or sense-antisense fusions, and a further 19 %  are eliminated by the junctional, genomic, and transcriptomic filters indicating likely false positives .
the remaining 61 %  are not aligned, indicating that their anchors map suboptimally to refseq.
furthermore, plotting the distribution of the number of reads supporting each structure identified by circrna_finder only, by ptesfinder only, and by both methods , revealed that the vast majority of structures identified by circrna_finder alone are supported by a single read.
this is in sharp contrast to structures identified by both methods, or by ptesfinder alone.
while these single-read structures may include bona fide low frequency circrnas, they are also likely to contain false positives caused by suboptimal mapping, consistent with the lower specificity of circrna_finder with our simulated data.fig.
<dig> comparisons with real rnaseq data & published results.
a approximately 64 %  of ptes transcripts utilising  <dig> refseq  splice sites were identified by both circrna_finder and ptesfinder from srr <dig>  read exclusion criteria for ptes transcripts identified by circrna_finder only, when analysed by ptesfinder  distribution of read numbers supporting ptes transcripts identified by circrna_finder only, by ptesfinder only, and by both   ptesfinder identified over 50 %  of transcripts reported in memczak et al.
<cit> .
e the majority of the  <dig> structures reported by memczak et al.
<cit>  but not identified by ptesfinder are intronic or intergenic.
f exclusion criteria for reads presented as evidence for exonic structures in memczak et al.
<cit>  which were not reported by ptesfinder
runtimes for ptesfinder were also 25–35 % lower than for the memczak method, and 50–82 % lower than for mapsplice , but by far the best runtimes were achieved by circrna_finder which utilises the star aligner  <cit> .
these were, however, achieved at higher computing memory cost .
we then used ptesfinder to analyse rnaseq reads previously mined in two further studies  <cit> .
consistent with the above, it identified 13 % more distinct structures from leukocyte and hek <dig> data than were reported by memczak et al.
<cit>  , and  <dig>  % more structures than reported by salzman et al.
<cit>  from leukocyte data ).
as both structures and supporting reads were reported by  <cit> , it was possible to re-analyse the  <dig> structures identified exclusively by their method using ptesfinder.
this established that none correspond to structures which ptesfinder is designed to identify :  <dig>  are derived from intronic, and intergenic regions, and of the  <dig> reads supporting the remaining  <dig> genic structures,  <dig> were excluded by ptesfinder due to low map quality  or multiple map locations ,  <dig> reads were excluded by ptesfinder filters, and the remaining  <dig> possessed at least  <dig> exon boundary which did not map to known splice junctions .
again, while some of these latter reads will undoubtedly correspond to bona fide ptes structures , further blat analysis established that  <dig> mapped in a linear fashion to  <dig> annotated pseudogenes .
approaches to ptes discovery involve a compromise between the ability to detect all potentially rearranged transcripts, and the ability to identify artefacts generated as a result of the sequence and structural complexity of eukaryotic genomes, and of current library construction methods.
it is now clear that the majority of transcripts with re-arranged exon order utilize known exon junctions  <cit>  which are processed by the spliceosome  <cit> .
as a result, methods which utilise existing transcript annotation from the genome under study, such as ptesfinder and those employed by  <cit> , benefit from the reduced noise inherent in this approach and are suited to quantitative analyses of ptes structures that can be characterized using existing annotations.
the use of known/experimentally verified splice sites does reduce the misidentification of template switching artefacts or unspliced pseudogenes as bona fide ptes transcripts.
however, it does mean that not all rearranged transcripts will be identified.
although a recent analysis of human data unconstrained by existing annotation suggests that circrnas which function as mirna sponges are rare  <cit> , discovery of transcripts which do not utilise known splice sites  requires a genome-wide approach unconstrained by existing annotations.
such approaches are, however, inherently more susceptible to artefacts.
the analyses presented above illustrate both the problem of false positive structures, the trade off between sensitivity and specificity in all methods designed to identify rearranged transcripts, and the utility of multiple filters designed to target distinct populations of known artefacts.
CONCLUSIONS
a major challenge in ptes identification is to discriminate between bona fide ptes structures and a wide variety of false positives with distinct origins.
currently, no method which has been used for ptes discovery explicitly excludes all known classes of false positive reads.
to that end, we have developed ptesfinder to identify both linear and circular ptes transcripts from high throughput rnaseq data.
compared to publicly available methods recently used in circrna discovery, ptesfinder achieves higher specificity and sensitivity, and generates output tailored for downstream comparative analyses of transcript abundance, making it an appropriate tool to investigate these rnas within complex mammalian genomes.
availability and requirements
project name: ptesfinder
project home page:http://ptesfinder-v <dig> sourceforge.net/
operating system: linux
programming language: shell, java  <dig>
other requirements: bowtie  <dig> .
<dig> &  <dig> .
<dig>  bedtools  <dig> .0
license: osi-approved open source
any restrictions to use by non-academics: none
additional files
additional file 1:
distinct ptes structures identified from dataset mined in memczak et al.
<cit> .
additional file 2: table s <dig>  analyses of rnaseq data from human fibroblast cells.
additional file 3: figure s <dig>  example reads excluded by filters.
a) reads filtered out by genomic filter for mapping better to pseudogenes & segment duplicated regions b) reads excluded by the transcriptomic filter for having  <dig> % alignment to a canonical splice between exons  <dig> and  <dig> of hnrnph1 c) reads excluded by applying the junctional filters, segment pid and jspan .
additional file 4: table s <dig>  circrna transcripts published by memczak et al.
<cit>  with 100 % overlap to annotated pseudogenes and excluded by ptesfinder.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
mj, de and msk conceived this project; msk and oi developed the software; oi, aa and hma performed in silico analyses of ptes transcripts; oi, msk and mj wrote the paper; de provided guidance and manuscript editing.
all authors read and approved the manuscript.

 recently, there has been an increased interest in a novel class of transcripts where the exon order differs from that found in the genome. once considered cloning artefacts or products of aberrant splicing, it is now established that the majority of these molecules represent circular rna species, although some linear transcripts have been reported. thousands of these novel transcripts have now been identified in a variety of eukaryotic cells, many are conserved across species, suggesting functional relevance, and two have been shown to harbour numerous mirna binding sites and act as mirna sponges to modulate gene expression. recent reports also implicate circrnas in synaptic development and some have expression patterns that correlate with diseases and may act as biomarkers for ageing. despite these reports, the function of the vast majority of these transcripts has not been established. the defining feature of these transcripts at the sequence level is the presence of a splice junction with exons in an order inconsistent with their position in the genome. as this feature alone does not enable inference of structure or mechanistic origins, we use the term post-transcriptional exon shuffled transcripts to collectively describe this population of rna molecules. recent reports have shown that the vast majority of these transcripts emanate from known genes, utilise known splice junctions, and that their biogenesis competes with splicing of canonical transcripts. transcripts arising from ptes specifically exclude chimeric rnas without known splice junctions, and a class of circular rnas comprised of spliced introns, cirnas. many computational methods for identifying chimeric rna molecules from high-throughput rna sequence data have been described. the majority of these are designed to detect heterotypic trans-splicing and fused genes, so are not suitable for ptes detection, or require post processing steps to analyse reads supporting different types of splicing events. recently, however, a number of programs for ptes discovery have been described and used primarily for circrna characterisation. most analyse reads which fail to fully align to a reference sequence, and split these into two or more segments which are then independently aligned to define rearrangements. some make use of gene annotation to guide discovery, while others adopt an unbiased genome-wide approach to capture structures which do not utilise known splice junctions or are non-genic. in addition, the occurrence of ptes can be inferred when two paired end reads map to the transcriptome in a configuration that is not consistent with a linear transcript. the identification of ptes exon junctions within rnaseq data is, however, confounded by known artefacts. false positives can arise from template switching during cdna synthesis, from genes with duplicated exons, from transcription read-through between genes in close proximity due to weak termination signals, and from closely related genes within duplicons or tandem arrays. although experimental enrichment has been combined with informatic approaches to define bona-fide circrnas, many classes of false positive structures are not directly excluded by existing identification methods. for example, reads defining of the most abundant human circrnas reported by memczak et al. map with high sequence identity to the reference sequence and include which are indistinguishable from linear refseq entries. furthermore, a recent experimental analysis of previously identified ptes trancripts concluded that many are template switching artefacts, and template switching predominantly leads to rearragements where the breaks do not occur at splice junctions.fig. examples of intragenic false positives. schematic diagrams showing inferred structure and key sequence relationships for of the most abundant circrnas reported in. in each case, the inferred structure shares 100 % identity to a linear transcript spanning the defining exon-exon junction. within the top hsa_circ_ and show similar patterns of identity to multiple genomic locations. blue – inferred donor exon, red – inferred acceptor exon, black – upstream or downstream refseq exon sharing 100 % identity to donor or acceptor exon at junction. approximate chromosomal locations are shown here we present a method, ptesfinder, that identifies putative ptes structures by mapping rnaseq reads to sequence models generated using existing transcript annotation. it then applies a series of mapping and alignment filters to systematically remove known classes of false positives. it does not make use of paired end mapping information as the lack of intervening sequence precludes such filtering and may affect specificity. we first describe the implementation of this method, and then investigate the effects of different filtering criteria. the program requires certain user adjustable parameters; we therefore also explore systematically the choice of these parameters. finally, through analysis of real and simulated data, we compare ptesfinder to other publicly available methods which have been used to identify circrna transcripts in both cell lines and tissues. implementation pipeline for ptes discovery ptesfinder requires as input files: rnaseq data in fastq format, genomic reference in fasta format, and an annotated transcriptome reference in bed format. the pipeline is split into three phases: a discovery phase to identify putative ptes structures within rnaseq data and define ptes transcript models, an evaluation phase to assess these ptes models, and a filtering phase to exclude potential false positives.fig. overview of ptes discovery pipeline. the workflow includes three major phases: discovery phase, evaluation phase and filtering phase. putative ptes structures discovered using 20 bp anchor reads are evaluated by aligning full fastq reads to the models. the filtering phase includes stringent criteria designed to systematically exclude all known classes of false positive structures discovery phase short sequences from both ends of each fastq read are aligned to the transcriptome reference using bowtie with tolerance for a single mismatch. pairs of anchors from the same read that map to the same gene and in the same orientation, but which map in inverted order with respect to their order in the sequencing read, are then identified. this excludes reads emanating from fused genes and sense-antisense template switching artefacts. retained anchor pairs are then used to determine the exon junctions which define putative ptes events and create sequence models of the inferred products. these constructs are generated by concatenating the last 65 bp of the 5′ exon and the first 65 bp of the 3′ exon. the segment size of is used by default, with the full exon sequence used if an exon is smaller than 65 bp. this parameter is adjustable to accommodate various rnaseq read lengths, and we recommend that it be at least 10 bp shorter than the read length to ensure that only reads mapping across ptes defining junctions are processed in subsequent filtering steps. evaluation phase all the original reads are then re-mapped to ptes constructs generated in the discovery phase. this serves two purposes. first, as rnaseq reads can be short, this enables reads with putative ptes exon junctions close to their termini to be accurately mapped. second, it enables read mapping scores obtained using the ptes constructs to be directly compared to scores obtained from genomic and transcriptomic alignments for filtering purposes. optionally, evaluation can also be ‘guided’ by supplying constructs of previously discovered ptes structures, effectively bypassing the discovery phase. filtering phase to eliminate potential false positives originating from the genome under investigation, all the original reads are mapped to both genomic and transcriptomic references. the number of edits required for alignment, and the number of perfectly aligned base pairs, are used to remove reads which align as well or better to either of these reference sequences than to the ptes constructs. to reduce template switching artefacts, which have heterogeneous junction points within short regions of often imperfect sequence homology, reads which do not align perfectly to the exon junctions which define ptes are also removed using junctional filters. first, a user adjustable minimum junction span parameter is applied to ensure that there are no mismatches or indels within ‘n’ nucleotides either side of the junction position, where n is an even integer. second, to eliminate reads with regions of low quality alignment, a user adjustable segment percent identity parameter is also applied independently to the segments on either side of the ptes junction, such that for a read to be retained both must meet or exceed the specified pid when aligned to the ptes construct. these user adjustable filters rely on alignment summaries provided by the nm field, md field and cigar in the sam files. the output includes the coordinates of the exon end involved in the junctions, a descriptor of the ptes and the number of reads supporting the structure. this is presented in bed format. a second file contains additional information, read counts of all canonical exon junctions from transcripts where a ptes structure has been identified, to facilitate comparison with ptes counts. assessment of pipeline and comparisons to other methods rnaseq data from jeck et al. were analysed at various combinations of jspan and pid. all analyses were performed with and without genomic and transcriptomic filters applied to enable reads discarded by each filtering criterion to be identified. the numbers of ptes structures identified and supporting reads were also recorded. to assess sensitivity and specificity, simulated datasets were generated using all published ptes structures within circbase.org. for each dataset, ptes junctions were randomly selected along with canonical junctions, and constructs were generated for each junction by concatenating the full sequence of both exons involved in each case. 100 bp simulated reads with random start positions within each construct were then generated. scripts published by memczak et al., ciri v., circrna_finder, and mapsplice v.. used in, were compared to ptesfinder by analysing leukocytes cell line rnaseq data, fibroblasts rnaseq data, and simulated data. for each simulation, transcripts correctly identified by each method were determined by comparing genomic coordinates of identified transcripts with the genomic coordinates of transcripts expected to be recovered from within each dataset. the numbers of correctly identified ptes transcripts, incorrectly identified ptes transcripts, ptes transcripts incorrectly excluded, and canonical junctions correctly excluded, were used to estimate sensitivity: tp /, specificity: tn /, and false discovery rate: fp /.  ptesfinder uses established rnaseq tools to identify putative ptes structures, and then systematically excludes known classes of false positive structures by applying genomic, transcriptomic and junctional filters. as an initial assessment of ptesfinder function, rnaseq data from human fibroblast total rna which has previously been mined for circrnas, were analysed both with and without the application of the genomic and transcriptomic alignment filters. reads recovered during analysis, together with alignment edit distances of reads identified by each filter applied seperately, are shown in fig. 3a.fig. summary of reads excluded by filters. a edit distance distribution of reads filtered out by genomic, transcriptomic and junctional filters. inset: seven structures are supported by reads with 100 % alignment but are excluded by the transcriptomic filter. b 30 % of reads filtered out support a false positive structure from  s rrna and are excluded in this plot. c venn diagram showing number of reads excluded by filters. majority of false positive reads are excluded by all three filters. each filter also excludes a distinct population of false positive reads filters target overlapping populations of reads but none is redundant from a total of over million reads analysed, approximately  % have shuffled coordinates with respect to exon position. however, of these only map to ptes sequence constructs generated during the discovery phase, indicating that most of the reads with rearranged anchor pairs do not map to single genes and/or known exon junctions. approximately 85 % of the reads which map to ptes constructs are subsequently removed by the genomic, transcriptomic and junctional filters, with the majority being identified by more than one filter. for instance, over 98 % of reads excluded by the genomic filter are also excluded by the transcriptomic filter, and 60 % of all filtered reads are identified by all three. most of these have high edit distances indicative of low quality alignment. despite this, the genomic, transcriptomic and junctional filters uniquely exclude ~  %, ~  %% and  % of reads mapping to ptes models respectively, indicating that none is wholly redundant. the subset of reads identified specifically by the junctional and transcriptomic filters are defined by low edit distances of between and, although a small number of reads excluded by the transcriptome filter map perfectly to putative ptes constructs with nm = . fig. 3a also reveals a bimodal distribution of mapping qualities for reads excluded by all three filters with peaks at nm =  and nm =  upon manual analysis, most of the excluded reads with nm =  were found to support a false positive structure from  s rrna. comparable rrna derived structures have been identified previously and filtered manually. in fig. 3b, reads supporting this structure have been removed to show the underlying distribution of mapping quality scores. reads excluded by specific filters have different origins to investigate the activity of specific filters further, the mapping co-ordinates of reads removed by the genomic filter were first compared to the co-ordinates of annotated pseudogenes and segmental duplications. this established that ~74 % of reads excluded by the genomic filter had superior alignments to segmental duplications, and ~12 % had superior alignments to pseudogenes. the reads identified by the genomic filter but not by the transcriptomic filter were also found to be enriched for reads derived from segmental duplications and pseudogenes. we next used blat to manually investigate the reads excluded specifically by the transcriptome filter which mapped perfectly to putative ptes constructs. these support putative ptes structures from genes. however, blat analysis established that they all also mapped contiguously with ~100 % identity to the transcriptome due to high sequence identity between neighbouring exons. for example, reads which support a putative single exon ptes structure map with ~100 % identity to exons and of the canonical hnrnph transcript due to high sequence identity between these neighbouring exons. as a result, these reads cannot be taken as supporting evidence for ptes. it is noteworthy that such structures will pass any qualitative filter criterion requiring only unambiguous mapping to ptes constructs, illustrating the value of the transcriptome filter. finally, manual analysis of a subset of the reads identified only by the junctional filters established that these support structures with distinct patterns of suboptimal mapping, such as low alignment quality specific to only one of the two exons in the structure, and low sequence identity specifically at the junction, the latter being consistent with the expected pattern of alignment for template switching artefacts. as one further assessment of the filters, we analysed rnaseq data derived from fibroblast rna which had been pre-digested with rnase r. this selectively removes linear rnas, and enriches for circrnas, and has been shown to significanty increase the recovery of ptes reads. however, we would anticipate that this would also selectively remove false positives derived from pseudogenes and segmental duplications which mimic ptes structures, without necessarily reducing reverse transcription artefacts such as template switching. only ~12 % of reads from the rnaser digested sample which map to ptes sequence constructs are identified by the genomic and transcriptomic filters, compared to 69 % in the undigested sample. furthermore, only 17 % of these map to segmental duplications, compared to 74 % in the undigested sample. in contrast, the proportion of reads excluded by the junctional filters is considerably higher after rnaser digestion, consistent with expectation. pid has greater impact than jspan to investigate the impact of varying the user defined jspan and pid parameters which comprise the junctional filter, the same data was re-analysed using different combinations of these parameters, both with and without the genomic and transcriptome filters applied. this established that varying the pid has a greater impact than varying the jspan, with reads filtered with maximal pid and lowest jspan compared to only reads filtered with the maximal jspan at lowest pid. furthermore, varying the pid between 60 % and 75 % has little impact at any jspan value, but above 75 % there is a linear relationship with the number of reads filtered. as the default junctional filter parameters failed to identify some reads excluded by the other filters, this analysis was repeated using only these reads to establish the jspan and pid parameters required to identify them. over 99 % of these reads are excluded with the most stringent junctional filter parameters. furthermore, the vast majority are filtered with a pid of 85 %, suggesting this is a logical setting for this parameter. the jspan setting only has a major impact at low pids.fig. effect of varying junctional filter parameters. a number of reads passing filter at different combinations of jspan and pid. b percentage of reads only excluded by transcriptome and genomic filters at default settings, which are filtered at different combinations of jspan and pid specificity, sensitivity and comparison with other methods to assess the sensitivity and specificity of the pipeline and compare it to other methods, simulated reads were generated from previously identified ptes and associated canonical transcripts, and analysed at various read depths of coverage using default parameters. in addition to assessing ptesfinder for de novo ptes discovery, the use of constructs of previously reported structures for guided discovery was also assessed, as were four publicly available methods which have previously been employed to identify circrna transcripts: mapsplice v. used in, ciri v., circrna_finder and the method used by memczak et al.. results from simulated datasets are presented in fig. 5a-c, and illustrate that sensitivity varies considerably with coverage, and between methods. at read coverage of the sensitivity of ptesfinder is below. this can be attributed to ptes junctions occurring within the terminal 20 bp of reads, as the low tolerance for mismatches during anchor mapping will result in their elimination. however, sensitivity reaches >90 % at coverage of or higher for both guided and unguided analyses, with guided ptesfinder being equally or more sensitive than all other methods at all read depths. strikingly, the sensitivity of mapsplice is low, remaining below at all read depths. in contrast, specificity is over for all methods at all read depths, although ptesfinder achieves the highest specificities averaged across all depths for both de novo and guided ptes discovery, with all canonical junction reads being correctly identified as such within the simulated data. only the memczak method has similar specificity when averaged across all read depths.fig. sensitivity and specificity in comparisons to other methods. a sensitivity and specificity of ptesfinder and other publicly available methods analysed using simulated data. c mean false discovery rate % of all methods averaged across all read depths analysed to compare performance using real data we first reanalysed data from jeck et al. using all methods. to allow direct comparison to ptesfinder, the number of putative circrna structures identified which utilise refseq splice sites was recorded for all other methods, as the total numbers include structures from intergenic and intronic regions of the genome. for all samples analysed, ciri consumed >90gb of memory, resulting in incomplete analyses. it was therefore not analysed further. of the remaining methods, ptesfinder identified on average 15 % more structures than the memczak method and ~70 % more than mapsplice. the latter is consistent with our finding that mapsplice, which was used in their analysis, has low sensitivity at all depths of sequence coverage. however, circrna_finder reported the highest number of putative circrna transcripts from both exonic and non-exonic regions of the genome, reporting approximately 31–42 % more structures with refseq co-ordinates than ptesfinder.table  number of ptes transcripts identified from human fibroblast samples using four methods memczaka mapsplicea circrna_findera acircrnas utilizing two refseq annotated splice sites in brackets to investigate the origins of the refseq related structures identified exclusively by circrna_finder, reads defining these structures from sample were re-analysed using ptesfinder. of reads re-analysed, approximately 20 % are defined as mutilocus or sense-antisense fusions, and a further 19 % are eliminated by the junctional, genomic, and transcriptomic filters indicating likely false positives. the remaining 61 % are not aligned, indicating that their anchors map suboptimally to refseq. furthermore, plotting the distribution of the number of reads supporting each structure identified by circrna_finder only, by ptesfinder only, and by both methods, revealed that the vast majority of structures identified by circrna_finder alone are supported by a single read. this is in sharp contrast to structures identified by both methods, or by ptesfinder alone. while these single-read structures may include bona fide low frequency circrnas, they are also likely to contain false positives caused by suboptimal mapping, consistent with the lower specificity of circrna_finder with our simulated data.fig. comparisons with real rnaseq data & published results. a approximately 64 % of ptes transcripts utilising refseq splice sites were identified by both circrna_finder and ptesfinder from srr read exclusion criteria for ptes transcripts identified by circrna_finder only, when analysed by ptesfinder distribution of read numbers supporting ptes transcripts identified by circrna_finder only, by ptesfinder only, and by both ptesfinder identified over 50 % of transcripts reported in memczak et al.. e the majority of the structures reported by memczak et al. but not identified by ptesfinder are intronic or intergenic. f exclusion criteria for reads presented as evidence for exonic structures in memczak et al. which were not reported by ptesfinder runtimes for ptesfinder were also 25–35 % lower than for the memczak method, and 50–82 % lower than for mapsplice, but by far the best runtimes were achieved by circrna_finder which utilises the star aligner. these were, however, achieved at higher computing memory cost. we then used ptesfinder to analyse rnaseq reads previously mined in two further studies. consistent with the above, it identified 13 % more distinct structures from leukocyte and hek data than were reported by memczak et al., and  % more structures than reported by salzman et al. from leukocyte data ). as both structures and supporting reads were reported by, it was possible to re-analyse the structures identified exclusively by their method using ptesfinder. this established that none correspond to structures which ptesfinder is designed to identify: are derived from intronic, and intergenic regions, and of the reads supporting the remaining genic structures, were excluded by ptesfinder due to low map quality or multiple map locations, reads were excluded by ptesfinder filters, and the remaining possessed at least exon boundary which did not map to known splice junctions. again, while some of these latter reads will undoubtedly correspond to bona fide ptes structures, further blat analysis established that mapped in a linear fashion to annotated pseudogenes. approaches to ptes discovery involve a compromise between the ability to detect all potentially rearranged transcripts, and the ability to identify artefacts generated as a result of the sequence and structural complexity of eukaryotic genomes, and of current library construction methods. it is now clear that the majority of transcripts with re-arranged exon order utilize known exon junctions which are processed by the spliceosome. as a result, methods which utilise existing transcript annotation from the genome under study, such as ptesfinder and those employed by, benefit from the reduced noise inherent in this approach and are suited to quantitative analyses of ptes structures that can be characterized using existing annotations. the use of known/experimentally verified splice sites does reduce the misidentification of template switching artefacts or unspliced pseudogenes as bona fide ptes transcripts. however, it does mean that not all rearranged transcripts will be identified. although a recent analysis of human data unconstrained by existing annotation suggests that circrnas which function as mirna sponges are rare, discovery of transcripts which do not utilise known splice sites requires a genome-wide approach unconstrained by existing annotations. such approaches are, however, inherently more susceptible to artefacts. the analyses presented above illustrate both the problem of false positive structures, the trade off between sensitivity and specificity in all methods designed to identify rearranged transcripts, and the utility of multiple filters designed to target distinct populations of known artefacts.  a major challenge in ptes identification is to discriminate between bona fide ptes structures and a wide variety of false positives with distinct origins. currently, no method which has been used for ptes discovery explicitly excludes all known classes of false positive reads. to that end, we have developed ptesfinder to identify both linear and circular ptes transcripts from high throughput rnaseq data. compared to publicly available methods recently used in circrna discovery, ptesfinder achieves higher specificity and sensitivity, and generates output tailored for downstream comparative analyses of transcript abundance, making it an appropriate tool to investigate these rnas within complex mammalian genomes. availability and requirements project name: ptesfinder project home page:http://ptesfinder-v sourceforge.net/ operating system: linux programming language: shell, java other requirements: bowtie. &. bedtools.0 license: osi-approved open source any restrictions to use by non-academics: none additional files additional file 1: distinct ptes structures identified from dataset mined in memczak et al.. additional file 2: table s analyses of rnaseq data from human fibroblast cells. additional file 3: figure s example reads excluded by filters. a) reads filtered out by genomic filter for mapping better to pseudogenes & segment duplicated regions b) reads excluded by the transcriptomic filter for having % alignment to a canonical splice between exons and of hnrnph1 c) reads excluded by applying the junctional filters, segment pid and jspan. additional file 4: table s circrna transcripts published by memczak et al. with 100 % overlap to annotated pseudogenes and excluded by ptesfinder. competing interests the authors declare that they have no competing interests. authors’ contributions mj, de and msk conceived this project; msk and oi developed the software; oi, aa and hma performed in silico analyses of ptes transcripts; oi, msk and mj wrote the paper; de provided guidance and manuscript editing. all authors read and approved the manuscript. 
uniquely, it systematically excludes potential artefacts emanating from pseudogenes, segmental duplications, and template switching, and outputs both ptes and canonical exon junction counts to facilitate comparative analyses. with high sensitivity and specificity, user-adjustable filters that target known sources of false positives, and tailored output to facilitate comparison of transcript levels, ptesfinder will facilitate the discovery and analysis of these poorly understood transcripts. ptesfinder also identifies between 13 % and  % more structures, compared to publicly available methods recently used to identify human circular rnas. these have been identified in a wide variety of tissues and cell types from many eukaryotes, and are now known to be mostly circular, cytoplasmic, and non-coding. accurate identification of these transcripts can, however, be difficult due to artefacts from a wide variety of sources. here, we present a computational method, ptesfinder, to identify these transcripts from high throughput rnaseq data.
 transcripts, which have been subject to post-transcriptional exon shuffling, have an exon order inconsistent with the underlying genomic sequence. these have been identified in a wide variety of tissues and cell types from many eukaryotes, and are now known to be mostly circular, cytoplasmic, and non-coding. although there is no uniformly ascribed function, several have been shown to be involved in gene regulation. accurate identification of these transcripts can, however, be difficult due to artefacts from a wide variety of sources.  here, we present a computational method, ptesfinder, to identify these transcripts from high throughput rnaseq data. uniquely, it systematically excludes potential artefacts emanating from pseudogenes, segmental duplications, and template switching, and outputs both ptes and canonical exon junction counts to facilitate comparative analyses. in comparison with four existing methods, ptesfinder achieves highest specificity and comparable sensitivity at a variety of read depths. ptesfinder also identifies between 13 % and  % more structures, compared to publicly available methods recently used to identify human circular rnas.  with high sensitivity and specificity, user-adjustable filters that target known sources of false positives, and tailored output to facilitate comparison of transcript levels, ptesfinder will facilitate the discovery and analysis of these poorly understood transcripts. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords rnaseqcircrnaptessplicingmrna processingsoftwarehttp://dx.doi.org/ /501100000268biotechnology and biological sciences research council bb/j014516/1izuogu osagie g. http://dx.doi.org/ /501100000275leverhulme trust rpg-2012-795jackson michael s. http://dx.doi.org/ /100004440wellcome trust wt089225maelliot david j. issue-copyright-statement© the author 2016 