Guillaume Rizk
email: guillaume.rizk@inria.fr

Anaïs Gouin

Rayan Chikhi

Claire Lemaitre
email: claire.lemaitre@inria.fr

MindTheGap : integrated detection and assembly of short and long insertions

Motivation: Insertions play an important role in genome evolution. However, such variants are difficult to detect from short read sequencing data, especially when they exceed the paired-end insert size. Many approaches have been proposed to call short insertion variants based on paired-end mapping. However, there remains a lack of practical methods to detect and assemble long variants. Results: We propose here an original method, called MINDTHEGAP, for the integrated detection and assembly of insertion variants from re-sequencing data. Importantly, it is designed to call insertions of any size, whether they are novel or duplicated, homozygous or heterozygous in the donor genome. MINDTHEGAP uses an efficient k-mer based method to detect insertion sites in a reference genome, and subsequently assemble them from the donor reads. MINDTHEGAP showed high recall and precision on simulated datasets of various genome complexities. When applied to real C. elegans and human NA12878 datasets, MINDTHEGAP detected and correctly assembled insertions longer than 1 kb, using at most 14 GB of memory.

INTRODUCTION

Structural variants (SVs) are large-scale structural changes in the genome. They have been typically defined in opposition to point mutations, which are single nucleotide polymorphisms (SNPs) and short insertions or deletions (indels). SVs therefore include insertions, deletions and inversions of genomic sequences. Recent research has shown that they play an important role in evolution and diseases [START_REF] Lemaitre | A map of human genome variation from population-scale sequencing[END_REF][START_REF] Stewart | A comprehensive map of mobile element insertion polymorphisms in humans[END_REF]. However, SVs are challenging to discover using present-day sequencing approaches, as they generally span genomic regions that are longer than the reads. Computational methods have been designed to extract evidence of SVs from sequencing data using two types of analyses: paired-end mapping of reads to a reference genome and copy-number estimation using read depth [START_REF] Medvedev | Computational methods for discovering structural variation with next-generation sequencing[END_REF][START_REF] Alkan | Genome structural variation discovery and genotyping[END_REF]. * to whom correspondence should be addressed

Definition of insertion variants

In this work we will focus on insertion variants: sequences that are present at one site (position) in the donor genome but are absent from the reference genome at this site. We divide insertions into three mutually exclusive types: (i) novel insertions in the donor genome that have no match in the reference, (ii) duplicated insertions, that are found at two or more sites in the donor and a strict subset of those in the reference, and (iii) transpositions, that are sequences in the reference that moved to a different site in the donor. Duplicated insertions include mobile element insertions (MEI), for which databases of known sequences have been created to facilitate discovery [START_REF] Stewart | A comprehensive map of mobile element insertion polymorphisms in humans[END_REF].

All three types of insertions are difficult to detect using short reads. Different techniques are used to detect insertions that are short (shorter than the reads), medium (of size between read length and insert size) or long (of size exceeding insert size). In the next two sections, we review techniques used to identify insertion sites, and techniques used to reconstruct insertion sequences.

Identification of insertion sites

As short insertions are likely to be fully contained in several reads, mapping donor reads to a reference genome enables simultaneous discovery of the sites and contents of insertions [START_REF] Depristo | A framework for variation discovery and genotyping using next-generation dna sequencing data[END_REF][START_REF] Li | The sequence alignment/map format and samtools[END_REF][START_REF] Ye | Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads[END_REF][START_REF] Albers | Dindel: Accurate indel calls from short-read data[END_REF]. In this context, results are sensitive to mapping parameters and may be degraded in low-coverage or low-complexity regions of the reference. Although the discovery of short indels has been an extensively studied problem, a recent article has observed considerable differences between the results of popular tools [START_REF] Pabinger | A survey of tools for variant analysis of next-generation genome sequencing data[END_REF].

Sites of medium-sized insertions can be detected by analyzing mapping positions of paired reads. General SV calling tools call insertions sites by clustering neighboring read pairs that have a shorter insert size than expected, e.g. BreakDancer and GASV [START_REF] Chen | Breakdancer: an algorithm for highresolution mapping of genomic structural variation[END_REF][START_REF] Sindi | A geometric approach for classification and comparison of structural variants[END_REF]. NovelSeq [START_REF] Hajirasouliha | Detection and characterization of novel sequence insertions using paired-end next-generation sequencing[END_REF] and SOAPindel [START_REF] Li | Soapindel: efficient identification of indels from short paired reads[END_REF] detect sites of long, novel insertions by clustering paired reads for which one mate is unmapped.

Alternatively, tools based on read coverage can detect duplicated insertions of any length by finding reference segments that have higher read depth than expected. While insertion sites cannot be determined by this method alone, the Reprever [START_REF] Kim | Reprever: resolving lowcopy duplicated sequences using template driven assembly[END_REF] software identifies low-copy duplicated insertions by combining paired-end c Oxford University Press 2014.

1

mapping with read depth analysis. Finally, several methods detect sites of mobile element insertions using collections of known transposable element sequences, by searching for read pairs where one mate is mapped to a known element and the other to a unique part of the reference genome [START_REF] Hormozdiari | Next-generation variationhunter: combinatorial algorithms for transposon insertion discovery[END_REF][START_REF] Ewing | Whole-genome resequencing allows detection of many rare line-1 insertion alleles in humans[END_REF][START_REF] Stewart | A comprehensive map of mobile element insertion polymorphisms in humans[END_REF].

Reconstruction of inserted sequences

While short insertions are easy to reconstruct (as seen in Section 1.2), to the best of our knowledge, only a few methods are capable of handling medium or long insertions. They are based on global or local de novo assembly of reads that are potentially involved in an insertion.

SOAPindel [START_REF] Li | Soapindel: efficient identification of indels from short paired reads[END_REF], Scalpel [START_REF] Narzisi | Accurate detection of de novo and transmitted indels within exome-capture data using micro-assembly[END_REF] and TIGRA [START_REF] Chen | Tigra: A targeted iterative graph routing assembler for breakpoint assembly[END_REF] select paired reads for which one of the mates maps nearby an insertion site. The other mates are used to assemble separately each inserted sequence. This approach can only reconstruct insertions that are shorter than twice the insert size. NovelSeq [START_REF] Hajirasouliha | Detection and characterization of novel sequence insertions using paired-end next-generation sequencing[END_REF] reconstructs novel insertions (of any size) by assembling all unmapped reads, then aligning the extremities of assembled sequences to all predicted insertion sites. Parrish et al proposed to extend this approach to duplicated insertions by performing a global assembly of all reads that are either unmapped, discordantly paired or mapped to high-coverage regions [START_REF] Parrish | Assembly of non-unique insertion content using next-generation sequencing[END_REF].

Cortex var [START_REF] Iqbal | De novo assembly and genotyping of variants using colored de bruijn graphs[END_REF] builds a colored de Bruijn graph from the reference genome and all donor reads. Insertions appear in the graph as bubbles (sets of paths between two nodes), where one short path corresponds to the reference genome, and longer paths correspond to inserted sequences. Theoretically, this approach enables the discovery of insertions regardless of their size and type. However, due to practical limitations, Cortex var only finds a restricted class of bubbles: those that (i) contain exactly two paths and (ii) all intermediate nodes have exactly one in-neighbor and one out-neighbor.

To summarize, available tools are highly specialized and lack the versatility to detect and assemble insertions of any size and any type. SOAPindel, Scalpel and Cortex var are practically limited to short insertions, Reprever is limited to low-copy duplicated sequences and Novelseq is limited to novel insertions.

Our contribution

We propose a new tool, MINDTHEGAP, for detecting and assembling insertions. MINDTHEGAP has several novel features that are not found in other tools. First, an mapping-free site detection algorithm has been designed to detect insertions of any size. Second, an improved method for insertion assembly enables the reconstruction of long insertions of all three types. Third, a memory-efficient data structure enables high scalability.

We evaluated MINDTHEGAP on simulated and real Illumina sequencing data. Among 1 kbp simulated homozygous insertions, a large fraction were found and correctly assembled (recall values between 65-98.4% precision over 97%). Simulated heterozygous 1 kbp insertions proved to be more challenging to assemble (60% recall for C. elegans, 35% for human chromosome 22), however precision remained high (resp. 93% and 89%). We assembled long insertions using MINDTHEGAP on an actual whole-genome human dataset, which required only 14 GB of memory.

METHODS

The input of MINDTHEGAP is a set of reads and a reference genome. The software performs three steps: (1) construction of the de Bruijn graph of the reads, (2) detection of insertion breakpoints on the reference genome (find module) and (3) local assembly of inserted sequences (fill module). Both the detection step and the assembly step rely solely on the constructed graph.

The output of the second step is a set of putative insertion positions on the reference genome, whereas the output of the last step is, for each insertion site, one or several assembled sequences.

de Bruijn graph construction

The de Bruijn graph is a directed graph over all distinct k-mers in the reads. An edge is present when two k-mers share an exact (k -1)-overlap. The graph is constructed using the algorithms implemented in the Minia assembler [START_REF] Chikhi | Space-efficient and exact de bruijn graph representation based on a bloom filter[END_REF][START_REF] Salikhov | Using cascading bloom filters to improve the memory usage for de brujin graphs[END_REF]. Minia encodes the graph using a Bloom filter and an additional hash table to suppress false positives. The data structure supports two operations: (i) membership queries for k-mers that are neighbors of existing k-mers in the graph, and (ii) traversal of the graph from an existing k-mer. These operations are respectively used in Section 2.2 (insertion site detection) and Section 2.3 (local assembly).

Find module : detection of insertion sites

MINDTHEGAP detects insertion sites by scanning the reference genome and testing membership of reference k-mers in the de Bruijn graph. Homozygous and heterozygous insertions are handled using two different methods.

Homozygous insertions

The general case for detecting homozygous insertions can be modeled as follows. Let Sr be a sequence (the reference). For a position j in the reference, the k-mer at position (jk + 1) (resp. j + 1) is called the left (resp. right) flanking k-mer. Let S d (the donor) be a copy of Sr where a sequence I has been inserted between the nucleotides at position i and i + 1 (the insertion site, see Figure 1). For each position in the reference genome, a binary character records whether the k-mer starting at this position is present ('1') or absent ('0') in the donor reads. Depending on the context, the reference genome will correspond to the string of nucleotides or to the string of binary characters. Let a gap be a substring in the reference genome equal to 0 n (formed by repeating '0' n times), for n > 0, that is immediately flanked by '1' characters. In most cases, a homozygous insertion site at position i has a gap of size k -1 starting at position i-k +2 (all k -1 k-mers overlapping the insertion site are absent in S d ). We refer to this situation as a canonical insertion site (see Figure 1-A).

A gap may be shorter than k -1, for instance when the prefix of the inserted sequence I exactly matches the prefix of the sequence to the right of the insertion site (see Figure 1-B). More generally, if the longest common prefix of I and Sr[i + 1 . . .] is of size r 1 and the longest common suffix of I and Sr[. . . i] is of size r 2 , then the size of the gap is kr 1 -r 2 -1. It is important to note that when r 1 or r 2 are greater than 0, with only sequences S d and Sr at hand, it is not possible to localize precisely the insertion site, as it can be at any location in [ir 2 . . . i + r 1 ]. We refer to such sites as fuzzy sites. Homozygous insertion sites are called when gaps of size in the range [k -1 -r, k -1] are detected, with r being a user-defined parameter indicating the largest allowed repeat at the insertion.

The size of the gap is an important criterion to detect homozygous insertion sites, since other types of variants also yield gaps. SNPs create gaps of size exactly k, deletions of length d yield gaps of size k + d -1.

Variants that are separated by less than k nucleotides yield longer gaps. In fact, only new junctions between existing sequences can yield gaps of size < k, which is the case for insertion events, but also for inversion or translocation sites. Finally, gaps of various sizes may also appear due to insufficient read coverage or non-uniqueness of k-mers inside the reference genome. These effects are controlled by the value of k, which is a parameter of our method. be detected. The left flanking k-mer of a heterozygous insertion site has at least two out-neighbors in the de Bruijn graph: one neighbor in the reference sequence and at least one other neighbor that is a prefix of the inserted sequence. Similarly, the right flanking k-mer has at least two in-neighbors with similar properties (see Figure 1-C). As in the homozygous case, small repetitions at the extremities of inserted sequences slightly alter the pattern. The left flanking k-mer may overlap the right flanking k-mer in the reference genome. MINDTHEGAP detects heterozygous insertion sites by scanning the reference genome and testing neighborhoods of putative left and right flanking k-mers whose distance from one another is comprised between kr and k, r being the same user-defined parameter as for homozygous insertions, indicating the largest allowed repeat at the insertion.

Heterozygous insertions

Note that heterozygous SNPs and deletions yield similar patterns, but the left and right flanking k-mers are further separated from each other (k + 1 nucleotides apart for SNPs and 1-bp deletions, k+d-1 nucleotides apart for deletions of size d). However, heterozygous inversions and translocations do exhibit identical patterns. Also, inexact repetitions in the reference genome create branching k-mers, that may yield by chance the same pattern as a heterozygous insertion. To reduce this effect, we apply an additional filter: L .R.

.R. the k -1 suffix (resp. prefix) of the right (resp. left) flanking kmer must have less than h occurrences in the reference genome. When h is set to 1, this prevents the detection of patterns that may be generated by repetitions in the reference genome alone, in absence of any sequence variants.

Fill module : assembly of inserted sequences

The third step of MINDTHEGAP is called the fill module. Starting from a known insertion site represented by flanking kmers (L, R), the module performs de novo assembly to attempt to reconstruct the inserted sequence between L and R. In a nutshell, a graph of contigs is constructed by performing breadth-first traversal of k-mers, starting from L. The traversal is halted when graph becomes too complex. Then, all the contigs in the graph are searched for the presence of R. All paths between L and the contigs containing R are enumerated, and one or more putative insertion sequences are returned. More specifically, insertions are assembled using the algorithm of Minia [START_REF] Chikhi | Space-efficient and exact de bruijn graph representation based on a bloom filter[END_REF][START_REF] Salikhov | Using cascading bloom filters to improve the memory usage for de brujin graphs[END_REF]. Assembly is performed by traversing the graph from a given starting k-mer in a breadth-first fashion. A consensus sequence (contig) is generated by skipping over certain motifs, such as bubbles (putative short variants) and tips (putative errors). This Minia assembly procedure stops whenever a contig cannot be unambiguously extended.

A graph of contigs is constructed for each insertion site (L, R), as follows. First, an initial contig c L is constructed by calling the Minia assembly procedure from the L k-mer. Given a contig c (initially c = c L ), the four putative neighbors of the last k-mer of c are examined. If no neighbor is present, indicating that c could not be extended, then no further action is performed for this contig. Otherwise, if two or more neighbors are present in the data structure, new contigs will be constructed starting from each of these neighbors. Directed edges will be inserted from c to these new contigs. This process goes on to construct the contig graph in breadthfirst order until a maximum number of contigs (parameter n, usually set to 100) is reached, or a maximal depth (parameter i, usually set to 10 kbp and computed by counting nucleotides in contigs) is reached.

An exhaustive search is performed to find occurrences of R within all contigs in graph, as an exact match (default behavior) or up to a constant number of mismatches. All possible paths between L and R are exhaustively enumerated (i.e. putative insertions). If all paths spell pair-wise identical sequences (minimum identity of 80%), then one of them is returned. Otherwise, the insertion site is considered to be unsuccessfully assembled and all paths are returned. The fill module is performed bi-directionally, i.e. should the (L, R) insertion site yield no path, then the module attemps to assemble the (rc(R), rc(L)) insertion, where rc() denotes the reversecomplement operation.

Evaluation protocol

Simulated datasets

To evaluate MINDTHEGAP, we generated artificial read datasets and reference genomes based on real genomes. First, we simulated sequencing reads for a real genome (the donor). Then, another genome (the reference) was obtained by simulating non-overlapping deletions from a copy of the donor genome. Deletion locations were sampled uniformly along the sequence. These deletions correspond to homozygous insertions in the donor. To simulate heterozygous insertions, reads were sampled in equal numbers from the donor and the reference genomes. Sequencing was simulated with Wgsim from the Samtools package [START_REF] Li | The sequence alignment/map format and samtools[END_REF] using the following parameters: paired-end mode with 2 × 100 bp reads, an insert size of 300 bp (std=50) and a base error rate of 0.01. Coverage was set to 40x for homozygous datasets and 60x for heterozygous datasets.

Three different genomes were used: Escherichia coli K12 (4.6 Mbp), Caenorhabditis elegans (100.3 Mbp) and the human chromosome 22 (35 Mbp without N bases). For each of them, simulated datasets were generated with homozygous or heterozygous deletions of varying sizes.

Assessment of results

Positions of found breakpoints are compared to positions of introduced deletions in the genome. A breakpoint is considered as true positive (TP) if its location is at most 10 bp from a generated deletion position. This margin is meant to take into account fuzzy sites, for which breakpoints are not necessarily found at the exact position of the corresponding deletions (see Methods). For each TP breakpoint, a global alignment between the assembled inserted sequence and the real sequence of the deletion is then performed with needle from the EMBOSS tool suite. We consider the filled sequence as TP if the alignment shows more than 90% of identity. Finally, the recall is the number of TP filled sequences over the number of simulated insertions, and the precision is the number of TP filled sequences over the number of filled insertions.

Real sequencing data

Paired-end sequencing data from C. elegans strain N2 were downloaded from SRA (accession SRX026594). This dataset is composed of 33.8 M Illumina 2x100 bp read pairs (insert size of 350 bp), representing roughly 70x of coverage on the 100.3 Mb reference sequence of C. elegans (downloaded from NCBI version WBcel235). Since we did not find any validated dataset of known large insertions for this genome, we simulated insertion variants following the protocol of simulated data: 1000 regions of a given size (here 1-100 bp or 1 kb) were deleted in the reference genome, corresponding to homozygous insertion variants in the N2 donor genome.

Sequencing data of human individual NA12878 from DePristo et al [START_REF] Depristo | A framework for variation discovery and genotyping using next-generation dna sequencing data[END_REF], consisting of 2.8 G Illumina 2x101 bp read pairs, was downloaded from EBI 1 . The human genome reference assembly (NCBI36hg18) was downloaded from UCSC Genome Browser. A set of predicted or validated large insertions were obtained from supplementary material of Kim et al [START_REF] Kim | Reprever: resolving lowcopy duplicated sequences using template driven assembly[END_REF]. This set contained 30 validated insertions from the study of Kidd et al [START_REF] Kidd | A human genome structural variation sequencing resource reveals insights into mutational mechanisms[END_REF] and 44 which were predicted by Kim et al [START_REF] Kim | Reprever: resolving lowcopy duplicated sequences using template driven assembly[END_REF] based on the alignments of 40 kb sequenced fosmids from this individual to the hg18 reference genome. These are long insertions (median size of 5 kb), 68 are longer than 1 kb, 4 are longer than 10kb. Among these, 61 are predicted as novel insertions, 7 as duplicated and the remaining ones have an unknown status.

RESULTS

Results on simulated datasets

MINDTHEGAP was applied on several simulated datasets to precisely estimate its recall and precision. This enabled to quantify the impact of different levels of genome complexity, to 1 ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/ technical/working/20101201_cg_NA12878/NA12878. hiseq.wgs.bwa.raw.bam independently evaluate each module and modes (detection versus assembly, homozygous versus heterozygous) and to analyze the range of insertion sizes MINDTHEGAP is able to detect and assemble.

High recall and precision in homozygous mode For insertions of 1 kb, MINDTHEGAP recovered between 65% and 98.4% of the simulated insertions, depending mainly on the complexity of the studied genome (Table 1). Almost all predicted homozygous insertions are true positives, resulting in high precision (consistently above 97%). Table 1 shows that almost all insertion sites were detected by the find module in homozygous mode. However, 19-35% of detected insertions could not be assembled by the fill module.

Varying insertion lengths

Figure 3 shows that MINDTHEGAP can detect and assemble insertions of any size. We observed that the performance of the find module is independent of the size of the insertions: recall of the find module never fell below 98.5% (data not shown), without any false positive, even for the human chromosome 22 dataset. However, lower recalls are due to the fill module failing to assemble longer insertions. For small insertions (<100 bp), MINDTHEGAP obtained high recall and precision for all simulated datasets.

Only 650 over 1000 insertions of 1 kb could be assembled in the chromosome 22, and among these, 646 showed more than 90% identity with the original deleted sequences. This was likely due to the high repeat content of this chromosome. We observed that the insertions MINDTHEGAP fails to assemble generally correspond to complex graph of contigs, containing many exact repeats longer than (k -1).

Heterozygous mode

To evaluate the heterozygous mode of MINDTHEGAP, we simulated datasets with only heterozygous insertions (see Section 2.4). Our analysis in Methods showed that heterozygous insertion sites were likely to be more difficult to detect and distinguish from genomic repetitions than heterozygous insertions sites. Table 2 shows that for the human and C. elegans simulated datasets, both recall and precision are significantly below those in homozygous mode. Further investigation showed that the low recall is due to poor performance of the find module. We found that the results in this module were sensitive to the values of parameters r (maximal repeat size at fuzzy sites) and h (maximal number of occurrences of flanking kmers in the reference genome).

k to a higher value and r and h smaller values (here: k = 51, r = 2, h = 1) enabled to reach a precision around 97%, at the cost of a noticeably lower recall. However, using a high k is detrimental to the fill module, due to read coverage being halved in heterozygous insertions. Table 2 shows that on these datasets, the fill module assembles significantly more insertions with k = 31.

Comparison with SOAPindel

On insertions of size 1-100 bp, SOAPindel shows similar recall and precision to MINDTHEGAP (Figure 3). However, SOAPindel is limited in the size of detectable insertions, depending on the insert size of the reads: given our simulation parameters, we observed that SOAPindel recall decreased for insertions larger than 175 bp, and the largest insertion detected was of length 189 bp. Noticeably, the performance of SOAPindel was independent of the genotype of insertions.

Table 1. Precision and recall results for MINDTHEGAP in homozygous mode on simulated and real datasets. Simulated insertions of size 1000 (homozygous). The number of deletions simulated in the reference genome appears in the column 'N sim.'.

Evaluation on a real sequencing dataset of C. elegans

To evaluate the impact of real reads and a real donor genome with some degree of polymorphism in the reference genome, MINDTHEGAP was run on a C. elegans strain N2 read dataset against the reference genome containing simulated deletions. This is to simulate homozygous insertion variants in the donor genome. Note that additional insertions variants are likely to exist C. elegans strain. Thus, the number of FP could be evaluated, as the true set of insertions present in these reads is unknown. For 1 kb insertion variants, 81.1% were correctly predicted and assembled by MINDTHEGAP (Table 1). Compared to the simulated dataset on the same simulated insertions, the find module missed more insertion sites, whereas the fill module had a better recall of inserted sequences. The first observation could be explained by small polymorphism near the insertion breakpoints that generated longer gaps (see Methods), whereas the second by a higher read coverage in this dataset.

Additionally, we compared MINDTHEGAP and SOAPindel on this dataset with 1-100 bp simulated insertions. Recall values were similar for both tools, respectively of 89% and 91%.

Application on real insertions of human individual NA12878

To evaluate the ability of MINDTHEGAP to recover real insertions in real data, we executed it on a human individual NA12878 dataset containing 2.8 G 100 bp reads. As the coverage was high, parameter k was empirically set to 63 and t to 5 (kmers with less than 5 occurrences were discarded). Predictions were then compared to a set of 74 large insertions predicted by alignment of fosmid sequences to the reference hg18 genome (see Section 2.4). 20 insertion sites were recovered by the find module. No heterozygous insertions were predicted. We set r = 15, which enabled to find twice more sites than with r = 5. This suggests that real insertions contain longer repeated sequences at their breakpoints than expected in a random simulation. By analyzing paired-end reads that mapped near each fosmid-predicted breakpoint, we could infer the genotypes: only 23 breakpoints could be confidently assigned to a homozygous genotype (i.e. with less than 5 read pairs spanning the breakpoint). The find module recovered 11 of them. Of the remaining 12 likely homozygous sites, the breakpoints of 8 of them were included in a large gap (≥ k) in the reference binary string. This suggests that these sites were close to other form of polymorphism, which would explain why MINDTHEGAP did not detect them.

Among the 20 detected insertions by the find module, the fill module succeeded in reconstructing correctly 2 inserted sequences of sizes 4137 bp and 6729 bp, with respectively 99.9% and 99.8% identity to the fosmid sequences. This corresponds to a recall of 18%, when comparing to the 11 true homozygous insertions that were detected by the find module. This recall value is similar to one obtained on simulated insertions of 5 kb with the same read dataset (22 %, data not shown).

Time and memory performance

Figure 4 shows the total runtime and maximal memory used by MINDTHEGAP and SOAPindel on the real C. elegans dataset. The machine used for all tests is a 12-core Intel E5-2640 @ 2.50GHz with 192 GB of memory. For MINDTHEGAP, the breakdown of the 3 steps index, find, and fill shows that the major part of running time is spent on the index step. For SOAPindel, the time required to map the reads to the reference with bwa is included in the total time, however SOAPindel alone remains slower than MINDTHEGAP as a whole. SOAPindel used 8 threads and MINDTHEGAP only one.

Importantly, even though MINDTHEGAP stores in memory the whole de Bruijn graph of the C. elegans read dataset, its memory peak (0.7 GB) is 6 times lower than SOAPindel. On the NA12878 dataset with 2.8 billion reads, MINDTHEGAP also proved to scale efficiently: the index / find / fill steps respectively took 32 / 6 / 7 hours, with peak memory usage of 6 / 14 / 6 GB.

DISCUSSION

MINDTHEGAP is the first integrated method to detect and assemble insertion variants of any size and any type, using modest computing resources. The find module of MINDTHEGAP differs from most other existing methods by not relying on read mapping. Instead, the de Bruijn graph of reads is compared against the reference sequence, which enables fast and low-memory analysis. However, one current limitation of the find module that it fails to detect insertions when other polymorphism occurs near the insertion site. Improvements to waive this limitation are under development, based on a more detailed analysis of gaps longer than k. Furthermore, the method could also be used to output SNPs and other types of structural variants.

Long insertion variants are challenging to detect and assemble; thus there is a shortage of tools to compare MINDTHEGAP with. We compared our results with SOAPindel, which is a popular indel detection software limited to short insertions. The NovelSeq software [START_REF] Hajirasouliha | Detection and characterization of novel sequence insertions using paired-end next-generation sequencing[END_REF] is designed to find and assemble large insertions, and therefore would have been another candidate for comparison. However, despite several attempts and reaching out to the author, we were unable to run the software successfully on any of our datasets (the novelseq cluster step ran indefinitely). NovelSeq relies on a complex pipeline, and we conjecture that it may be tailored to specific data types. While most other insertion detection methods require to run external software, MINDTHEGAP is stand-alone and is therefore easy to use. If needed, the modular organization of MINDTHEGAP allows users to replace the find module with the results of a classical insertion detection based on paired-end mapping. The fill module could also be used as a de novo assembly finishing step, i.e. gap-filling between adjacent contigs in scaffolds, although we did not evaluate its performance for this task.

One important design choice for the fill module is to perform assembly with all the k-mers in the read dataset. This enables to assemble not only novel insertions, but also duplicated insertions and transposition events. Classification of assembled insertions into the different event types is not done by MINDTHEGAP, but can be done by re-mapping insertions to the reference genome. One drawback of considering all reads during insertion assembly is that the de Bruijn graph becomes more complex to analyze. An important future work will be to improve the recall of the fill module by using paired-end reads information to guide traversal of contig graphs. As repeated regions are notoriously difficult to assemble, we anticipate that our approach might not be effective for mobile element insertions. However, there exists methods tailored to the assembly of MEI, based on local assembly with recruitment of mate reads.

Our tests on the NA12878 dataset showed there is room for improvement: only 2 long homozygous insertions were successfully assembled out of 23 predicted ones. We postulate that (i) polymorphism or repetitions near the insertion sites hinder detection by the find module, and (ii) the complexity of the human genome makes de novo assembly of large contigs difficult. As no other tool was able to assemble long insertions, we could not assess whether our results were due to weaknesses in our method, or to specificities of this particular dataset (complex insertion sequences or mis-predicted insertions).

ACKNOWLEDGEMENT

  While heterozygous insertions sites do not yield gaps, flanking k-mers at these sites still exhibit features that can GTATTACTATGCTATCTATTATTTA GTATTACTATG.insertion..

Fig. 1 .

 1 Fig. 1. A. Canonical insertion site. The site is detected by its specific signature: the (k-1) k-mers spanning the insertion site are missing in the sequence with the insertion (here k=6), these k-mers are represented as red segments. B. fuzzy insertion site. Insertion ends with the same nucleotides (TG) present on the left of the site. In dashed lines, an alternative insertion site. C. Heterozygous insertion site. Flanking k-mers (in black) surrounding an heterozygous site respectively have two right branching kmers (for the k-mer on the left of the site) and two left-branching k-mers (for the right k-mer).

Fig. 2 .

 2 Fig. 2. Fill module. A graph of contig is constructed from the left flanking kmer L, in a breadth first search order. Construction stops when a maximum number of nodes is reached, or when a branch becomes too deep.The right flanking kmer R is searched within all nodes, finally all paths (in blue) between L and R are outputted as putative insertions.

FindFig. 3 .

 3 Fig. 3. Results of MINDTHEGAP and SOAPindel for several insertion sizes and several genome complexities. SOAPindel results are shown only for insertions of 1-100 bp (first two shaded bars in of each genome section), since it could only detect insertions smaller than 189 bp. Best results of SOAPindel were obtained with k parameter set to 31 (shown here) rather than 51. MINDTHEGAP best results were obtained with k set to 31 for E. coli datasets and 51 for C. elegans and human chromosome 22.

Fig. 4 .

 4 Fig. 4. Time (real time in minutes reported by the unix command time) and peak of memory used by MINDTHEGAP (with k = 51) and SOAPindel (parameter k = 51) on the C. Elegans real sequencing dataset (SRX026594). Peak of memory of SOAPindel approach was reached by the SOAPindel software itself (not bwa).

Table 2 .

 2 Precision and recall results for MINDTHEGAP in heterozygous mode on simulated datasets, containing each 1000 simulated heterozygous insertions of size 1000 bp. Parameter r was set to 2, and assembled insertions smaller than 5 bp were filtered out.

					Find module Fill module k = 51 Fill module k = 31
		Recall (%) Precision (%) N sim. TP	FP	TP	FP	TP	FP
	C. elegans dataset		59.9	93.4	1000 807	11	310	80	599	42
	Human chromosome 22 dataset	35.5	89.0	1000 816	28	226	8	355	44
				4.2GB
		Total 562
		SOAPindel : 232
	Total 69	BWA : 330	0.7GB
	Find : 8 Fill : 9
	Index :52
	MindTheGap	SOAPindel	MindTheGap SOAPindel
	Time (minutes)	RAM (GB)

Funding: This work was supported by the ANR (French National Research Agency), ANR-12-BS02-0008 Colib'read project, ANR-12-EMMA-0019-01 GATB project and ANR-11-BSV7-005-01 SPECIAPHID.

The authors are grateful to Raluca Uricaru for her help and advice.


Anaïs Gouin

Fabrice Legeai

Pierre Nouhaud

Annabel Whibley

Jean-Christophe Simon
email: jean-christophe.simon@rennes.inra.fr

Claire Lemaitre
email: claire.lemaitre@inria.fr

Whole genome re-sequencing of non-model organisms: lessons from unmapped reads

Keywords:

Unmapped reads are often discarded from the analysis of whole genome resequencing, but new biological information and insights can be uncovered through their analysis. In this paper, we investigate unmapped reads from the re-sequencing data of thirty-three pea aphid genomes from individuals specialized on different host plants. The unmapped reads for each individual were retrieved following mapping to the Acyrthosiphon pisum reference genome and its mitochondrial and symbiont genomes. These sets of unmapped reads were then cross-compared, revealing that a significant number of these unmapped sequences were conserved across individuals.

Interestingly, sequences were most commonly shared between individuals adapted to the same host plant, suggesting that these sequences may contribute to the divergence between host plant specialized biotypes.

Analysis of the contigs obtained from assembling the unmapped reads pooled by biotype allowed us to recover some divergent genomic regions previously excluded from analysis and to discover putative novel sequences of A. pisum and its symbionts. In conclusion, this study emphasizes the interest of the unmapped component of re-sequencing datasets, and the potential loss of important information. We here propose strategies to aid the capture and interpretation of this information.

Introduction

Next Generation Sequencing (NGS) and whole genome re-sequencing is nowadays commonly used to identify genomic variants that underlie phenotypic variations, genetic diseases, adaptation or speciation in natural populations. Typically, the reads are mapped against a reference genome and the genotypes (i.e. SNP and structural variant calls) are based on these mapped reads [START_REF] Altshuler | A map of human genome variation from population-scale sequencing[END_REF][START_REF] Nielsen | Genotype and SNP calling from next-generation sequencing data[END_REF]. In addition to universal caveats regarding unknown insertions and/or genomic contamination, which can be overlooked in pure mapping approaches, nonmodel organisms may suffer from the poor quality of the nuclear reference genome and incomplete symbiont or organellar genomes. Moreover, mapping is constrained by the level of divergence between the reads and the available reference sequence [START_REF] Sousa | Understanding the origin of species with genomescale data: modelling gene flow[END_REF]. The resulting ascertainment bias could be problematic, especially when studying adaptation or speciation processes, as genomic regions of interest are expected to display important levels of divergence. These different issues produce a non-negligible fraction of unmapped reads, whose sequences are generally disregarded in favor of the mapped reads in the subsequent steps of the analysis, despite potentially containing useful information. This study offers one strategy for mining the unmapped reads in order to extract relevant biological knowledge, leading to advice and recommendations for other re-sequencing projects.

We investigated this question in the context of a large scale resequencing project on the pea aphid species complex. The pea aphid Acyrthosiphon pisum is a phytophagous insect which feeds on host plants in more than 20 Fabaceae genera. This species forms a complex of sympatric populations, or biotypes, each specialized on one or a few legume species [START_REF] Simon | Host-based divergence in populations of the pea aphid: insights from nuclear markers and the prevalence of facultative symbionts[END_REF][START_REF] Via | Specialized host plant performance of pea aphid clones is not altered by experience[END_REF]. Peccoud et al (2009a) showed that these biotypes include at least eight partially reproductively isolated host races and three cryptic species, forming a gradient of specialization and differentiation potentially through ecological speciation. This complex of biotypes started to diverge between 8,000 and 16,000 years ago, with a burst of diversification at an estimated 3,600-9,500 years (Peccoud et al. 2009b).

In addition, the pea aphid is associated with an obligatory endosymbiont, Buchnera aphidicola, which is found in specialized cells called bacteriocytes and provides its host with essential amino acids. The pea aphid also harbors several facultative symbionts whose distribution is strongly correlated with plant specialization of their hosts [START_REF] Ferrari | Population Genetic Structure and Secondary Symbionts in Host-Associated Populations of the Pea Aphid Complex[END_REF][START_REF] Henry | Horizontally Transmitted Symbionts and Host Colonization of Ecological Niches[END_REF][START_REF] Simon | Host-based divergence in populations of the pea aphid: insights from nuclear markers and the prevalence of facultative symbionts[END_REF] and it has been posited that some of these symbionts could play a role in plant adaptation, although clear evidence is still lacking [START_REF] Mclean | Effects of bacterial secondary symbionts on host plant use in pea aphids[END_REF][START_REF] Tsuchida | Host plant specialization governed by facultative symbiont[END_REF].

This study was carried out on thirty-three aphid re-sequenced genomes from 11 different plant-adapted biotypes. The reads were mapped against the A. pisum reference genome, its mitochondrial genome and its known obligate (B. aphidicola) and facultative symbiont genomes. The A. pisum genome (530 Mb) was assembled using a combination of sequencing technologies (International Aphid Genomics Consortium, 2010;www.aphidbase.com).

Although a second version of the A. pisum reference genome has since been released (International Aphid Genomics Consortium, 2010), the genome assembly remains highly fragmented (23,924 scaffolds) and it has not been subjected to the same level of scrutiny and finishing as the genomes of model organisms such as Drosophila. Moreover, symbiont genome sequences may not be well characterized for this species, and genomic divergence is expected to be important within the whole complex. As a result, a sizeable portion of the reads was not mapped.

In this paper, we scrutinized these unmapped reads by performing crosscomparisons between the sets, assembling the reads by biotype and analyzing the resulting contigs. We used tools developed for NGS such as ABySS [START_REF] Simpson | ABySS: A parallel assembler for short read sequence data[END_REF] and Compareads [START_REF] Maillet | Compareads: comparing huge metagenomic experiments[END_REF], and more classical ones such as the BLAST suite of tools [START_REF] Altschul | Basic local alignment search tool[END_REF]. This analysis revealed that meaningful biological information is contained in the unmapped reads and could help to recover some divergent genomic regions previously excluded from analyses and to discover putative novel sequences of A. pisum and its symbionts.

Material and Methods

NGS data

Thirty-three pea aphid genomes were paired-end re-sequenced using the Illumina HiSeq 2000 instrument with around 15x coverage for each genome. The individuals belonged to different populations each referred to as a biotype due to their adaptation to a specific host plant. In this study, eleven biotypes were each represented by three individuals (Table S1 in Supplementary Material). Reads were 100bp long, sequenced in pairs with a mean insert size of 250 bp and between 32.5 and 59.2 million read pairs (42.5 million on average) were obtained for each individual (see Supplementary Material).

Reads were mapped using Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] with default parameters (up to 10 mismatches per read, or fewer if indels are present -command-line in Supplementary Material) to a set of reference genomes. We also tested another popular mapper, BWA [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler Transform[END_REF], but the percentage of unmapped reads was higher than for Bowtie2 (on average over the 33 individuals, 6.1% vs. 3.7% for BWA and Bowtie2, respectively). The reference set comprised the published pea aphid A. pisum reference genome (IAGC, 2010), and its mitochondrial genome along with the genome of its primary bacterial symbiont and several secondary symbiont genomes reported for the pea aphid (Hamiltonella defensa, PAXS or X-type, Regiella insecticola, Rickettsia sp., Rickettsiella sp., Serratia symbiotica, Spiroplasma sp., Wolbachia sp., [START_REF] Oliver | Facultative Symbionts in Aphids and the Horizontal Transfer of Ecologically Important Traits[END_REF][START_REF] Russell | Uncovering symbiont-driven genetic diversity across North American pea aphids[END_REF]. When available we took the reference genome sequence of the symbiont associated with the pea aphid (i.e0 Hamiltonella defensa 5AT

Extraction of unmapped reads

Fragments for which both reads of the pair did not map to the reference genomes were extracted from the BAM file (mapping result file) using Samtools features [START_REF] Handsaker | Discovery and genotyping of genome structural polymorphism by sequencing on a population scale[END_REF]. In order to check the quality of the unmapped reads, Prinseq [START_REF] Schmieder | Quality control and preprocessing of metagenomic datasets[END_REF] was used.

Sequences were trimmed if, working from the 3' end of the read, base quality dropped below 20 within a window of 10 nucleotides. Read pair information was not preserved and only sequences of at least 66 nucleotides in length were retained for the analysis. Quality-trimmed single-end unmapped read sets were used as the input to the pipeline.

Pipeline for the analysis of unmapped reads

The analysis pipeline, shown in Figure 1, was composed of three major stages: i) pairwise comparisons between unmapped read sets, ii) de novo assembly of pooled sets of reads and iii) analysis of the assembled contigs.

Pairwise comparisons between the read sets were computed in order to identify biologically relevant signals and to define groups of individuals based on the quantity of similar reads. The second part consisted of the assembly of the common reads within previously defined groups. The contigs larger than 1 kb were then analyzed in terms of size, read coverage and similarity to reference genomes.

Figure 1:

Comparison of unmapped reads

Compareads [START_REF] Maillet | Compareads: comparing huge metagenomic experiments[END_REF] was used to compare the read content of the trimmed unmapped read set in each individual in a pairwise manner: this software can find similar reads between two sets of reads in an assemblyfree manner. To be considered a match, a read of set A needs to share at least 2 non-overlapping kmers of size 33 with at least one read of set B.

This comparison thus gives two percentages of similarity between sets A and B: the percentage of reads of A similar to reads of B and vice versa. For all pairwise comparisons, a symmetric similarity score was also provided, computed as follows: The 33 samples were classified based on this similarity measure, using R (version 2.15) software with the maximum distance for the distance matrix and the complete linkage method for hierarchical clustering (function heatmap.2 from gplots package).

AinterB+ BinterA N A +N B ,

Assembly

We pooled common unmapped reads from the three individuals that belonged to the same biotype, i.e. reads present in at least one pairwise comparison between individuals of a biotype were all concatenated in one fastq file. The de novo assembler ABySS [START_REF] Simpson | ABySS: A parallel assembler for short read sequence data[END_REF] was used to assemble these common unmapped reads for each biotype. The size of the k-mer for the De Bruijn graph was set to 31.

To calculate the contig coverage statistics, the sets of unmapped reads were re-mapped to the obtained contig sequences using Bowtie2 (default parameters) and the number of mapped reads was obtained using Samtools.

In accordance with the mean coverage observed in the main dataset (i.e.

where the pea aphid nuclear genome had an average coverage of 15x in each individual), and since reads from two to three individuals were pooled at this step, we considered that contigs with coverage ranging from 20x to 60x likely issued from the pea aphid (nuclear-like coverage), whereas contigs with higher coverage were considered more likely to derive from the symbionts (symbiont-like coverage) or repetitive sequences.

Comparison and analyses of contigs

BLASTClust was used to assess whether large homologous contigs (longer that 1 kb) could be found in different biotypes. A match was retained between two sequences if they were 80% identical over at least 90% of each sequence length. The contigs were then assayed by BLASTn search against the pea aphid reference genomes (nuclear, mitochondrion and symbionts) in order to ascertain their origin. Contigs with hits with an e-value less than 1e-50 were considered to represent highly divergent region of the A. pisum or its known symbiont genomes, i.e. assumed to contain reads that could not be mapped during the first mapping step.

De novo assembly and characterization of an aphid symbiont genome

We performed the following analyses to assemble the genome of a bacterial symbiont detected in the unmapped read set of the individual Vc3.

First, starting from the full read set of Vc3 (40.3 million read pairs), reads were filtered according to their k-mer coverage to obtain only the reads originating from the targeted genome and thus avoid simultaneously assembling the whole nuclear pea aphid genome. Given that the targeted genome had an average read depth in Vc3 of around 600x, only reads for which 68% of the length was covered by 31-mers present at least 100 times in the dataset were retained, using readFilter (P. Peterlongo et al.

unpublished) a custom software based on k-mer counts performed by the DSK software [START_REF] Rizk | DSK: k-mer counting with very low memory usage[END_REF]. Reads that could be mapped to the B.

aphidicola or mitochondrial genomes were removed since their coverage levels would otherwise lead them to be retained in this read set. Only read pairs that remained intact following these filtering steps were included, and these pairs (which totaled 8.8 million read pairs) were assembled using

SPAdes [START_REF] Bankevich | SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing[END_REF], which has been reported to perform well with bacterial genomes [START_REF] Magoc | GAGE-B: an evaluation of genome assemblers for bacterial organisms[END_REF]. Several k-mer sizes were combined in SPAdes (31,41,63,81,89), with default values employed for the other parameters. We kept contigs longer than 500 bp and removed those aligning with the non-Spiroplasma reference genomes. Alignments were performed with the global aligner Mummer [START_REF] Kurtz | Versatile and open software for comparing large genomes[END_REF]. We used GeneMarkS+ [START_REF] Besemer | GeneMarkS: a self-training method for prediction of gene starts in microbial genomes. Implications for finding sequence motifs in regulatory regions[END_REF] to predict proteins in the remaining contigs. These proteins were then compared to the NR database (version 22/01/2014) using BLASTp.

Identification and analysis of potentially divergent regions of the reference genome

To delineate potentially divergent regions of the reference genome which were present in the most divergent biotype (L. pratensis), contigs obtained from the unmapped reads of this biotype were aligned against the nuclear pea aphid genome with Mummer. The regions matching the reference genome with more than 80% identity and larger than 500bp were retained for further analyses.

In these regions, several metrics were computed, including the read depth at multiple mapping stringencies and SNP calling statistics. Read depth was computed first from the initial mapping obtained with Bowtie2 and also following mapping with Stampy [START_REF] Lunter | Stampy: A statistical algorithm for sensitive and fast mapping of Illumina sequence reads[END_REF], an aligner which is reported to perform well when mapping to a divergent reference.

SNP calling statistics were collated from the results of the GATK (DePristo et al, 2011) pipeline applied to the complete dataset of 33 genomes. This pipeline consisted of PCR duplicate removal, indel realignment, base quality recalibration, and genotyping with the UnifiedGenotyper. We used the number of "undefined" calls, i.e. polymorphic positions in the genome 265 for which the genotype could not be determined by UnifiedGenotyper, as a proxy for alignment success. Finally, the gene content of these regions has been established using the version 2.1 of the official gene set of the pea aphid provided by AphidBase [START_REF] Legeai | AphidBase: a centralized bioinformatic resource for annotation of the pea aphid genome[END_REF].

Results

Mapping to reference genomes confirms variation in symbiotic composition between individual host genomes

The coverage of the A. pisum nuclear genome was 14.3x on average (min=10.6x and max=19.96x) while its mitochondrial genome was covered 946.0x on average (min=257.09 and max=3245.60x) and its obligate symbiont genome, 748.8x on average (min=138.08x and max=1509.03x).

The coverage of the facultative symbiont genomes depended strongly on the individual host and varied from 0x to 117.7x. Observed variation for symbiont genome coverage among pea aphids was strongly linked to the infection status of the hosts. Indeed, when we compared expected symbiotic composition based on PCR detection tests and results of mapping, we obtained a good match in most cases: the presence of a given symbiont as detected by diagnostic PCR was confirmed by >2x coverage of reads that mapped against the reference genome (Table S2 S2). Note also that no reads mapped to the Wolbachia genome, confirming the absence of this symbiont in our selection of A. pisum genotypes.

A non-negligible fraction of reads does not map

For a given individual, there were between 0.6 and 7 million pairs of reads (mean = 1.3 million) where both reads did not map to any of the reference genomes (nuclear genome, mitochondrion or known symbionts). This constituted an average of 3.7 % of the initial read sets. Moreover, most of these reads were of good quality, as shown in Figure 2, since few reads were removed (about 17 %) by quality trimming (see Methods).

Figure 2 :

A direct analysis of these read datasets did not allow to characterize the unmapped reads in comparison to the mapped ones, in terms of sequence complexity (Shannon entropy) or signal for repeats (no enrichment of sequence matches with small RNAs targeting pea aphid transposable elements). However, the small size of the reads makes such direct analyses difficult and limits the sensitivity of such characterization.

We can also see in Figure 2 that the fraction of unmapped reads varied between individuals. In particular, the individual Vc3 showed an atypically large amount of unmapped reads with more than 14 million reads representing 18.5% of the initial read set for this individual. For some biotypes, the fraction of unmapped reads was very similar across all individuals, perhaps implying a common cause of mapping failure.

However the fraction of reads did not seem to be correlated with the divergence of the individuals (or biotypes) with respect to the reference genome. The absence of this correlation suggests that the failure to map is not a simple consequence of inappropriate mapping parameters, as if mapping were too stringent we would expect to obtain a correlation between the unmapped fraction and biotype divergence from the reference.

Unmapped reads contain biologically meaningful information

Each set of unmapped reads was compared to all other sets using Compareads. Across the 1056 (33x32) pairwise comparisons, the percentage of common reads between 2 individuals varied greatly, from 6 to 95% with an average value of 50% (see Figure S3 of Supplementary Material). For all but one individual, there was at least one other individual with which it shared 50% of its reads.

Interestingly there was a significant difference when comparing individuals of the same biotype, where on average 70% of reads were shared between individuals, versus comparisons between individuals of different biotypes, which on average shared 48% of reads (p-value <10 -16

for the Welch 2-sample test). This trend was confirmed by the hierarchical classification of individuals based on the pairwise similarity scores computed from the read set intersections (see Methods). Indeed, we can see in Figure 3 that individuals belonging to the same biotype were largely clustered together.

Figure 3:

One extreme case is the L. pratensis biotype, which is known to be the most divergent biotype and is considered a cryptic species (Peccoud et al, 2009a). It showed a very specific profile on the heatmap with strong similarity within this biotype (yellow group on Fig. 3): a L. pratensis individual shared on average 72 % of its unmapped reads with another L.

pratensis individual, whereas only 23 % were shared with an individual of another biotype.

These results show that the sets of unmapped reads contain sequence information specific to biotype or group of individuals, and therefore may contain valuable sequences for biological analyses.

Where do these sequences come from?

In order to get longer and more readily interpretable sequences, we assembled them conjointly by biotype, using the assembler ABySS. Pools of unmapped read sets were used as inputs to obtain sufficient coverage for good quality assemblies. Since the individual classification accorded with the biotype composition and because individuals from the same biotype are genetically closer than those from other biotypes (Peccoud et al, 2009a), we pooled unmapped reads that were shared between at least 2 of the 3 individuals that belonged to the same biotype. By removing reads uniquely present in a single individual, putatively low coverage sequences were excluded, limiting one potential source of noise to the assembly process.

Overall, 94 Mb of contig sequences, each ranging from 100 bp (shorter contigs were filtered) to 35.6 kb, were assembled. On average, 45 % of the unmapped reads could be remapped to the assembled contigs. The average N50 was low (around 428 bp), but we obtained more than 11 800 contigs larger than 1 kb (see Table 1).

The subsequent analysis considered contigs larger than 1kb in more detail. Coverage of the contigs varied considerably, with 57 % of them having a nuclear-like coverage, i.e. between 20 and 60x (see Material and Methods), consistent with an origin from the pea aphid nuclear genome. On the other hand, 14 % of contigs had coverage greater than 60x which would be consistent with an origin from bacterial symbionts, the mitochondrion, or repeated sequences. Contigs with coverage lower than 20x (29%) could correspond to sequences from other microbes (including unreported symbionts) that are in low abundance in the aphid host.

Table 1 Alignment of the contigs to the set of reference sequences can also suggest a potential genomic source. Overall, 63 % of the contigs had a significant blast hit to one of the reference genomes, with the large majority matching with the nuclear pea aphid genome (89 %).

As was found in the coverage analysis, the BLAST analysis revealed biotype-specific trends (Fig. 4). Both approaches can thus be applied to classify the contigs into one of the two main origins: either symbiotic or nuclear. Moreover, the attributions by coverage and BLAST are largely consistent, with a concordant origin for 93 % of the contigs with an origin assigned by both methods. sativum and M. lupulina biotype contig sequences predominantly showed significant similarity to reference symbiont genomes (see Figure 4). In line with this symbiont status these contigs had a high coverage (more than 140x on average). The detected similarity was due to one particular symbiont genome: Rickettsia sp. endosymbiont of Ixodes scapularis (and included in the reference genome set). However, very few reads had mapped initially to this genome, which had an overall coverage of only 3x for the P. sativum and M. lupulina biotypes, and was absent from all other biotypes. This suggested that the chosen reference genome for Rickettsia was too distant from the actual pea aphid symbiont. By comparing these contigs to other Rickettsia species we identified R. bellii as a more closely related species.

This closer relationship was also confirmed by a phylogenetic analysis of 16S ribosomal RNA genes from all available Rickettsia species having their complete genome sequenced (Fig. S1 in supplementary material).

Substituting this genome as a reference resulted in an improved coverage for both P. sativum and M. lupulina biotypes and confirmed the presence of this facultative symbiont in individuals that showed negligible coverage when their reads were mapped to R. ixodes (Table S2 in supplementary material). However, the discrepancy between this coverage level and that observed in the assembled contigs, which was over 2-fold higher, suggests that Rickettsia from A. pisum may diverge significantly from R. bellii which requires to characterize its genome.

Unlike the P. sativum and M. lupulina biotypes, the V. cracca biotype contigs showed almost no similarity to the reference symbionts despite a coverage signal which averaged 602x and was thus consistent with symbiont origin. When aligning these contigs to the NR nucleic acid database, we found few matches and typically only low similarity scores but noted that these hits were enriched in sequences from the Mollicute group and, more precisely, the Spiroplasma genus. This implied that some individual genomes of the pea aphid contained a Spiroplasma symbiont whose genome is distant from available Spiroplasma genomes in the public databases. This hypothesis was confirmed by the fact that three genotypes (of which two V. cracca) of the pea aphid were positive for Spiroplasma infection based on PCR-specific detection (Table S2 in supplementary material). Therefore, Spiroplasma was likely present in at least three individuals, but in high abundance in only one: Vc3. Indeed, most of the V.

cracca unmapped reads came from this single V. cracca individual, which had 5 times more unmapped reads than the average (more than 14 million reads, see Figure 2). Based on the hierarchical classification in Figure 3, Vc3 was grouped with the individuals Vc1 and Tp3 (in agreement with PCR results), with 90.5 and 95% of its unmapped reads being similar to these two individuals respectively. The high abundance of reads from this uncharacterized source in Vc3 led us to attempt the de novo assembly of this Spiroplasma genome. The assembly was performed with SPAdes, after having extracted only putative "Spiroplasma" reads from the full Vc3 read set (see Material and Methods). The final assembly contained 509 contigs longer than 500 bp (2442 bp on average), totaling 1.2 Mb of sequence.

While at the nucleotide level, these contigs showed weak similarity to available Spiroplasma genomes, at the protein level their relationship with Mollicute (and mainly Spiroplasma) proteins was confirmed for 546 annotated genes (on contigs summing to 780 kb). Moreover, the assembly size, low GC content (24%) and the fragmented assembly were consistent with known Spiroplasma genome features: the genome size varies from 1.4 to 1.9 Mb, GC content is around 26 % and the genomes contain lots of repeated sequences and viral elements that make the assembly task harder [START_REF] Carle | Partial Chromosome Sequence of Spiroplasma citri Reveals Extensive Viral Invasion and Important Gene Decay[END_REF][START_REF] Lo | Comparative genome analysis of Spiroplasma melliferum IPMB4A, a honeybee-associated bacterium[END_REF]. Additionally, a phylogenetic analysis of its 16S RNA gene confirmed its membership to the Spiroplasma genus and the absence of any close relative with a complete genome available in the databases (Fig. S2 in supplementary material).

Finally, when unmapped reads were re-mapped to the partially assembled genome of Spiroplasma isolated from A. pisum, individuals which had been found to be positive for Spiroplasma by a PCR-based diagnostic assay but which had a negligible coverage when their reads were mapped to S. melliferum, registered a high coverage on contigs from the A. pisum-derived Spiroplasma (up to 1185x for some individuals, Table S2, supplementary material). In one case, the sensitivity of the sequence analysis may have exceeded that of the PCR test since individual Lc3 was PCR-negative for Spiroplasma but recorded on re-mapping, suggesting a possible infection under the threshold of PCR detection.

Sequences of nuclear origin

All biotypes possessed contigs with a putative nuclear origin, as shown on Figure 4. Some of these contigs were similar between several biotypes or even between all biotypes. We clustered the contigs together using BlastClust and obtained overall 10.1 Mb of distinct sequences having a nuclear-like coverage, of which 4.2 Mb had no similarity to the reference genome of A. pisum. Some of these are likely to be insertion polymorphisms, whereas the 8.6 kb that are shared in at least 8 biotypes could represent pea aphid sequences missing from the current reference assembly either due to error or to deletions in the individual genome that was used to build the reference genome.

Aside from these common sequences, the L. pratensis biotype was particularly enriched in sequences with a putative nuclear origin (Figure 4).

Most of its contig sequences had a significant blast hit to the nuclear reference genome (2.4Mb (69.8 %) of total contig length) and a nuclear-like coverage (86 % of total length), suggesting that these contigs were assembled from reads that were too divergent to map in the first place.

1137 regions (covering 1001 kb) that exhibit similarity to a L. pratensis contig over at least 500bp were then delimited on the reference genome, using the global aligner Mummer. The analysis of read coverage in these regions uncovered two types of region: "low-coverage" regions in which very few reads had mapped (coverage less than 30x for the three L.

pratensis individuals combined, 377 regions summing to 337 kb), and "normal-to-high-coverage" regions (760 regions, 663 kb). While the latter could be explained by one or several divergent copies not present in the reference genome, the former are likely to be regions that are too divergent in all the L. pratensis genomes, in which we may miss important biological information. Indeed, we observed a high proportion of undefined SNP calls for L. pratensis samples in these "low-coverage" regions. On average each L. pratensis individual had 61% of undefined calls, whereas this percentage never exceeded 20% in these same regions for other biotypes (with an average of 13%). These are high values compared to the proportion of undefined calls over the whole genome (on average, 9% for L. pratensis samples and 3.7% for other biotypes). Moreover, half of the regions showed more than 50% of undefined calls for all 3 L. pratensis individuals. This supports the assertion that SNP information is lost because of unmapped reads.

When using a more sensitive mapping approach with Stampy, some of these missed SNPs could be recovered. For the three L. pratensis genomes, overall, 64% of the initially unmapped reads were re-mapped onto the set of reference genomes. Among these rescued reads, 0.66% mapped to the "lowcoverage" regions which was more than expected knowing that these regions of interest represent only 0.06% of the whole genome. This sensitive mapping enabled recovery of, on average, 60% of undefined SNP calls, with 12.5% of regions completely resolved (i.e. with no undefined SNPs). However, for 54% of the regions, the total coverage (Bowtie2 + Stampy) still did not reach normal levels and remained below 30x.

Discussion and conclusion

Whilst approaches for mining unmapped read sets for specific purposes have been described, for example for pathogen discovery (e.g. [START_REF] Kostic | PathSeq: Software to identify or discover microbes by deep sequencing of human tissue[END_REF], this portion of reads is typically disregarded in re-sequencing projects. The sources of unmapped reads are various: they may derive from characterized or uncharacterized symbionts, bacterial, viral or eukaryotic pathogens, highly divergent genomic regions, genomic insertion sequences or library contaminants. The relative proportions of these contributions can vary, and factors such as reference genome quality and the genetic distance between reference and target can play a major role. Therefore, in non-model systems, both the contribution of unmapped reads to the dataset and the likelihood that these reads are a reservoir of useful biological information are increased. However, identifying the unmapped reads and ascribing them to specific sources is not trivial. We have here proposed a novel approach to rescue some of this potentially lost information and have explored the unmapped read sets in the context of thirty-three re-sequenced genomes from biotypes of the pea-aphid species complex.

The direct pairwise comparisons of read sets, before assembly, enabled the rapid identification of similar read sets and highlighted atypical samples and biotypes. Moreover, since the coverage of each individual alone was too low to expect a good quality assembly, merging samples in order to achieve sufficient coverage was necessary for de novo assembly quality. However, selecting and merging only reads common to a single biotype or population, would preclude the identification of other interesting sequences specific to one genotype or to a combination of individuals of different biotypes.

Therefore a more in-depth analysis of the pairwise comparisons followed by the assembly of particular combinations of read sets could be interesting to conduct and may help to uncover unexpected links between individuals.

The assembly phase generated longer sequences than the pre-processed read sets, and these can be more efficiently analyzed and compared to sequence databases. However, although bacterial sequences, such as the ones obtained from Rickettsia and Spiroplasma, could be relatively easily assembled and led to large contigs, we observed that the remaining contigs were usually very short (N50 around 400 bp) and, probably one of the consequence is that a large fraction of the unmapped read sets could not be remapped on the shortest contigs. The recovery of short contigs may be influenced by our methods: extracting and assembling only unmapped read pairs would mean that we assemble only regions of high divergence which may be interspersed in the genome with less divergent regions that are well served by the mapping. In this case, we would predict that samples from the more divergent populations would have, on average, larger contigs of nuclear origin. This is supported in our case by the most distant biotype from the reference, L. pratensis, which shows the greatest proportion of large contigs with similarity with the nuclear genome and de facto the highest number of remapped reads (53.2%).

The final step of our approach was to align the contigs against the reference genomes (nuclear and symbionts) with less stringent similarity criteria than those used during the first mapping step of our process. This, together with the average read coverage of contigs, allowed us to ascribe a putative origin (nuclear or symbiont) of most of the larger contigs. For contigs of symbiont origin, this revealed notably the mis-specification of a reference genome and identified a closer representative species. Without this analysis, we would have concluded from the first mapping that this symbiont was absent (or at very low abundance) from all individuals. Moreover, this revealed the presence in three individuals of a symbiotic bacterium of the genus Spiroplasma, that has been previously reported for the pea aphid [START_REF] Fukatsu | Spiroplasma symbiont of the pea aphid, Acyrthosiphon pisum (Insecta: Homoptera)[END_REF]) but never sequenced and for which we produced a first draft assembly of its genome. Again, the presence of this symbiont would have been completely missed with the first mapping.

In addition, this analysis allowed us to highlight specific parts of the nuclear genome which are enriched in the unmapped read set. These are large regions which are either absent from the reference genome, or show high divergence to the corresponding reference sequence such that each of the read pairs originating from it cannot be mapped. The latter explanation seems to be the most frequent in our dataset. This highlights the major drawback of classical comparative genomics approaches relying on a reference genome. The regions of the reference genome with important genomic divergence for some individuals will contain fewer mapped reads from these individuals and ultimately little divergence will be detected, leading to an erroneous interpretation. We confirmed this consequence of unmapped reads by observing an increased level of unassigned genotypes in these particular parts of the reference genome for the most divergent biotype.

This mapping issue could lead to the loss of valuable biological information or biases in the analysis of genomic variation. Careful calibration of mapping parameters to better handle sequence mismatches between reads and the reference genome can reduce the fraction of reads that cannot be mapped. We explored this by using the Stampy aligner to reprocess the unmapped reads, and could recover 64%. Whilst this offers an improvement on the original Bowtie mapping, most of the observed regions with missing genotype information remained unresolved and it is important to note that relaxing these settings will increase false positive mapping and also increases the time and computing resources required to process the datasets.

Here, our approach helped to recover those divergent regions and, having applied this strategy, the biological signals and functions of these regions can then be interrogated. In the case of the pea aphid dataset, the genic content of the regions will be investigated with a view to determining whether they are enriched in genes involved in host-plant adaptation (e.g. receptors and enzymes). More generally, recovery of these regions enabled them to be subjected to further study, for example to identify signatures of positive selection.

Supplementary information is available at Heredity's website.

Table 1: Contig statistics. For each biotype, the number of unmapped reads in million (n reads) used for the assembly is indicated along with several statistics describing the properties for two contigs length cut-offs (100bp and 1kb), namely the number of obtained contigs (nb), their cumulative length (assbl. Mb), the percentage of reads (% reads) that could be mapped to the contigs and the N50 value.

Legends of figures

  Figure 4:

Figure 1 :

 1 Figure 1: Global overview of the pipeline followed for the analysis of

Figure 2 :

 2 Figure 2: Percentage of unmapped reads (unmapped by pair) for each

Figure 3 :

 3 Figure 3: Hierarchical classification of the sets of unmapped reads. Each

Figure 4 :

 4 Figure 4: Analysis of contigs larger than 1kb in terms of blast matches and

  Figure 1

Acknowledgements

This work was supported by ANR Blanc SPECIAPHID (ANR 11 BSV7 005 01) and by the Région Bretagne (Adaptro). We warmly thank Jean Peccoud for aphid collection and fruitful discussions, and Lucie Mieuzet for DNA extractions.


V Joloboff

S P Wang

Y D Deng

Fast approximately timed simulation

Keywords: virtual prototyping, virtual platform, performance estimation, cycleaccurate simulation, approximately timed simulation, Instruction Set Simulation, ISS, dynamic binary translation, hot path, computer architecture, processor model, System-on-chi, pipelines, software

In this paper we present a technique for fast approximately timed simulation of software within a virtual prototyping framework. Our method performs a static analysis of the program control flow graph to construct annotations of the simulated program, combined with dynamic performance information. The static analysis estimates execution time based on a target architecture model. The delays introduced by instruction fetch and data cache misses are evaluated dynamically. At the end of each block, static and dynamic information are combined with branch target prediction to compute the total execution time of the blocks. As a result, we can provide approximate performance estimates with a high simulation speed that is still usable for software developers.

Introduction

Performance of embedded software is important and must meet the requirements. Many embedded application developers use virtual prototyping framework including a hardware simulator to run the application software over a virtual platform. The simulation community has forged the terms Timeless, Loosely Timed (LT), Approximately Timed (AT) and Cycle Accurate (CA) to characterized simulators, although there is no formal definition of each category. Timeless simulators are purely functional and do not provide any notion of clock, LT simulators provide a notion of time, but the timing provided may be fairly inaccurate, AT simulators are supposed to provide performance estimates within a reasonable error margin and CA simulators (e.g. RTL level) are used to simulate exactly digital circuits. There exists very fast Timeless or LT simulators such as BOCHS [START_REF] Lawton | Bochs: A Portable PC Emulator for Unix/X[END_REF], or QEMU [START_REF] Bellard | QEMU, a Fast and Portable Dynamic Translator[END_REF] or SimSoC [START_REF] Helmstetter | SimSoC: A SystemC/TLM integrated ISS for full system simulation[END_REF] that can run the application software to validate its functionality and possibly test real time software using timers. But such simulators do not provide good enough timings to evaluate the software performance. On the other hand CA simulators are too slow to be used by application software developers. The idea of Approximately Timed simulation is to provide estimates that are as close as possible to the hardware performance, but at a simulation speed that is still an order of magnitude faster. An Instruction-Set Simulators (ISS) is typically used to emulate the behavior of the target on a host machine to execute embedded running on the target processor. In this paper we investigate an approach to provide a fast Approximately Timed ISS that provides good precision estimates. It consists in developing a higher abstraction model of the processor that still execute ISS instructions, but in parallel maintains some architecture state to measure the delays introduced by cache misses and pipe-line stalls, although the pipe-line is not really simulated.

Related work

There has been much work done around Instruction-Set Simulation to reach high simulation speed. Two popular instruments to describe higher abstraction level of hardware models are SystemC [START_REF]SystemC Language Reference Manual (IEEE Std 1666-2005[END_REF] and Transaction Level Modeling [START_REF] Cai | Transaction level modeling: an overview[END_REF][START_REF] Donlin | Transaction level modeling: flows and use models[END_REF], widely used in industry and academia and many efforts have been done to improve Instruction-Set-Simulator (ISS) speed. Whereas the early ISS's used interpretive simulation [START_REF] Sutarwala | Insulin: An instruction set simulation environment[END_REF], most recent ISS's have used some kind of dynamic binary (cached) translation to accelerate simulation [START_REF] Reshadi | Instruction set compiled simulation: a technique for fast and flexible instruction set simulation[END_REF][START_REF] Nohl | A universal technique for fast and flexible instruction-set architecture simulation[END_REF][START_REF] Scott | Retargetable and reconfigurable software dynamic translation[END_REF][START_REF] Shi | An intermediate language level optimization framework for dynamic binary translation[END_REF]. Regarding performance estimate, given the slowness of CA simulation, attention has turned towards sampling. With sampling, only a few chunks of code sequences are selected and analyzed. Then statistical methods are used to generate performance estimates. Popular representatives of sampling methods are Simpoint [START_REF] Hamerly | How to use simpoint to pick simulation points[END_REF], SMARTS [START_REF] Wunderlich | Smarts: accelerating microarchitecture simulation via rigorous statistical sampling[END_REF] and EXPERT [START_REF] Liu | Expert: expedited simulation exploiting program behavior repetition[END_REF]. Sampling techniques can provide fairly accurate performance prediction with a high probability but may also generate very large data files. Another approach consists in statistical workload generation and simulation after collecting data from one complete simulation [START_REF] Nussbaum | Modeling superscalar processors via statistical simulation[END_REF][START_REF] Rao | Hlspower: Hybrid statistical modelling of the superscalar power-performance design space[END_REF]. Because cache misses are the major performance bottlenecks, specific attention has been devoted to cache performance and cache behavior prediction. For example, statistical methods have been successfully used for cache analysis in [START_REF] Berg | A statistical multiprocessor cache model[END_REF]. Yet another approach consists in doing compile time static analysis of the code for predicting cache behavior [START_REF] Vera | Let's study whole-program cache behaviour analytically[END_REF] based on cache miss equations. This method is very fast, but limited in the scope of programs that can be analyzed, and of course requires access to the source program.

More recently, the approach of code annotation has been used, with several variants and our method is affiliated to it as well. The idea consists in annotating the basic blocks of the code back with additional performance or power consumption data that is later used to compute non-functional properties. Such static data is obtained from a model of the target platform in some way. It can be obtained from intermediate compiler representation, or from the real target code of the application [START_REF] Hwang | Cycle-approximate retargetable performance estimation at the transaction level[END_REF][START_REF] Bouchhima | Automatic instrumentation of embedded software for high level hardware/software co-simulation[END_REF][START_REF] Stattelmann | Fast and accurate source-level simulation of software timing considering complex code optimizations[END_REF]. The application code can then be either host compiled [START_REF] Gerstlauer | Abstract system-level models for early performance and power exploration[END_REF] using an abstract RTOS [START_REF] Krause | Combination of instruction set simulation and abstract RTOS model execution for fast and accurate target software evaluation[END_REF], together with the back annotation, or run into a virtual prototype, possibly combining both [START_REF] Gao | Multiprocessor performance estimation using hybrid simulation[END_REF]. Purely static data cannot encompass dynamic behaviour of the hardware, such as cache and branch prediction. Our method consists in annotating the target code with information statically obtained from the binary target application code, based upon an abstract processor architectural model, and then execute that target code with an ISS to obtain the dynamic information. Moreover, we consider an approximation of that dynamic information in order to speed up simulation by introducing usage of hot paths and by constructing tables of bounded size (linear with code size) to store dynamically obtained information that will not need to be recomputed each time. When the table has reached its maximum size, its contents are used statically and are not recomputed.

Approximately timed simulation

Our method is using full system simulation of the target program, with an ISS. During the simulation, cache misses and hits are accurately detected and the control flow graph (CFG) in basic blocks of the program is constructed. A static analysis of the basic blocks is achieved for predicting their performance only once. After that, the simulator collects the performance of the individual basic block. Our work has started from the open source SimSoC Simulator [START_REF]Simsoc open source software[END_REF] that includes SystemC using the TLM standard for communications among simulation models. SimSoC uses a dynamic translation schema whereby the code to be executed is analyzed and the Control Flow Graph of the simulated program is dynamically constructed. After some simulation time, a large part, if not all, of the application basic blocks are stored into the cache. The SimSoC simulator as released in open source is Loosely Timed, using a naïve approach: after executing N instructions, the clock is increased by the time of executing these instructions based on static configuration parameters. Our work has focused on creating an Approximately Timed version of the Power architecture ISS, leveraging off from the existing basic block structure.

Modern processors have complex architectures and can execute a certain number of instructions per clock cycle. There are however several cases where the instructions flow is disrupted, introducing delays in the computation. The major causes of delays are because (i) there are cache misses (data or instruction), or (ii) the pipe-line is stalled or (iii) there are wait states because of communication with peripherals. The latter delays can be captured by TLM transactions. In order to make an AT simulator, our idea is to simulate enough of the processes causing the processor delays, yet we are not simulating the exact hardware. We use instead a model with which the delays can be computed with good approximation while maintaining fast simulation. Our method consists in evaluating these delays with the following approach.

First, we approximate the delays created by the instruction cache misses. We have implemented a cache simulator that does not simulate the specific hardware cache in detail, but reproduces enough of the cache behavior to tell whether there is a cache hit or miss. The instruction buffer is also simulated so that we can compute whether or not the pipe line is fed with instructions. Second, we have built a model of the architecture that makes it possible to evaluate delays without reproducing in detail the hardware. Third, we evaluate with a high precision only the most frequently executed code (the hot blocks), and use a lower precision for the code that is rarely used (the cold blocks), and finally we rely on TLM interface to obtain I/O delays. We perform a static analysis of basic blocks only once to construct the annotations to the code, assuming that future execution of these blocks will take the same number of cycles regarding the pipe line status. Only the data cache hit/misses need to be analyzed dynamically. On average, a program spends a lot of time to execute a small portion of its code, the hot paths. Since accurate performance estimate is costly, an idea to speed up AT simulation is to measure the time spent in infrequently used code (the cold blocks) with a low precision, but fast to obtain, and conversely measure with a high precision the hot blocks. The LT method is used for a block until it exceeds some threshold, then performance estimation switches to the accurate method. This introduces an additional error margin but accelerates simulation and the higher error margin introduced on cold blocks is acceptable. As the threshold is a configurable parameter, the developers can always obtain more accurate simulation by setting the threshold to 1.

The main SimSoC simulator structure is maintained in the AT version. The ISS still executes instructions on the basis of basic blocks and all of the instructions of a block are executed in one step, but at the same time an annotation is constructed. The instructions can be categorized into three categories (a) those that may not stall the pipeline, or those that may, either because they (b) have a dependency on missing data, or (c) influence the branch prediction, which can only be the last instruction of the basic blocks. Therefore we can perform a static analysis of a block to know where the pipe line would stall and compute the delays, to annotate the block with this information. For example, if the pipe-line has K-issues specialized for some instructions, it is necessary to watch for the execution of the basic block once to see which instructions will be parallelized and which instructions will inevitably stall. This simulation does not need to reflect the full pipe line architecture, it is sufficient to understand the dispatching of instructions into the pipe line issues and the stalling factors. After one execution of a basic block, the timing is known and an annotation is constructed.

Given the instruction cache size and cache line size of the target architecture, one can to know whether there is a pre-fetch miss in a block. The simulator checks at the entrance of each basic block whether or not there will be a fetch miss inside that block. As many blocks are small compared to the (instruction) cache line, there may be none. If there is a fetch miss, it is still necessary to simulate the Instruction Buffer process as it has the effect to reduce cache miss delays. Only if the instruction buffer becomes empty while the cache is refilled, then it results into a delay, but in many cases, the instruction buffer process is compensating for the cache miss. This is architecture dependent, as the instruction buffer may be filled at different rates but this can be configured with parameters. Simulating the interconnect behaviour to compute the fetch miss would slow down simulation. However, one can approximate that a delay is constant for all occurrences of a cache miss at a particular position in the basic block. Hence, in long loops over a series of basic blocks, it is not necessary to recompute the delay each time. One can maintain for a basic block a table of fetch miss delays that gets computed during the first occurrences of the loop, and later re-used as constants. This table has a maximum size of the size of the basic block... Note also that when a new cache line of N instructions is entered, then it is not necessary to check for a fetch miss for the following N-1 instructions as it is known to execute sequential instructions that already are in the cache.

A question is the state of the instruction buffer at the entrance of a basic block. There are several possibilities, depending upon the branch prediction algorithm. At the end of the previous basic block, a branch was taken or not. If the branch prediction was wrong, the instruction buffer has been flushed. If the branch prediction was correct, some architecture would have filled correctly the instruction buffer, some others would flush anyway. In this work, we have not considered modelling entirely the branch target prediction and we make the pessimistic assumption that the instruction buffer is empty at the entrance of the basic block. Our performance estimate in this case is worse than reality, which is useful for worse case analysis. In practice, it does not add a considerable error margin as this error only occurs when the branch prediction for not taken was correct and there is a cache miss at the beginning of the block.

Many modern processors, including our case study processor, have a branch prediction unit. These units serve two purposes, to calculate the effective address of the branch in advance, and to fill the cache ahead of time with the appropriate instructions. During simulation, the address of next instruction after the block is known. Thus, one can check, using a model of the processor branch target prediction, whether the predicted target is correct or not, and then possibly compute how many cycles the branch instruction will take.

Regarding the data cache, it must remain fully dynamic as the cache behavior is dynamic, and the application itself may dynamically reset cache options. The SimSoC ISS knows from the target address map whether the address is a real memory access or a mapped I/O. When it is real memory, it uses the Direct Memory Interface (DMI) of the TLM standard to accelerate memory access. We have added a cache algorithm that does not perform actual data caching, but emulates enough of the cache behavior to compute whether there is a cache hit or miss, in cache levels 1 or 2. Again, a pessimistic constant time delay is used for a cache hit or miss at each level, so that the ISS does not need to issue memory transactions and can keep using the DMI interface.

Validation and performance measurements

For this experiment we have considered an embedded processor, namely a Power Architecture EZ model. This processor has a pre-fetch instruction cache, a two issues pipe line, a branch prediction unit, but the pipe line has a small number of stages. The simulation is driven by the SimSoC framework. We have not been able to run the SPEC benchmark programs used by other performance prediction reports because these benchmarks need OS support to construct reporting files not supported by SimSoC Embedded Linux simulator. We can however run open source benchmarks like the libcrypto cryptographic library and various sorting algorithms. Also we have constructed specific small tests to test special cases, that our partner ST Microelectronics could run on their cycle accurate simulator or real hardware, to help us validate our model. As a result, our current system generates in general accurate performance estimation within the error that we had set as our own goal. We have one error due to the fact that, in the current implementation as of this writing, we have not modelled yet the store buffer and this error shows up clearly in the specific benchmarks for store instructions. The impact of running the performance models in addition to the full system simulation has degraded performance, which was expected, but the performance loss is an average of about 33%.

Conclusion

We have presented here a Fast Approximately Timed method of embedded systems simulation. We have added to an existing ISS an abstract model of a processor so that it annotates each basic block with the number of cycles it will take to execute it. This model includes a model of the pre-fetch system, a model of the pipe line and a model of branch prediction. To construct the annotation; it executes the block semantics in one short step, and simultaneously evaluates the time to execute that block. On purpose, for the sake of simulation speed, our model is not totally accurate. It purportedly makes approximations by using constant values instead of issuing transactions to obtain a realistic and it does not maintain micro architecture state across basic blocks. In addition, we have two models, a simplified model and an approximate model. The simplified model is used for cold paths in the code and the more complex, slower, approximate model is used only for hot paths. We obtain an accurate timing only for the most frequently executed code, but this code accounts for a large part of the application timing. Because it uses static analysis of basic blocks to construct the annotation, which is done a limited number of times, it results in reasonably fast simulation, usable by software developers. We have measured our method to be effective to simulate a Power EZ model. Resulting performance estimates are within expectation. This method has the advantage of simplicity; no other external tool is required. The performance estimation is delivered at the end of the simulation session. The same tool can be used by software developers to validate the functionality of application code within the full system simulation of SimSoC framework and obtain performance estimates.

Acknowledgements

We thank ST Microelectronics for their support in this project for comparing simulation results with cycle accurate platforms.