Vendrell-Mir et al. Mobile DNA (2019) 10:53 https://doi.org/10.1186/s13100-019-0197-9 METHODOLOGY Open Access A benchmark of transposon insertion detection tools using real data Pol Vendrell-Mir1, Fabio Barteri1, Miriam Merenciano2, Josefa González2, Josep M. Casacuberta1* and Raúl Castanera1* Abstract Background: Transposable elements (TEs) are an important source of genomic variability in eukaryotic genomes. Their activity impacts genome architecture and gene expression and can lead to drastic phenotypic changes. Therefore, identifying TE polymorphisms is key to better understand the link between genotype and phenotype. However, most genotype-to-phenotype analyses have concentrated on single nucleotide polymorphisms as they are easier to reliable detect using short-read data. Many bioinformatic tools have been developed to identify transposon insertions from resequencing data using short reads. Nevertheless, the performance of most of these tools has been tested using simulated insertions, which do not accurately reproduce the complexity of natural insertions. Results: We have overcome this limitation by building a dataset of insertions from the comparison of two high- quality rice genomes, followed by extensive manual curation. This dataset contains validated insertions of two very different types of TEs, LTR-retrotransposons and MITEs. Using this dataset, we have benchmarked the sensitivity and precision of 12 commonly used tools, and our results suggest that in general their sensitivity was previously overestimated when using simulated data. Our results also show that, increasing coverage leads to a better sensitivity but with a cost in precision. Moreover, we found important differences in tool performance, with some tools performing better on a specific type of TEs. We have also used two sets of experimentally validated insertions in Drosophila and humans and show that this trend is maintained in genomes of different size and complexity. Conclusions: We discuss the possible choice of tools depending on the goals of the study and show that the appropriate combination of tools could be an option for most approaches, increasing the sensitivity while maintaining a good precision. Keywords: Benchmark, Transposable elements, Polymorphism, Transposon insertion, Resequencing Background insertion close to genes can generate new gene expression Transposable elements (TEs) constitute a very important patterns. In addition, TEs, and in particular LTR- fraction of eukaryotic genomes, and their ability to trans- retrotransposons and MITEs (Miniature Inverted Trans- pose, excise and produce complex genomic rearrangements posable Elements), have been shown to contain transcrip- make them a key source of genomic diversity. Previous tion factor binding sites, which can be mobilized by work done over the last decades has uncovered their enor- transposition rewiring new genes into pre-existing tran- mous potential as gene regulators, a role that TEs play scriptional networks [5, 12, 20]. As a consequence, TEs through a variety of genetic and epigenetic mechanisms have the potentiality to generate important genomic and [12, 43]. Certain TEs, such as Long Terminal repeat (LTR)- transcriptional variability, and the interest in these elements retrotransposon carry their own promoters, and their has drastically increased in the last years. Due to their repetitive nature and their sequence * Correspondence: josep.casacuberta@cragenomica.es; diversity, the annotation of TEs is more complex than raul.castanera@cragenomica.es that of protein coding genes. Nevertheless, thanks to the 1 Centre for Research in Agricultural Genomics CSIC-IRTA-UAB-UB, Campus UAB, Edifici CRAG, Bellaterra, 08193 Barcelona, Spain development of tools such as Repeatmasker (http:// Full list of author information is available at the end of the article www.repeatmasker.org) and sophisticated pipelines such © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 2 of 19 as REPET [16], methodologies of TE detection and an- repetitive genomes is uncertain. This is especially rele- notation in assembled genomes are today robust. The vant as S. cerevisiae contains only 51 full LTR- availability of high-quality reference genomes coupled retrotransposons in the whole genome [8], whereas in with the exponential increment of resequencing data has most plant and animal genomes the LTR- boosted our capacity to evaluate intraspecific variability. retrotransposon load is several orders of magnitude By obtaining accurate maps of genetic variation, charac- higher. Also, a recent study focused on simulated but terizing the genetic basis of phenotypic variance is now also real human AluY, L1 and SVA families revealed possible at a genome-wide scale thanks to association huge differences in the ability of seven tools to detect studies (GWAS). Until now, most of the efforts have been TIPs [41]. In spite of the importance of these families focused on analyzing the variability at the nucleotide level for human research, they do not represent the diversity (SNPs, single nucleotide polymorphisms), as there are ro- of the TE landscape of other animals and plants, which bust algorithms to perform variant calling. However, TEs is far more complex, with many families from different generate an important part of the genetic variability orders being potentially active, and where the amount of present in a particular species. Moreover, the timing of oc- truncated non-autonomous elements greatly outnum- currence of TE and SNP mutations is different, as the bers the active copies. former can amplify in bursts generating a great amount of In plants, TEs are at the origin of important agronomic diversity in short periods of time, whereas SNP mutation traits, such as apical dominance in maize [45], the skin rates are more constant in time. Therefore, the identifica- and flesh colors in grape [28] and blood oranges [4]. Dif- tion of Transposon Insertion Polymorphisms (TIPs) is of ferent efforts have been made recently to identify TIPs high interest. Nevertheless, our capacity to accurately that could be responsible for important variability in identify TIPs using re-sequencing data is hampered by the plants. Carpentier et. al [7] screened the presence of 32 structural complexity of TEs. rice LTR-retrotransposon families in the 3000-rice gen- In the last few years, many laboratories have developed ome dataset and uncovered more than 50,000 TIPs, most bioinformatic tools to look for TIPs and have started to of them occurring at a very low frequency, which is indi- analyze their impact in intra-species variability, including cative of recent activity. Besides LTR-retrotransposons, crop plants [7, 10, 42]. There are two main approaches MITEs are probably the most prevalent group of transpo- that can be used to detect TIPs in whole-genome se- sons in plants, including rice, where they have experienced quence data: i) inference from discordant read-pair map- recent massive amplification bursts [10, 35]. MITEs are pings, and ii) clustering of ‘split’ reads sharing common structurally very different from LTR-retrotransposons, as alignment junctions [2, 15]. Most of the recently devel- they are non-autonomous, usually non-coding, and rela- oped tools incorporate both methodologies, and in some tively small. They are of particular interest because they cases TIPs have been experimentally validated [27]. tend to integrate close to genes and may carry regulatory Moreover, in some cases the authors have evaluated domains [20], having the potential to create or rewire their sensitivity and precision (also known as positive regulatory networks [12]. In the present study, we have predictive value) [11, 24]. However, in most cases these taken advantage of the existence of several high-quality as- evaluations were performed by generating simulated in- sembled genomes of different rice varieties to create a vali- sertions that are randomly placed in the genome, and dated dataset of natural LTR-retrotransposon and MITE then used to compare with tool predictions. Simulated insertions obtained by direct comparison between the as- insertions are far from representing the complexity of sembled genomes (Nipponbare and MH63), that we have “natural” TIPs, as many of their features are difficult or used to benchmark the performance of 12 TIP calling impossible to mimic accurately (i.e.: element degener- tools. Moreover, we have also analyzed the sensitivity of ation, nested insertions, insertion preferences, etc.). As a the best performing tools to detect experimentally vali- consequence, the benchmarks done with simulated data dated TIPs in Drosophila and humans. Our results evi- tend to overestimate performance of the tools analyzed dence that tool performance is in general lower than [21]. An example of such benchmarks is the one re- estimated by previous simulations, and highly variable de- ported by the developers of McClintock, a pipeline that pending on sequencing coverage and TE type. Also, we integrates six tools [36] (Table 1). In their study, the au- show that an appropriate combination of tools can in- thors provided a detailed comparison of their compo- crease the sensitivity of predictions while maintaining high nent’s performance in sensitivity and positional accuracy precision levels. based on simulated LTR-retrotransposon insertions, which also includes some real resequencing data, in the Results yeast Saccharomyces cerevisiae. In spite of the interest of Tools selected for benchmarking such comparative analysis, the direct translation of these We selected 12 of the most widely used tools for the de- results to other eukaryotic models with bigger and more tection of TIPs (Table 1). Among them, four were Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 3 of 19 Table 1 Tools selected for the benchmark of TE insertions Tool Target Prediction Input Output Perceived difficulty Manual format Installation Input preparation RelocaTE2 Non-reference All families fastq gff file Easy Easy https://github.com/JinfengChen/RelocaTE2 insertions Jitterbug Non-reference All families Bam gff file Medium Medium https://github.com/elzbth/jitterbug insertions Retroseq a Non-reference All families Bam vcf file Easy Difficult https://github.com/tk2/RetroSeq/wiki/ insertions RetroSeq-Tutorial ITIS Non-reference Single-family fastq Bed file Easy Medium https://github.com/Chuan-Jiang/ITIS insertions MELT Reference and non- Single-family Bam vcf file Easy Medium http://melt.igs.umaryland.edu/manual.php reference insertion PopoolationTE2 Reference and non- All families fastq Tool-specific Easy Easy https://sourceforge.net/projects/ reference insertions popoolation-te2/ Teflon Reference and non- All families fastq Tool-specific Medium Medium https://github.com/jradrion/TEFLoN reference insertions Trackposon Reference and non- Single-family fastq Bed file Easy Easy http://gamay.univ-perp.fr/~Panaudlab/ reference insertions TRACKPOSON.tar.gz TEMP a Reference and non- All families Bam Tool-specific Easy Difficult https://github.com/JialiUMassWengLab/ reference insertions TEMP/blob/master/Manual TE-locate a Reference and non- All families Sam Tool-specific Easy Difficult https://sourceforge.net/projects/te-locate/ reference insertions a PopoolationTE Reference and non- All families fastq Tool-specific Easy Difficult https://sourceforge.net/p/popoolationte/ reference insertions wiki/Workflow/ ngs_te_mapper a Reference and non- All families fastq Bed file Easy Difficult https://github.com/bergmanlab/ reference insertions ngs_te_mapper McClintock Reference and non- All families fastq Bed file Easy Difficult https://github.com/bergmanlab/mcclintock reference insertion a These tools were run as part of the McClintock pipeline. Perceived difficulty refers to McClintock and not the original methods Installation: Easy = available in Conda, or automatic / semi-automatic installation. Medium = Needs several dependencies or specific versions of packages that need manual installation. Input preparation: Easy = can be run using common formats (ie fasta, bed) without the need of specific formatting. Medium = Needs specific formatting. Difficult = Needs very specific formatting specifically designed to detect non-reference insertions cluster. This can be sometimes complex due to the num- (not present in the reference genome) (RelocaTE2 [11], ber of different dependencies, and especially due to the Jitterbug [21], Retroseq [27] and ITIS [24]), and eight specificity of the input file preparation. In this regard, were able to detect reference (present in the reference we found that RelocaTE2, PopoolationTE2 and Trackpo- genome) and non-reference insertions (MELT [18], son were the less problematic tools (Table 1). One possi- Popoolation TE2 [29], Teflon [1], Trackposon [7], TEMP bility that would make the installation of these tools [48], TE-locate [37], Popoolation TE [30], and ngs_te_ much easier would be to have them integrated in an en- mapper [32]. Tools specifically designed to detect pres- vironment such as Conda. This is a possibility that fu- ence/absence of reference TE insertions in re-sequenced ture developers should take into account. genomes (i.e.: T-lex 3) [3] were not benchmarked here. In addition to their different targets, some of the tools LTR-retrotransposon and MITE landscape in Nipponbare were family-specific (meaning that they run with one TE and MH63 genomes family at a time only), whereas most of them are able to In order to perform a benchmarking exercise that could detect insertions from all the families in the same run be representative of as much as possible TIP detection (broad-spectrum). Five out of the 12 tested tools were in eukaryotes, we decided to use rice as a model as it run as components of McClintock, a pipeline that com- has a genome of 430 Mb, which is relatively large and bines the use of several TIP detection tools and stan- complex in terms of TE landscape, and that has already dardizes their outputs into the commonly used BED been considered as being as close as possible to a repre- format (Table 1). sentative genome for angiosperms [7]. Moreover, there The first difficulty that the user has to face is properly are several good-quality assemblies and short-read data- installing and making the tools run, often in a computer sets of rice varieties available [23, 47]. In terms of the Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 4 of 19 TEs to be detected, we concentrated on LTR- Annotation of standard transposon insertion datasets for retrotransposons and MITEs as, in addition to be the tool benchmarking most prevalent TE types in plant genomes, they are The most straightforward way of identifying an insertion functionally and structurally very different. Indeed, polymorphism “in silico” when two high quality assem- whereas LTR-retrotransposons are relatively long ele- bled genomes are available (as it is here the case), is by ments (typically several Kb-long) and contain many aligning orthologous loci. To identify the Nipponbare structural features relatively easy to detect (e.g.: long orthologous loci to those that in MH63 contain a TE in- LTRs at their extremities, coding capacity for several sertion, we mapped the flanking regions of each MH63 well conserved enzymatic activities), MITEs are short full-length LTR-retrotransposon and MITE insertion (typically 100–800 nt), are non-coding and do not con- against the Nipponbare genome. As sequence diversity tain structural features (except for short inverted-repeats and structural differences between the two genomes may in most cases) allowing for structural detection. complicate this analysis, we tested different flanking se- We used a combination of structural and quence lengths and found that 500 nt was the one that homology-based approaches to annotate a high- allow to identify more reference and non-reference in- quality dataset of 3733 and 3787 full-length LTR- sertions (Additional file 6: Figure S1). By inspecting the retrotransposons in Nipponbare and MH63 (Minghui distance between the two mapped flanks, we could as- 63) assemblies, respectively (Table 2). These elements sign the orthology status to the locus (ie, empty site or contain intact Target Site Duplications (TSDs), Long full site). Using this approach, we were able to assign an Terminal Repeats as well as coding domains. All of orthology status to 86% of the MITE loci, but only to them were clustered at 80% similarity over 80% 41% of the LTR-retrotransposons loci. This was probably length to obtain families and we derived a consensus due to the difficulty to identify the orthologous loci of for each family. RepeatMasker was then run with insertions siting in repetitive sequences, which is much such consensuses to identify all the LTR- more frequent for LTR-retrotransposons than for retrotransposon copies of the genome (including frag- MITEs. Therefore, although this strategy seems the ments and non-autonomous elements) related to the more straightforward, it has clear limitations. Moreover, characterized families. A similar strategy was used to as defining the precise TE-genome junctions for non-full identify ~ 46,000 full-length MITEs, as well as ~ 200, length elements (ie, degenerated or partial elements, 000 partial MITE copies (see methods section). which are the vast majority of LTR-retrotransposons, Whereas full-length LTR-retrotransposons represent a Table 1) is challenging, we could not use this strategy to very small proportion of the total number of LTR- analyze the possible polymorphisms arising from non- retrotransposon copies detected, (3%, Table 2), full- full-length LTR-retrotransposons. To overcome those length MITEs represent an important fraction (23%). limitations and increase the dataset of curated insertions, The distribution along the chromosomes of the two we developed a strategy aimed at complementing the transposon groups is also different, with LTR- TIPs dataset with TIPs predicted with the 12 tools ana- retrotransposons being more abundant in the centro- lyzed here (Table 2), which were individually validated. meric and pericentromeric regions and MITEs To this end we ran the different TIP-prediction tools populating evenly the rest of the chromosome using MH63 paired-end reads mapped to the Nippon- (Fig. 1). bare reference genome. We divided the Nipponbare genome in 500 nt windows and mapped the windows Table 2 Annotation of LTR-retrotransposons and MITEs in rice containing predicted insertions (red boxes, Fig. 2) to the assemblies MH63 genome. An inspection of the aligned sections TE Classification Nipponbare MH63 allowed determining whether the predicted insertion LTR-all a 131,905 117,362 corresponded to a reference (shared) or non-reference LTR full-length b 3733 3787 (MH63 specific) insertion or if it should be considered a LTR- Gypsy 1354 1303 false positive (Fig. 2b). Indeed, in case of reference (shared) insertions, the Nipponbare and the correspond- LTR- Copia 944 759 ing MH63 sequences would perfectly align, showing that c LTR- Unclassified 1435 1725 the sequence, which contains a TE insertion is MITE-all (1) 211,732 191,113 conserved in both genomes (Fig. 2b, left); in case of a MITE full-length d 45,963 46,725 non-reference (MH63 specific) insertion, the alignment a Repeatmasker fragments. Includes both intact and truncated elements will be split by an insertion in the MH63 sequence cor- b c High confidence elements containing intact LTRs, TSDs and coding domains responding to an annotated TE (Fig. 2b, right); and in Intact elements whose poor coding domain conservation doesn’t allow proper classification case where the two sequences show a continuous align- d Elements spanning more than 80% of its family consensus length ment in the absence of an annotated TE insertion in Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 5 of 19 A 10 20 MITEs 0 30 asker Teflon ll-length q rtions te Retrose E2 E2 TEloca epeatM lationT relocaT f. Inse MITE fu Popoo MITE R Non-re Predictions B 10 20 0 30 LTR- retrotransposons sker ngth Teflon q rtions te Retrose eatMa E2 E2 TEloca ll-le lationT relocaT f. Inse LTR fu ep LTR R Popoo Non-re Predictions Fig. 1 Density of MITEs (a) and LTR-retrotransposons (b) along the rice chromosome 5 (window size = 50 Kb). Black circles represent centromeres. Track 1 shows the density of all elements annotated in the chromosome by RepeatMasker. Track 2 shows the density of full-length elements. Track 3 shows the density of validated non-reference insertions (MH63-specific insertions) in the benchmarking standard. Tracks 4–8 show the density of non-reference predictions of five tools Nipponbare, this will indicate that the TE prediction is a performed, and we decided to restrict the dataset of false positive (Fig. 2b, middle). After running all tools, LTR-retrotransposons to a single chromosome (chr5). adjacent windows corresponding to TIP predictions of This strategy combined the power of detection of the same category were merged to produce a final read-based methods (useful for uncovering polymor- dataset. LTR-retrotransposon insertions are frequently phisms derived from both full and degenerated ele- more complex than MITEs (i.e.: length, tendency to ments), with the reliability of the validation based on form nested insertions and extremely high amount of alignments between high-quality assembled genomes. By truncated and degenerated elements, Table 2). Because using this combined approach, we increased the number of this, it was difficult in many cases to automatically of validated non-reference MITE insertions from 1898 validate the insertions. Therefore, manual inspection of to 3117 whereas for LTR-retrotransposons (chr5) the the alignments of LTR-retrotransposons TIPs was amount of non-reference insertions in our validated Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 6 of 19 dataset of TE insertions that were not detected by the evaluated tool were counted as False Negatives (FN). Most of the tools analyzed here are able to detect insertions from all the families in the same run (broad- spectrum). Some of these tools are able to detect refer- ence and non-reference insertions, whereas others only detect non-reference insertions. The programs use dif- ferent strategies to identify these two types of insertions, and consequently we analyzed their performance separately. Detection of reference insertions by broad-spectrum tools We observed that whereas the precision detecting MITE and LTR-retrotransposon reference insertions was very high for both types of elements, the sensitivity levels of most of the tools were much higher for MITEs (Fig. 3). For MITEs, the sensitivity of most tools increased with coverage and tended to stabilize at 20-40X coverage (Fig. 3a). Teflon had consistently the best sensitivity and overall performance (F1-score) in the detection of refer- Fig. 2 Individual validation of predicted insertions. Black boxes ence MITE insertions even at low coverage, reaching a represent TE annotations in Nipponbare IRGSP (green rectangle) and sensitivity of 74% at 10X with an almost 100% precision MH63 (blue rectangle) assembled genomes. Examples of shared (Fig. 3a). All tools showed precision levels higher than (reference) and MH63-specific (non-reference) insertions are shown 99% at all coverages, and all tools except ngs_te_mapper in a. Insertions predicted by each tool (shown as arrows in b) were intersected with windows of 500 bp spanning the entire Nipponbare yielded a sensitivity higher than 60% at 40X (Fig. 3a, IRGSP genome, and windows having an intersection (red boxes, b) Additional file 3: Table S2). By contrast, the sensitivity at were aligned to MH63 genome. True positive reference insertions 5X was low in general, with Teflon being the only tool (TP ref.) were those having full-length alignments with an MH63 reaching more than 50% (Fig. 3a). region where a MITE or LTR-retrotransposon was annotated. False Regarding the detection of reference LTR- positives (FP) have high-quality alignments (MAQ > 30) to regions were no MITE or LTR-retrotransposon was present. True positive retrotransposons, the general tool performance was non-reference insertions (TP non-ref) alignments were those having much lower than for MITEs (Fig. 3b). In this case, TE- a spliced alignment in which the two hits were separated by a locate reached the maximum sensitivity followed by region that overlaps with a MITE or LTR-retrotransposon annotated Teflon and was only slightly higher than 50% (Fig. 3b), in MH63 and the other tools remained below 40% sensitivity. The sensitivity of TE-locate was higher than 50% in all the dataset increased from 22 to 239 (Additional file 2: Table coverages, whereas in Teflon, PopoolationTE2 and Popoo- S1). The result was a high-quality dataset of True lationTE it increased with coverage (Fig. 3b). When we Positive (TP) and False Positive (FP) reference and non- focused only on the detection of full-length LTR- reference insertions (Additional file 2: Table S1). In retrotransposons, the performance of all tools increased addition, there were predicted insertions that did not considerably, reaching a maximum sensitivity of 85.4% match neither with TP nor FP (i.e.: cases that did not fit (Fig. 3c). TE-locate was again the best performer showing in the scenarios described in Fig. 2b). We analyzed the a sensitivity over 80% for all the coverages. We excluded specific cases of unclassified non-reference insertions the predictions of TEMP for reference insertions, as this and found that 86% of these LTR-retrotransposon pre- tool is based on the detection of absences assuming the dicted TIPs and 92% of such MITE TIPs overlapped presence as default, which leads to an overestimation of with other transposons annotated in the reference. the number of insertions, especially at a very low These cases were not used for downstream analyses, as coverage. most tools specifically indicate in their manuals that they cannot properly detect nested insertions. In order to Detection of non-reference insertions by broad-spectrum evaluate the performance of each tool, we intersected tools the windows corresponding to the TE insertions pre- All the benchmarked tools are able to detect non- dicted by the tool (both reference and non-reference TE reference insertions, a task that is more challenging than insertions) with those of the curated dataset to identify detecting reference insertions, as the former are not TP and FP (Fig. 2b). Insertions present in the curated present in the reference genome to which the reads are Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 7 of 19 Fig. 3 Performance of broad-spectrum tools in the detection of reference insertions of MITEs (a), all LTR-retrotransposons (b) and full-length LTR-retrotransposons (c) Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 8 of 19 mapped. In this case sensitivity was strongly dependent a high number of false positives, which translates into on coverage (Fig. 4). Precision was very different for low precision levels (Fig. 4b). In general, the precision MITE and LTR-retrotransposon predictions, showing a detecting LTR-retrotransposons with respect to MITEs tendency to decrease at high coverage (Fig. 4). Regarding was much lower for all tools. Jitterbug was the only pro- MITEs, Teflon was the best performer followed by gram with a moderate precision (> 75%) across all cover- PoPoolationTE2 and Retroseq (Fig. 4a). These tools age levels, although its sensitivity was low (maximum of reached a sensitivity close to 75% (up to 75.6% in 40X 32.7% at 40X) (Fig. 4b). According to the F1-score, coverage for Teflon), whereas the rest of the tools had a PopoolationTE2 and Teflon were the best performers at much lower sensitivity (Fig. 4a). The precision was very low coverages (5X-10X), whereas at higher coverages high (> 95%) for most tools with the exception of TE- PopoolationTE2 and Jitterbug showed the best balance locate, which dropped from 92.5% in 5X to 75.6% in between sensitivity and precision (Fig. 4b). Differently to 40X. All the tools improved their performance when the what we previously did for reference insertions, we did coverage increased (except Jitterbug, which performed not compute the performance of the tools using only the best at 20X), with PopoolationTE2 and Retroseq full-length LTR-retrotransposons because they represent showing the steepest increase, especially between 5X and only a small fraction of the non-reference annotated 20X (Fig. 4a). insertions. Regarding LTR-retrotransposons, PopoolationTE2 The output of most tools contains information that achieved the highest sensitivity, reaching a maximum of can be used for filtering the putative insertions to 88.5% at 40X (Fig. 4b). Nevertheless, these tools yielded achieve more precise detection levels. We checked Fig. 4 Performance of broad-spectrum tools in the detection of non-reference insertions of MITEs (a) and LTR-retrotransposons (b). Relocate2 on LTR-retrotransposons at 40X was killed after 5 days running with 8 CPUs and 64GB of RAM Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 9 of 19 different filters for each program looking for gains in of the detected non-reference true and false positives for precision with a low cost in sensitivity. In some cases, the five better performing tools for MITE or LTR- such as Jitterbug, the precision was already very high, retrotransposon TIP predictions (40X), taking into ac- and the filtering was not needed. In others, the cost in count their sensitivity and precision. In spite of the dif- sensitivity was too high and the filtering was not consid- ference in the amount of predictions between MITEs ered useful. For the two best-performing tools, Popoola- and LTR-retrotransposons, the results showed very tionTE2 and Teflon, filtering did result in significative similar trends: 54% of TP were detected only by one tool gains in precision without an excessive cost in sensitiv- for both MITE and LTR-retrotransposon insertions ity. For PopoolationTE2 we applied a zygosity filter of (Fig. 6). As expected, the FP detected were tool-specific 0.7 (based on the fraction of reads supporting the inser- in the vast majority of the cases (90.2% were detected by tion) which led to a drop of sensitivity for both MITEs only one tool for MITEs and 98% for LTR- (from 76 to 63%) and LTR-retrotransposons detection retrotransposons). The number of insertions detected by (from 88 to 65%, Additional file 7: Figure S2), but with all tools was very low (1.3% of all TIPs detected for an increase of precision, which was particularly striking MITEs and 1.4% for LTR-retrotransposons). These re- for LTR-retrotransposons (from 28.9 to 91.9% at 40X). sults suggest that combining tools may increase the sen- For Teflon, a zygosity filter of 1 resulted in a drop of sitivity of the TIP detection, although this may come sensitivity for MITEs (from 78 to 61.5%) and LTR- with the cost of decreasing precision, as false positives retrotransposons (from 57.7 to 44.2%) but with import- are highly tool-specific. ant gain in precision for LTR-retrotransposons (from 15.2 to 70.8%), which was not significative for MITEs Combining tools to improve sensitivity (98.4 to 98.5%) (not shown). In summary, based on the Our previous results suggest that a combination of tools F1-score, filtering by zygosity greatly improved the over- could be useful to increase the sensitivity in identifying all performance of PopoolationTE2 and Teflon for LTR- non-reference transposon insertions. To this end, we retrotransposon detection, whereas the effect of this fil- combined the predictions of PopoolationTE2 (the overall ter on MITEs detection was much less pronounced due best performer) sequentially with up to four tools se- to the already high precision of the unfiltered results. lected based on their sensitivity and/or precision levels. As a general trend, the combination of tools led to Detection of non-reference insertions by family-specific higher sensitivity levels, reaching more than 90% for tools both MITEs and LTR-retrotransposons at 40X coverage Some tools have been designed to look only for TIPs of when combining five different tools (Fig. 7). However, a single TE family instead of all families at the same time the increase in sensitivity comes with a decrease in pre- (i.e., ITIS, MELT and Trackposon). In order to analyze cision, particularly clear for LTR-retrotransposons, that the performance of such tools, we used the largest MITE approaches 10% for 40X coverage when combining five and LTR-retrotransposon families, which contain 194 different tools. The results presented suggest that the (whole genome) and 22 (chr5) MH63-specific insertions, combination of two tools provided the best balance be- respectively (Additional file 7: Table S1). The analysis of tween sensitivity and precision. Specifically, the combin- MITE TIPs showed that ITIS and MELT did not per- ation of zygosity-filtered PopoolationTE2 and Teflon for form well and displayed low sensitivity and overall F1- MITEs reached 82.1% sensitivity and 97.4% precision at score levels (Fig. 5a). By contrast, Trackposon performed 40X. Regarding LTR-retrotransposons, the combination well, displaying up to 72.8% sensitivity with 93.1 preci- of zygosity-filtered PopoolationTE2 and Jitterbug sion at 40X coverage. In line with the results found for reached 75% sensitivity and 86.7% precision at 40X. broad-spectrum tools, sensitivity in the detection of As already mentioned, McClintock is an available LTR-retrotransposons was strongly dependent on the pipeline that combines several tools. Therefore, we com- coverage. Trackposon and MELT showed moderate sen- pared the performance of the combination of tools here sitivity levels at 40X (58.6 and 55.2%, respectively) proposed with that of the McClintock pipeline, which whereas ITIS reached a maximum of sensitivity of combines the use of Retroseq, TEMP, TE-locate, Popoo- 13.8%. Regarding precision, Trackposon was the best lationTE and ngs_te_mapper (we excluded RelocaTE performer with values ranging between 76.9 and 100% from the pipeline due to excessive running time). The (Fig. 5b). combination of tools here proposed (PopoolationTE2 and Jitterbug for LTR-retrotransposon insertions and Overlap between TIP prediction tools PoPoolationTE2 and Teflon for MITEs) yielded consist- As there is no tool showing 100% sensitivity, we asked ently a better sensitivity and much better precision and whether the predictions of the different tools were com- F1-scores than McClintock at all coverages (especially in mon or specific for each tool. We evaluated the overlap the case of LTR-retrotransposons, Fig. 8). The most Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 10 of 19 Fig. 5 Performance of family-specific tools in the detection of non-reference insertions of MITEs (a) and LTR-retrotransposons (b). Trackposon was run on 10 kb for LTR-retrotransposons windows as described in [7] important differences were found in precision at inter- dataset contained TIPs from 12 different transposon fam- mediate and high coverages. As an example, for MITEs ilies, including retrotransposons (LTR and LINE) and cut- at 40X PopoolationTE2-Teflon had 97.4% precision and-paste DNA transposons (TIR) experimentally vali- whereas McClintock had 83.8% (Fig. 8a). Regarding dated by Lerat et al. [31] Merenciano et al. [33] and Ullas- LTR-retrotransposons at 40X, PoPoolationTE2-Jitterbug tres et al. [46] (Additional file 4: Table S3). The human precision was 86.7%, whereas that of McClintock dataset consisted of 148 TIPs obtained from one human dropped to 9% (Fig. 8b). individual at a coverage of 20X [44]. This dataset consisted of TIPs related to ALU, SVA and LINE-1 retroelements. Evaluation of best-performing tools using Drosophila and In the analysis of human insertions, we also included human datasets MELT, as it is the best-established tool for the detection In order to evaluate whether the benchmarking results of human TE polymorphisms. The detection levels of using rice data could be extrapolated to data obtained PoPoolationTE2 and Teflon in Drosophila were moder- from other species, we benchmarked the best performing ately high (69.1% of the insertions, Table 3 and Additional tools (PoPoolationTE2, Teflon and Jitterbug) using PCR- file 5: Table S4), and substantially higher than Jitterbug validated TIPs from Drosophila and humans. The Dros- (44.4% of the insertions). Using the combination of the ophila dataset consisted of 81 TIPs from ten Drosophila three tools, we were able to detect 81.5% of the insertions. lines sequenced at an average coverage of 42X [22]. This These results are in high concordance with the sensitivity Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 11 of 19 Fig. 6 Venn diagrams representing the detection overlap in non-reference true positives and false positives for MITEs and LTR-retrotransposons levels found using rice data with LTR-retrotransposons Running time and MITEs, where PoPoolationTE2 and Teflon showed Computation time is a limiting factor when running TIP superior detection levels to Jitterbug (Fig. 4). Regarding detection tools in large datasets. Therefore, it is an im- the human sample, MELT was the best tool identifying portant criterion to take into consideration for selecting homozygous insertions (97.8%, Table 4), whereas PoPoo- the most appropriate tool for a specific experiment. We lationTE2 was the best detecting heterozygous insertions tested the time needed by the tools to finish the predic- (88.2%). Taking into account both kind of insertions, tion with a 10X dataset and 432 MITE families as input. PoPoolationTE2 outperformed MELT, displaying an aver- It is important to mention that three tools (Trackposon, age detection level of 90.5%. The detection rate of these ITIS and MELT) work on a per-family basis. In these two programs was higher on human data than in Drosoph- cases, the reported time was that needed to finish the ila or rice, where sensitivity levels rarely exceeded 70% prediction for a single family. By contrast, the remaining using 20X coverage (Fig. 4). The detection levels of Jitter- tools work with all the annotated TE families at the bug were similar to those found using Drosophila and rice, same time. According to our results, Trackposon was ranging from 47.8 to 51%. Teflon was unable to complete the fastest tool, with only 1.7 CPU hours needed to fin- the task and the process was killed after five running days. ish (Fig. 9). Among the general tools, ngs_te_mapper, Using the combination of tools, the detection rate in- TE-locate and PoPoolationTE2 were the fastest tools, creased only 3.4% for the human dataset, reaching up to with 8.6, 9.6 and 9.7 CPU hours needed to finish the 93.9% (Table 4). prediction for the 432 families. RelocaTE2 took the Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 12 of 19 Fig. 7 Performance of tool combinations in the detection of non-reference insertions in MITEs (a) and LTR-retrotransposons (b) largest amount of time to finish the prediction (59.1 transposon insertions, with polymorphic truncated or CPU hours) (Fig. 9). degenerated elements and elements inserted in highly repetitive regions, among other confounding effects, are Discussion also difficult to simulate. As a consequence, the bench- The use of real data is essential for an accurate marking using simulated data may be overestimating the benchmarking of TE insertion detection tools performance of the TIP prediction tools. Indeed, our re- There are several tools available to detect TIPs from sults show that, most of the tools here analyzed have a short-read resequencing data, and some efforts have lower sensitivity than previously reported. For example, been made to validate the performance of such tools [36, RelocaTE2 and TEMP were previously benchmarked on 41]. However, their benchmarking has been essentially simulated rice data, and the sensitivity of both tools was based on simulated TE insertions and simulated short estimated to be higher than 99% at 10X [11]. On the reads. It is challenging to perfectly simulate sequencing contrary, our results using a dataset of real insertions errors, local coverage variations, biases due to GC con- and real short-read data show that both programs per- tent or other genome specific biases that real short-read form very different, with TEMP having a maximum sen- datasets contain. Similarly, the heterogeneity of real sitivity of only 13.3% for MITE detection and Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 13 of 19 Fig. 8 Performance comparison between McClintock pipeline and our proposed tool combinations for MITEs (a) and LTR-retrotransposons (b). PoPoolationTE2 and Teflon are filtered by zygosity as explained in the text (cutoffs of 0.7 and 1, respectively) RelocateTE2 showing a 35.6% sensitivity. Similarly, we that the sensitivity and precision previously reported for previously reported a sensitivity of close to 90% for Jit- TIPs detection tools, determined using simulated data, terbug, a program developed in our laboratory, using are probably overestimated and that the real perform- real short reads on simulated TE insertions [21]. Our re- ance of these tools is probably lower. We think that the sults now show that for the dataset analyzed (real TIPs performance levels of the different tools presented here and real short reads) the maximal sensitivity is of 32.7% are a much better estimation of their detection ability on (Fig. 4, LTR-retrotransposons), although it does so with real datasets. It is important to note, however, that de- a relatively high precision. Therefore, our results suggest pending on the genome to be analyzed, parameters used Table 3 Number of insertions detected by PoPoolationTE2, Jitterbug and Teflon using a validated Drosophila melanogaster dataset RAL-737 RAL-40 RAL-801 RAL-802 RAL-850 RAL-502 RAL-508 RAL-491 RAL-235 RAL-21 TOTAL % Validated insertions 17 16 9 7 4 5 7 5 4 7 81 PoPoolationTE2 12 5 9 5 3 3 6 5 3 5 56 69,1 Jitterbug 11 2 3 5 2 2 4 2 3 2 36 44,4 Teflon 12 6 9 4 3 4 5 4 4 5 56 69,1 Combination 15 6 9 7 3 4 7 5 4 6 66 81,5 Total number of insertions detected by each tool on each line is provided in Additional file 5: Table S4 Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 14 of 19 Table 4 Number of insertions detected by Jitterbug, MELT and their insertions. Because of that, we decided to bench- PoPoolationTE2 using a validated human dataset mark the different programs using two very different Tool Homozygous Heterozygous Total types of TEs that, in addition, are the most prevalent in Validated insertions 46 102 148 plants: MITEs and LTR-retrotransposons. The results PoPoolationTE2 44 (95,7%) 90 (88,2) 134 (90,5%) presented here show that, as expected, the analyzed tools do not detect different TE types with the same sensitivity Jitterbug 22 (47,8%) 52 (51,0%) 74 (50,0%) and precision. MITEs and LTR-retrotransposons repre- Teflon a – – – sent extreme examples based on their length and com- MELT 45 (97,8%) 84 (82,4%) 129 (87,2%) plexity, and the performance of the tools when used Combination 45 (97,8%) 94 (92,2%) 139 (93,9%) with other TEs will probably be in the range of this a Teflon was killed after 5 days running with 12 CPU and 300GB of RAM case-study. The analysis of the sensitivity of the best per- Total number of insertions detected: PoPoolationTE2 (ref and non-ref) = forming tools in detecting TIPs produced by different 186,038; Jitterbug (non-ref) = 624; MELT (non-ref) =1297 types of transposons (including LINEs, LTR- and especially on the quality of the annotation of the ref- retrotransposons and cut-and paste TIR transposons) in erence genome the performance of the programs may Drosophila and humans suggests that this is indeed the vary. All the programs benchmarked here are based on case. Our results indicate that MITEs are detected with the detection of discordant paired-end reads and/or better sensitivity and precision than LTR- split-reads at the junction of TE insertions. Among the retrotransposons. The difference is especially relevant in different confounding factors that can interfere with the the detection of non-reference insertions, where most detection process, the quality of the TE annotation of tools show low precision levels for LTR- the reference genome and in particular of the proper retrotransposons. In the present study, we ran all sam- definition of the TE-genome junctions, is an important ples in default mode or using the parameters described one. Therefore, it is important to work on refining the by the authors in the corresponding manuscripts or annotation of the TEs (or at least the more interesting manuals (Additional file 1). Nevertheless, we show that TE families for the purpose of the study) before search- precision can be increased by applying specific filters to ing for TIPs. the results. For example, we show that, for some pro- grams, LTR-retrotransposon detection can be drastically Tool performance varies depending on TE family improved by applying a zygosity filtering. Applying such Eukaryote genomes contain a high diversity of TE ele- filtering may be a good strategy when not intending to ments with very different copy numbers and functional study somatic insertions which should in most cases be and structural characteristics, which may impact on the heterozygous. The difficulties of detecting LTR- ability of TIP detecting programs to reliably identify retrotransposons come from the complexity of the Fig. 9 Running time of each tool to perform the detection of MITEs in a 10X dataset. Family-specific tools are marked with an asterisk. All tools were run using 8 CPUs and 64GB of RAM Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 15 of 19 elements and from the local regions where they insert. It performance of the tools is relatively independent of the is known that LTR-retrotransposons (especially those of genome used, and that the benchmarking here presented the Gypsy superfamily) tend to integrate in heterochro- could be of use for TIP analysis in many different matic regions enriched in other TEs. These repetitive re- systems. gions are likely a source of false positives that affects all This analysis also showed that the tools that per- the programs tested. These repetitive regions are, in fact, formed best on rice had an even better sensitivity on hu- difficult to annotate and polymorphisms within these re- man data. The difference of sensitivity was particularly gions may be challenging to detect even using long-read clear for PoPoolationTE2 and MELT. Although this data or when aligning good-quality assemblies. By con- could indicate a difference of the performance of these trast, MITEs tend to integrate close to genes [25] and tools in the two genomes, it could also be due to the their flanking regions are more likely to be unique in the particular nature of the human dataset. Indeed, the data- genome. The presence of non-repetitive TE flanks set of validated TIPs in humans contains insertions from greatly simplifies the detection of TIPs, as the probability TE families (LINE-1, ALU, SVA) that were detected in of finding multimapping reads in these regions is the first place using only one method, based on split- minimal. read and read-pair information [44] and therefore the Another important consideration related to the differ- sensitivity of the programs on this dataset could be over- ent TE families is the quality of the annotation. MITEs estimated. It is worth mentioning that the PCR-validated are easy to annotate and usually have well defined Drosophila and human insertions have been predicted boundaries. By contrast, LTR-retrotransposons form using a small number of tools in the original publica- nested insertions and are often degenerated. This makes tions, and therefore it includes only a subset of all the very difficult to accurately define their boundaries, and insertions present in these genomes. Moreover, the hu- as a consequence many chimeric elements are usually man and Drosophila datasets were validated by PCR, annotated. As already mentioned, an accurate TE anno- which could have introduced a bias in the TEs that were tation is essential to increase the capacity of the tools to included in these datasets. However, note that the num- identify TE insertions based on short-read data. In this ber of families included in the human and Drosophila context, it could be a good strategy to identify and re- validation datasets are similar or bigger than the ones in- move chimeric transposons from the annotation prior to cluded in the rice dataset and contain both full-length using any of these tools (ie, when working with consen- and truncated TEs. suses or with the actual annotation). A chimeric or nested transposon, for example an LTR-retrotransposon Sequencing coverage critically impacts TIP detection with a MITE inserted inside, will be targeted by reads Independently of the different performance found be- arising from the two elements, and other MITE inser- tween TE families, we found that coverage has a major tions of the same family present elsewhere in the gen- impact on tool performance for all the TE families ome could be wrongly identified as LTR- tested. In general sensitivity increases with increasing retrotransposons insertions by the TIP detection tools. coverage. Therefore, homogenization of sample coverage is essential when using TIPs prediction tools to quantita- Influence of the type of genome on the performance of tively compare the transposition rates between organ- the tools isms or populations. Some tools like PopoolationTE2 The ability of any of the tools to detect TIPs depend on have internal steps to carry out this task. Nevertheless, the nature of the transposon insertion itself and its for qualitative studies coverage homogenization is dis- flanking genome sequence, and none of them can detect couraged as down-sampling high-coverage datasets leads new transposon insertions in repetitive regions. There- to a smaller number of detected insertions. It is import- fore, in addition to the type of transposon generating the ant to note that the increase of sensitivity with increas- TIP, as already discussed, the performance of the tools ing coverage comes, in most cases, with a decrease in may depend on the genome under study. For this reason, precision. Therefore, depending on the goals of the study we have analyzed the sensitivity of the tools that per- a different level of coverage may be suitable. From the formed the best using rice data on Drosophila and hu- data presented here it seems that a coverage below 20X man data and compared their performance on the is probably not suited for most analyses, as the probabil- different datasets. The sensitivity of the different pro- ity of missing true insertions is very high. grams analyzed in Drosophila was very similar to the one obtained in rice. As the genomes of rice and Dros- Strategies to increase tool performance ophila are relatively different, the former being much The fact that an important fraction of the insertions de- bigger (430 Mb vs 175 Mb) and with a higher content of tected by the different tools are not shared supports the repetitive sequences (37% vs 20%), this suggests that the fact that combining different tools may increase the Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 16 of 19 quality of the results [36]. However, simply increasing relatively fast and show high sensitivity and precision in- the number of tools does not necessarily increase the dependently of the species and TE type analyzed. quality of predictions, due to the accumulation of tool- specific false positives (ie, the combination of five tools Conclusions yielded 95% of sensitivity but only 11.8% precision in Besides the important efforts of tool developers, our re- non-ref LTR-retrotransposon detection, Fig. 7). This is sults suggest that the identification of TIPs is still chal- due to the fact that whereas many true insertions are de- lenging. We propose here a number of approaches, such tected by several tools, most false positives are tool- as combining tools, which can be followed depending on specific (Fig. 6). Combining a limited number of well- the purpose of the study and the TE families to be ana- performing tools may be the best approach. Indeed, our lyzed, that can provide good results. However, it is im- results show that with the dataset used, the combination portant to note that in the best scenario (combining of PoPoolationTE2 and Jitterbug to detect LTR- optimal tools at best coverage, Fig. 7) and having a good retrotransposon insertions, or PoPoolationTE2 and Tef- TE annotation of the reference genome, the sensitivity lon to detect MITEs yielded superior TIP annotations could be around 70% with a precision of 80–90% for (better F1-score) than the tools alone. Also, the perform- non-reference insertions. These numbers may be enough ance of these tool combinations was better than that of for most studies, but it is important to keep in mind that the McClintock pipeline, especially regarding LTR- some insertions will be missed, especially when estimat- retrotransposons. In this sense, we recommend combin- ing insertion frequencies or when using TIPs for GWAS, ing tools based on their high precision and not only on for example. There are major limitations like the length their high sensitivity (ie, PoPoolationTE2 and Jitterbug). of the reads that can be resolved with current technolo- Nevertheless, there can be situations in which sensitivity gies (ie long-read sequencing) and will certainly improve has a priority over precision (ie, re-sequencing of a sin- in the following years. But there is still the need to de- gle individual, or interest only on a few families). In such velop new algorithms specifically designed to identify cases, running more tools can be an alternative and TIPs from long reads, to generate highly curated TE an- manual curation should be considered. notations of reference genomes and also more independ- ent benchmarks on real data to evaluate the performance of tools under different conditions. Selecting the appropriate tools for detecting TE insertions in resequencing data Methods Depending on the objective of the analysis, a family- Sequence data used specific tool could be more interesting than a broad- We used the available data for the japonica Nipponbare spectrum tool. For example, when tracking the effect of (GCA_000005425.2) and the indica MH63 (GCA_ certain treatment in a concrete set of elements. Another 001623365.1) assemblies, and the short-read resequen- important consideration is that the amount of storage cing of MH63 (SRX1639978), which were used to gener- needed is smaller in comparison to broad-spectrum ate the original assembly. tools, due to the smaller size of the alignment files. For such cases, a tool such as Trackposon could be a good MITE annotation option due to its fast speed, moderate sensitivity and MITE-hunter [19] was run on Nipponbare and MH63 high precision. Nevertheless, as a drawback, Trackposon assemblies to detect MITEs families, which were then does not report the exact insertion point and, which combined with the high-quality predictions available in could be a limitation for some studies. In those cases, PMITE database [9] (only families carrying TSD). Clus- MELT can be an interesting alternative, although it re- tering at 90% was performed to remove redundancy quires adjusting family-specific parameters to produce using cd-hit [17] and produce a final library. RepeatMas- high-quality results. This might be indeed the cause why ker (http://www.repeatmasker.org/) was run to annotate MELT did not perform well on the detection of rice all regions having significant homology with any of the MITEs. In general, it is possible that the tools analyzed MITE families. The annotations were further screened here, which were not specifically designed for MITEs to discriminate full-length elements (consensus length ± and LTR-retrotransposons, may work better for other 20%) from truncated hits. types of TEs or with modifications in the parameters used. Based on our results, if the objective of the study is LTR-retrotransposon annotation to analyze insertions of more than one family, and the LTR-retrotransposons were identified by running storage space is not a major limitation, using some of LTRharvest [14] on IRGSP and MH63 assemblies with the top broad-spectrum tools such as PoPoolationTE2 is default parameters. The internal conserved domains of probably a better option as those programs can also be these elements were obtained running hmmscan [13], Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 17 of 19 and only coding elements were retained for further ana- retrotransposons: True positives (TP): Insertions de- lyses. The identified elements were clustered with Silix tected by any tool matching with our curated dataset of [34] according to the 80–80 rule. All the elements in TPs. False positives (FP): Insertions detected by any tool each family were aligned with Mafft [26] and trimmed matching with our curated dataset of FPs. False nega- with Trimal [6]. Consensus sequences were built from tives (FN): Insertions present in our curated dataset of the alignments using the EMBOSS package [40]. TPs, not detected by the evaluated tool. These primary parameters were used for calculating the final bench- Determination of benchmarking standards marking ratios that have been previously used for asses- We took advantage of the availability of two high quality sing the performance of similar tools [41]. rice genome assemblies (IRGSP and MH63, the former used as reference) in order to obtain a curated dataset of Sensitivity = TP/ (TP+ FN). real “reference” (orthologous) and “non-reference” (spe- Precision = TP/ (TP + FP) cific to MH63) insertions as explained in Fig. 2. Mapping F1-score = 2 x [(Precision x Sensitivity) / (Precision of reference and non-reference windows to MH63 gen- + Sensitivity)] ome was performed using BBmap (https://sourceforge. net/projects/bbmap/). Intersections between annotations Supplementary information were done with BEDtools [38]. Supplementary information accompanies this paper at https://doi.org/10. 1186/s13100-019-0197-9. Drosophila and human benchmarking datasets The Drosophila dataset consisted of 81 TIPs from ten Additional file 1. Scripts used to run all TIP detection tools. (.sh) (SH 13 kb) Drosophila lines sequenced at an average coverage of Additional file 2 : Table S1. Insertion dataset used for benchmarking. Contains all the TP reference and non-reference windows, as well as all 42X [22], and validated by PCR by Lerat et al. [31], Mer- FP windows. (.xlsx) enciano et al. [33] and Ullastres et al. [46] (Additional Additional file 3 : Table S2. Numerical benchmark results. (.xlsx) file 4: Table S3). In Lerat et al. [31], TIPs were predicted Additional file 4 : Table S3. Drosophila melanogaster TE insertions using TIDAL [39] and PoPoolationTE2 [29] using 14 validated by PCR. TE names are provided in the table when authors gave European D. melanogaster pooled populations (average names to the non-reference insertions. All insertions were validated based on PCR band sizes. Validated insertion sites are provided when au- coverage of 90X). Briefly, validated TIPs were present in thors sequenced PCR bands evidencing the presence of a particular inser- the DGRP population and at least in one European tion. ND: not determined. (xlsx) population at > 10% frequency, not present in the Y Additional file 5 : Table S4. Total number of insertions detected by chromosome, and with a predicted length of < 6 kb to PoPoolationTE2, Jitterbug and Teflon in ten Drosophila lines. (xlsx) avoid problems with PCR amplification. In Ullastres Additional file 6 : Figure S1. Number of MH63 reference and non- reference insertions detected by direct comparison of 1000 LTR- et al. [46], TIPs were predicted by TIDAL in the DGRP retrotransposon flanking sites of different sizes from MH63 and Nippon- population [39]. Validated TIPs were inserted in regions bare genomes. (.pdf) with recombination rates > 0, and present in at least 15 Additional file 7 : Figure S2. Application of zygosity filtering to DGRP strains. Finally, in Merenciano et al. [33] TIPs PoPoolationTE2. PoPoolationTE2-F means that it was run and filtered at zygosity 0.7. PopoolationTE2-R corresponds to the raw results. (.png) were also predicted by TIDAL in the DGRP population [39] and all belonged to the roo family. Both full-length and truncated copies were validated, as no TE length fil- Acknowledgements ter was applied. Not applicable. The human dataset consisted of 148 TIPs obtained from a human individual (NA12891, SRA accession Authors’ contributions FB ran Jitterbug, MELT and produced the LTR-retrotransposon annotation. PV SRX207113) [44]. Original sequencing coverage of the ran all TIP detection tools and analyzed the data. RC produced the MITE an- human genome was down sampled to 20X. notation and the final dataset of TE insertions and analyzed the data to- gether with PV. MM and JG participated in the benchmarking using Drosophila and human datasets. JMC and RC conceived the study and wrote TIP prediction the manuscript. All the authors reviewed the final manuscript. All authors Predictions of transposon insertions were done using the read and approved the final manuscript. 12 tools shown in Table 2 using the default parameters and / or following the recommendations of the authors. Funding The scripts used for running each of the tools are shown This work was supported in part by grants from the Ministerio de Economia y Competitividad (AGL2016–78992-R). Fabio Barteri and Pol Vendrell hold a in Additional file 1. FPI (Formación de Personal Investigador) fellowship from the Spanish Ministerio de Economia y Competitividad. Raúl Castanera holds a Juan de la Evaluation parameters Cierva Postdoctoral fellowship from the Spanish Ministerio de Economia y Competitividad. JG is funded by the European Commission (H2020-ERC- We used the following parameters for evaluating the 2014-CoG-647900) and the Spanish Ministerio de Ciencia, Innovación y ability of each tool to detect MITEs and LTR- Universidades/AEI/FEDER, EU (BFU2017–82937-P). Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 18 of 19 Availability of data and materials 16. Flutre T, Duprat E, Feuillet C, Quesneville H. Considering transposable The datasets analyzed during the current study are available in the NCBI element diversification in de novo annotation approaches. Plos One. 2011; repository: 6(1):e16526. - Nipponbare Assembly: GCA_000005425.2 17. Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next- - MH63 assembly: GCA_001623365.1 generation sequencing data. Bioinformatics. 2012;28(23):3150–2. - Short-read resequencing data of MH63: SRX1639978 18. Gardner EJ, Lam VK, Harris DN, et al. The Mobile element locator tool - Human resequencing reads: SRX207113 (MELT): population-scale mobile element discovery and biology. Genome - Drosophila resequencing reads: PRJNA36679 Res. 2017;27(11):1916–29. 19. Han Y, Wessler SR. MITE-hunter: a program for discovering miniature inverted-repeat transposable elements from genomic sequences. Nucleic Ethics approval and consent to participate Acids Res. 2010;38(22):e199. Not applicable. 20. Hénaff E, Vives C, Desvoyes B, et al. Extensive amplification of the E2F transcription factor binding sites by transposons during evolution of Consent for publication Brassica species. Plant J. 2014;77(6):852–62. Not applicable. 21. Hénaff E, Zapata L, Casacuberta JM, Ossowski S. Jitterbug: somatic and germline transposon insertion detection at single-nucleotide resolution. BMC Genomics. 2015;16:768. Competing interests 22. Huang W, Massouras A, Inoue Y, et al. Natural variation in genome The authors declare that they have no competing interests. architecture among 205 Drosophila melanogaster genetic reference panel lines. Genome Res. 2014;24(7):1193–208. Author details 23. International Rice Genome Sequencing Project. The map-based sequence of 1 Centre for Research in Agricultural Genomics CSIC-IRTA-UAB-UB, Campus the rice genome. Nature. 2005;436(7052):793–800. UAB, Edifici CRAG, Bellaterra, 08193 Barcelona, Spain. 2Institute of 24. Jiang C, Chen C, Huang Z, Liu R, Verdier J. ITIS, a bioinformatics tool for Evolutionary Biology (CSIC-Universitat Pompeu Fabra), Passeig Maritim accurate identification of transposon insertion sites using next-generation Barceloneta 37-49, 08003 Barcelona, Spain. sequencing data. BMC Bioinformatics. 2015;16(1):72. 25. Jiang N, Wessler SR. Insertion preference of maize and rice miniature Received: 28 June 2019 Accepted: 17 December 2019 inverted repeat transposable elements as revealed by the analysis of nested elements. Plant Cell. 2001;13(11):2553–64. 26. Katoh K, Standley DM. MAFFT multiple sequence alignment software References version 7: improvements in performance and usability. Mol Biol Evol. 2013; 1. Adrion JR, Song MJ, Schrider DR, Hahn MW, Schaack S. Genome-wide 30(4):772–80. estimates of transposable element insertion and deletion rates in 27. Keane TM, Wong K, Adams DJ. RetroSeq: transposable element discovery Drosophila melanogaster. Genome Biol Evol. 2017;9(5):1329–40. from next-generation sequencing data. Bioinformatics. 2013;29(3):389–90. 2. Alkan C, Coe BP, Eichler EE. Genome structural variation discovery and 28. Kobayashi S, Goto-Yamamoto N, Hirochika H. Retrotransposon-induced genotyping. Nat Rev Genet. 2011;12(5):363–76. mutations in grape skin color. Science. 2004;304(5673):982. 3. Bogaerts-Márquez M, Barrón MG, Fiston-Lavier A-S, et al. T-lex3: an accurate 29. Kofler R, Gómez-Sánchez D, Schlötterer C. PoPoolationTE2: comparative tool to genotype and estimate population frequencies of transposable population genomics of transposable elements using Pool-Seq. Mol Biol elements using the latest short-read whole genome sequencing data. Evol. 2016;33(10):2759–64. Bioinformatics. 2019, btz727. 30. Kofler R, Orozco-terWengel P, De Maio N, et al. PoPoolation: a toolbox for 4. Butelli E, Licciardello C, Zhang Y, et al. Retrotransposons control fruit- population genetic analysis of next generation sequencing data from specific, cold-dependent accumulation of anthocyanins in blood oranges. pooled individuals. Plos One. 2011;6(1):e15925. Plant Cell. 2012;24(3):1242–55. 31. Lerat E, Goubert C, Guirao-Rico S, et al. Population-specific dynamics and 5. Cao Y, Chen G, Wu G, et al. Widespread roles of enhancer-like transposable selection patterns of transposable element insertions in European natural elements in cell identity and long-range genomic interactions. Genome Res. populations. Mol Ecol. 2019;28(6):1506–22. 2019;29(1):40–52. 32. Linheiro RS, Bergman CM. Whole genome resequencing reveals natural 6. Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. trimAl: a tool for target site preferences of transposable elements in Drosophila automated alignment trimming in large-scale phylogenetic analyses. melanogaster. Plos One. 2012;7(2):e30008. Bioinformatics. 2009;25(15):1972–3. 33. Merenciano M, Iacometti C, González J. A unique cluster of roo insertions in 7. Carpentier M-C, Manfroi E, Wei F-J, et al. Retrotranspositional landscape of the promoter region of a stress response gene in Drosophila melanogaster. Asian rice revealed by 3000 genomes4. Nat Commun. 2019;10(1):2. Mob DNA. 2019;10:10. 8. Carr M, Bensasson D, Bergman CM. Evolutionary genomics of transposable 34. Miele V, Penel S, Duret L. Ultra-fast sequence clustering from similarity elements in Saccharomyces cerevisiae. Plos One. 2012;7(11):e50978. networks with SiLiX. BMC Bioinformatics. 2011;12:116. 9. Chen J, Hu Q, Zhang Y, Lu C, Kuang H. P-MITE: a database for plant 35. Naito K, Zhang F, Tsukiyama T, et al. Unexpected consequences of a sudden miniature inverted-repeat transposable elements. Nucleic Acids Res. 2014; and massive transposon amplification on rice gene expression. Nature. 42(Database issue):D1176–81. 2009;461(7267):1130–4. 10. Chen J, Lu L, Benjamin J, et al. Tracking the origin of two genetic 36. Nelson MG, Linheiro RS, Bergman CM. McClintock: an integrated pipeline for components associated with transposable element bursts in domesticated detecting transposable element insertions in whole-genome shotgun rice. Nat Commun. 2019;10(1):641. sequencing data. G3. 2017;7(8):2763–78. 11. Chen J, Wrightsman TR, Wessler SR, Stajich JE. RelocaTE2: a high resolution 37. Platzer A, Nizhynska V, Long Q. TE-locate: a tool to locate and group transposable element insertion site mapping tool for population transposable element occurrences using paired-end next-generation resequencing. PeerJ. 2017;5:e2942. sequencing data. Biology. 2012;1(2):395–410. 12. Chuong EB, Elde NC, Feschotte C. Regulatory activities of transposable 38. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing elements: from conflicts to benefits. Nature Reviews. Genetics. 2017;18(2): genomic features. Bioinformatics. 2010;26(6):841–2. 71–86. 39. Rahman R, Chirn G, Kanodia A, et al. Unique transposon landscapes are 13. Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol. 2011;7(10): pervasive across Drosophila melanogaster genomes. Nucleic Acids Res. e1002195. 2015;43(22):10655–72. 14. Ellinghaus D, Kurtz S, Willhoeft U. LTRharvest, an efficient and flexible 40. Rice P, Longden I, Bleasby A. EMBOSS: the european molecular biology software for de novo detection of LTR retrotransposons. BMC open software suite. Trends Genet. 2000;16(6):276–7. Bioinformatics. 2008;9:18. 41. Rishishwar L, Mariño-Ramírez L, Jordan IK. Benchmarking computational 15. Ewing AD. Transposable element detection from whole genome sequence tools for polymorphic transposable element detection. Brief Bioinform. 2017; data. Mob DNA. 2015;6:24. 18(6):908–18. Vendrell-Mir et al. Mobile DNA (2019) 10:53 Page 19 of 19 42. Sanseverino W, Hénaff E, Vives C, et al. Transposon insertions, structural variations, and snps contribute to the evolution of the melon genome. Mol Biol Evol. 2015;32(10):2760–74. 43. Slotkin RK, Martienssen R. Transposable elements and the epigenetic regulation of the genome. Nat Rev Genet. 2007;8(4):272–85. 44. Stewart C, Kural D, Strömberg MP, et al. A comprehensive map of mobile element insertion polymorphisms in humans. PLoS Genet. 2011;7(8): e1002236. 45. Studer A, Zhao Q, Ross-Ibarra J, Doebley J. Identification of a functional transposon insertion in the maize domestication gene tb1. Nat Genet. 2011; 43(11):1160–3. 46. Ullastres A, Merenciano M, González J. Natural transposable element insertions drive expression changes in genes underlying Drosophila immune response. BioRxiv. 2019, 655225. 47. Zhang J, Chen L-L, Sun S, et al. Building two indica rice reference genomes with PacBio long-read and Illumina paired-end sequencing data. Sci Data. 2016;3:160076. 48. Zhuang J, Wang J, Theurkauf W, Weng Z. TEMP: a computational method for analyzing transposable element polymorphism in populations. Nucleic Acids Res. 2014;42(11):6826–38. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.