Jia Zhou
email: jia.zhou@cirad.fr

Christophe Proisy
email: christophe.proisy@ird.fr

Xavier Descombes
email: xavier.descombes@inria.fr

Guerric Le Maire
email: guerric.le_maire@cirad.fr

Yann Nouvellon
email: yann.nouvellon@cirad.fr

José-Luiz Stape
email: stape@ncsfnc.cfr.ncsu.edu

Gaëlle Viennois
email: gaelle.viennois@cirad.fr

Josiane Zerubia
email: josiane.zerubia@inria.fr

Pierre Couteron
email: pierre.couteron@ird.fr

Mapping local density of young Eucalyptus plantations by individual tree detection in high spatial resolution satellite images

Keywords: crown identification, object detection, stochastic point process, forests, remote sensing, Brazil

Local tree density may vary in young Eucalyptus plantations under the effects of environmental conditions or inadequate management, and these variations need to be mapped over large areas as they have a significant impact on the final biomass harvested. High spatial resolution optical satellite images have the potential to provide crucial information on tree density at an affordable cost for forest management. Here, we test the capacity of this promising technique to map the local density of young and small Eucalyptus trees in a large plantation in Brazil. We use three Worldview panchromatic images acquired at a 50 cm resolution on different dates corresponding to trees aged 6, 9 and 13 months and define an overall accuracy index to evaluate the quality of the detection results. The best agreement between the local densities obtained by visual detection and by marked point process modeling was found at 9 months, with only small omission and commission errors and a stable 4% underestimation of the number of trees across the density gradient. We validated the capability of the MPP approach to detect trees aged 9 months by making a comparison with local densities recorded on 112 plots of ~590 m² and ranging between 1360 and 1700 trees per hectare. We obtained a good correlation (r²=0.88) with a root mean square error of 31 trees/ha. We generalized detection by computing a consistent map over the whole plantation. Our results showed that local tree density was not uniformly distributed even in a well-controlled intensivelymanaged Eucalyptus plantation and therefore needed to be monitored and mapped. Use of the marked point process approach is then discussed with respect to stand characteristics (canopy closure), acquisition dates and recommendations for algorithm parameterization.

1.

Introduction

Local density in forest plantations corresponds to a planting density from which local mortality is substrated. It is one of the most important factors controlling the productivity of managed Eucalyptus forests (e.g. [START_REF] Walters | Saligna eucalyptus growth in a 15-year-old spacing study in Hawaii[END_REF]. It directly affects inter-tree competition for light, water and soil nutrients for the whole rotation period [START_REF] Gonçalves | Silvicultural effects on the productivity and wood quality of eucalypt plantations[END_REF]. As young trees are the most sensitive to adverse environmental conditions and inappropriate practices, mortality is likely to occur at early stages. There is thus a silvicultural challenge specifically to locate areas of excessive tree mortality or insufficient growth as soon as possible during the first year in order to decide on replanting accordingly. Areas of strong mortality or poor growth are generally conserved up to the end of the rotation and result in overall lower productivity [START_REF] Stape | The Brazil Eucalyptus Potential Productivity Project: Influence of water, nutrients and stand uniformity on wood production[END_REF].

Spatially explicit information about individual Eucalyptus trees may therefore be of great value for pinpointing problems that are often encountered even in intensive plantations due to inadequate management, heterogeneous soil conditions (e.g. [START_REF] Garau | Water stress tolerance in Eucalyptus globulus Labill. subsp. maidenii (F. Muell.) saplings induced by water restrictions imposed by weeds[END_REF], competition with herbaceous vegetation, seed genetic capacity, nutrients supply problems or other local damage [START_REF] Pallett | The relative contributions of tree improvement and cultural practice toward productivity gains in Eucalyptus pulpwood stands[END_REF]du Toit et al., 2005). For instance, trees that are rachitic in their early stages are later totally dominated by 'regular' trees, and by the end of the rotation their crown is expected to be under the canopy. Those trees are generally neglected in biomass estimates by companies.

Therefore, monitoring the actual number of regular-sized live trees in order to make comparisons with planned density turns out to be of prime importance in the management of extensive forest plantations. However, counting and positioning such young individual trees in the field is generally an unaffordable business when dealing with hundreds of hectares of plantation, as it is frequently observed in Brazil.

The potential of sub-metric spatial resolution remotely sensed data needs to be examined carefully with a view to count young trees at an affordable cost. Published studies have used airborne Light Detection and Ranging (Lidar) and optical satellite instruments. Both types of data can provide fine-scale information on canopy characteristics (e.g. [START_REF] Richardson | Strengths and limitations of assessing forest density and spatial configuration with aerial LiDAR[END_REF][START_REF] Hirschmugl | Single tree detection in very high resolution remote sensing data[END_REF]. The Lidar instruments have a good potential for single tree detection in open canopies (e.g. [START_REF] Tesfamichael | A semivariogram approach for estimating stems per hectare in Eucalyptus grandis plantations using discrete-return lidar height data[END_REF][START_REF] Vauhkonen | Comparative testing of single-tree detection algorithms under different types of forest[END_REF] or for a forest carbon inventory (e.g. [START_REF] Stephens | Airborne scanning LiDAR in a double sampling forest carbon inventory[END_REF] due to their ability to measure the vertical structure of forests by determining the distance between the sensor and a target within the canopy. However, the cost of airborne Lidar data collection remains high. In Brazil, it ranges from US$ 15/ha to US$ 200/ha depending on the location inside the country, the area covered and the signal density of points (E. Gorgens, pers. comm.). The additional time and cost required for data processing cannot be ignored [START_REF] Hummel | A comparison of accuracy and cost of LiDAR versus stand exam data for landscape management on the Malheur National Forest[END_REF]. Consequently, operational use of Lidar remote sensing over extensive Eucalyptus plantations in Brazil is still limited to one-off and well-funded research activities.

On the other hand, space-borne imagery acquired at 0.5 to 1-m resolutions and provided by recent satellites such as Geoeye, Ikonos, Worldview or Quickbird is now increasingly available at reasonable costs ranging from US$ 0.10 for archived images to a maximum of US$ 0.80/ha in tasking mode (e.g. e-GEOS, 2011). Fast delivery of data, daily revisit and viewing agility are also key characteristics of a forthcoming constellation of small satellites such as Pleiades. In sub-metric optical images, individual tree crowns are visible because of the contrasting sunlit and shadowed parts of the crown itself and also the gaps between neighboring trees. These data can therefore be used to study the canopy aspect of natural tropical forests by conducting a textural analysis and extracting texture indices that correlate with above-ground stand parameters [START_REF] Couteron | Predicting tropical forest stand structure parameters from Fourier transform of very high-resolution remotely sensed canopy figures[END_REF][START_REF] Barbier | The variation of apparent crown size and canopy heterogeneity across lowland Amazonian forests[END_REF][START_REF] Song | Estimating average tree crown size using spatial information from Ikonos and QuickBird images: Across-sensor and across-site comparisons[END_REF] and biomass [START_REF] Proisy | Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images[END_REF].

However, canopy texture approaches usually yielded less accurate estimations of stand density than other stand variables such as mean quadratic diameter (e.g. [START_REF] Ploton | Assessing aboveground tropical forest biomass using Google Earth canopy images[END_REF].

The possibility of automatically detecting individual tree crowns from optical images has been investigated by a variety of techniques, but they all work on pixel intensity to delineate tree crowns from their surrounding environment [START_REF] Korpela | Individual tree measurements by means of digital aerial photogrammetry[END_REF][START_REF] Erikson | Comparison of three individual tree crown detection methods[END_REF]. A first family of techniques includes local maxima detection (e.g. [START_REF] Pitkänen | Individual tree detection in digital aerial images by combining locally adaptive binarization and local maxima methods[END_REF][START_REF] Wulder | Comparison of airborne and satellite high spatial resolution data for the identification of individual trees with local maxima filtering[END_REF][START_REF] Gebreslasie | Individual tree detection based on variable and fixed window size local maxima filtering applied to IKONOS imagery for even-aged Eucalyptus plantation forests[END_REF], valley following or region growing (e.g. [START_REF] Gougeon | A Crown-following approach to the automatic delineation of individual tree crowns in high spatial resolution aerial Images[END_REF][START_REF] Culvenor | TIDA: an algorithm for the delineation of tree crowns in high spatial resolution remotely sensed imagery[END_REF], 3D template matching (e.g. [START_REF] Larsen | Optimizing templates for finding trees in aerial photographs[END_REF] or a combination of several of the above (Pouliot and King, 2005).

Another family of approaches applies stochastic geometry models to forest scenes.

These models can efficiently detect objects by considering both the noise in data and object variability. In remote sensing, marked point process modeling, a subfield of stochastic geometry, would appear to be well suited to the analysis of very high resolution images since it adjusts a random process by considering groups of pixels as parametric, a priori defined objects located in the geographic space [START_REF] Descombes | Effect of site management on leaf area, early biomass development, and stand growth efficiency of a Eucalyptus grandis plantation in South Africa[END_REF].

The different possibilities opened up by this approach include the counting of items as diverse as brain lesions [START_REF] Descombes | An object-based approach for detecting small brain lesions: application to Virchow-Robin spaces[END_REF], trees [START_REF] Perrin | Forest resource assessment using stochastic geometry[END_REF] and flamingos [START_REF] Descamps | An automatic counter for aerial images of aggregations of large birds[END_REF]. Although the review of [START_REF] Larsen | Comparison of six individual tree crown detection algorithms evaluated under varying forest conditions[END_REF] concludes that the above-mentioned individual tree detection techniques show decreasing accuracy with the heterogeneity of forest stands, the marked point process has been found to yield the best results over a forest plantation. The latter study was based on aerial photographs with pixel sizes ranging from 10 to 50 cm. To our knowledge, young stand density has not yet been assessed through individual tree detection from satellite optical imagery.

In this paper, we assess the potential of marked point processes (MPP) to detect individual tree as a contribution to the monitoring of young Eucalyptus plantations from three Worldview images acquired at the ages of 6, 9 and 13 months. We first used the results of visual detection by independent assessors to measure the performance of the MPP approach and its sensitivity to parameter values. We then validated the MPP estimations of local tree density using figures from a ground truth dataset of 112 plots.

Finally, a local density map was produced across the entire study area i.e. 170 ha. The usefulness and performance of MPP for mapping the early stages of Eucalyptus plantations are finally discussed.

Data

Main characteristics of the Eucalyptus plantation

This study took place in a plantation managed by the Duratex company. It is located in southeastern Brazil at approximately 22°58'S and 48°43.75'W (Fig. 1). The area is currently the subject of several research investigations aimed at monitoring and improving the growth of Eucalyptus trees (e.g. [START_REF] Laclau | Biogeochemical cycles of nutrients in tropical Eucalyptus plantations: Main features shown by intensive monitoring in Congo and Brazil[END_REF]le Maire et al., 2011a). Mean annual rainfall in the area is around 1300 mm, occurring mainly between where the land slopes gently down to a stream. This part of the stand has also a higher percentage of clay content and is more prone to encroachment by herbaceous vegetation [START_REF] Campoe | Stand-level patterns of carbon fluxes and partitioning in a Eucalyptus grandis plantation across a gradient of productivity[END_REF].

The study concerned a surface area of about 170 hectares, as shown in Figure 1. This stand is planted with a company-improved, highly productive clone of Eucalyptus grandis (W. Hill ex Maiden). Clone cuttings have been planted in rows approximately 3 meters apart, and the trees in the rows are about 2 m apart. This gives a target density of 1666 trees per hectare. Planting operations took place over a period of 3 to 4 weeks from the end of November 2009. The soil was fertilized and chemical weeding was carried out to remove competitors [START_REF] Garau | Water stress tolerance in Eucalyptus globulus Labill. subsp. maidenii (F. Muell.) saplings induced by water restrictions imposed by weeds[END_REF].

Unloading tracks clearly divide the area into 4 zones labeled A to D (Fig. 1). The four zones have specific plantation row azimuths (Fig. 2) of about 159°, 72°, 7° and 8° from the East. Despite uniform plant material, soil preparation and fertilization (le Maire et al., 2011b;[START_REF] Campoe | Stand-level patterns of carbon fluxes and partitioning in a Eucalyptus grandis plantation across a gradient of productivity[END_REF], each zone shows local productivity gradients related to local differences in topography, soil water, nutrient availability or insufficient weed control. This alters final planting performance by reducing growth or increasing mortality during juvenile stages in the least favorable contexts.

The growth of these trees averages 20-25 Mg of dry matter ha -1 year -1 , i.e. one of the highest forest productivities in the world [START_REF] Marsden | Relating MODIS vegetation index time-series with structure, light absorption and stem production of fast-growing Eucalyptus plantations[END_REF]le Maire et al., 2011a).

Tree crown extension was recorded in four inventory plots of 84 trees each (Fig. 1) in from May to October 2010 by measuring ground projections of both row and inter-row tree crown axes (diameters). Our measurements in May, August and October showed that crown radii were 0.80m ± 0.22m, 0.99m ± 0.21m and 1.1m ± 0.16m in the plantation row direction and 0.82m ± 0.24m, 1.06m ± 0.22m and 1.15m ± 0.16m in the orthogonal direction. Interestingly, the dissimilarity between the two directions proved weak over the growing sequence (less than 6%) and crown projection to ground was therefore approximated as circular discs. After October, no measurements were available but we assumed relative stable crown projections to ground due to the canopy closeness.

Validation data

The validation data included local tree densities collected in May 2010 on 112 plots located at the periphery of the plantation zones (Fig. 1). These plots corresponded to seven replications of clone trials involving sixteen different clones commonly used in commercial Eucalyptus plantations. The 16 clones were very different in terms of leaf properties (reflectance, transmittance) and canopy structural characteristics (leaf angles, leaf area index, crown diameters, height, etc.). On each of these 112 plots, 12 rows with 16 trees per row were planted on an approximately 3*2 m regular grid at the same time as the rest of the plantation. Within the central area of each plot, we measured the crown diameters of a total of 96 trees (12 trees per row). The area over which density was actually assessed in each plot was thus around 590 m 2 . In some plots, mortality occurred between the planting date and the inventory date, and some trees though surviving, had almost no growth. For the study, we considered as 'regular' only the trees that had a crown diameter above 20 cm in May 2010. The local density on the 112 plots of these regular trees ranged between 1360 and 1700 trees per hectare.

Remote sensing images and preprocessing

The study involved three high spatial resolution Worldview images acquired from May 2010 to January 2011 (Table 1). Only panchromatic data were used. The images were acquired at different viewing angles (Table 1) and were provided in a GeoTiff format at a spatial resolution of 50 centimeters. Sun elevation in the January image was 25° higher than for May and August. Sun-viewing angle values indicated an almost rear illumination (sun behind) of 20° for the May acquisition whereas illumination in August and January was from the side, with angles of about 95 and 105° (Table 1). A preliminary visual examination of the images (Fig. 2) showed how the January image contrasted sharply with the two previously acquired images as being slightly blurred in appearance with plantation rows no longer visible, especially in zones A and C. We interpreted this as a consequence of high sun elevation combined with the relative azimuth angle of the sensor while the distance between the crowns of Eucalypt trees in adjacent rows had narrowed since August (and crowns adjoined along a row). Indeed, rows in the zone B image were still clearly visible in the January image because the sunlight was coming from a near orthogonal direction of 80°.

The registration of the Worldview-2 images needed to be corrected in order to reach the pixel size precision required to superimpose the three images of the same area. To do this, we took the May image as the reference image for superimposition and located about 100 image-to-image landmarks spread over the stand and the neighboring region in both the August and January images. We then used a third order polynomial to reach a spatial superimposition error of less than 50 cm.

Lastly, we re-sampled the three superimposed images into 10 cm pixel images by applying a cubic convolution fit to keep the rectification smooth. This oversampling, as routinely accomplishable by all image processing software, was required to improve the flexibility of disc radius parameter adjustments to the crown in the MPP analysis.

Methods

The study comprised several operations (Fig. 3) applied to the preprocessed images that could be grouped into five main steps: tree positioning by visual analysis (3.1); automatic detection of trees by MPP (3.2); identification of the best sets of parameters by an analysis of detection quality (3.3); validation from ground-assessed densities before generalization to the whole plantation (3.4). These operations were implemented in Matlab® language.

Visual tree positioning

In our study we assumed that a careful visual analysis of the Worldview images would provide an adequate reference dataset of tree positions for comparison with MPP results (Fig. 3; step 2). This is reasonable provided that tree crowns are neither too small (compared to pixel size) nor sufficiently large to fill the inter-crown gaps along the rows.

We therefore visually analyzed the August image (9-month-old trees), which was immune to both difficulties and used the May image (6-month-old trees) only for additional comparisons. In addition, we restricted the tedious work of visual interpretation to a sample of training/validation zones given that not only the affordability but also the quality and objectivity of the interpretation process were likely to decrease with the interpreted area. We thus arbitrarily located and delineated a total of 18 square 50x50 m plots in zones A, B, C and D (Fig. 1) with the aim of interpreting a sample that featured all the different plantation row angles.

Visual detection was carried out by two experts working independently on GIS-based data and consisted in determining latitude and longitude projections of apparent crown centers for all identifiable trees located in each of the 18 plots. About 7,000 trees were located by each of the two experts with close agreement between the two for the number of trees in the plots (correlation coefficient 0.95). On the other hand, tree density across the 18 plots varied substantially, but independently of the zones as the number of visually detected trees ranged from 358 to 417, i.e. the maximum number of trees imposed by the plantation scheme for 2500 m².

Tree detection by the marked point process

The point process theory can be used to build random spatial patterns of points and is part of stochastic geometry. As the reader can find abundant literature on the subject from the mathematical point of view (e.g. [START_REF] Stoyan | Stochastic Geometry and Its Applications[END_REF][START_REF] Van Lieshout | Markov point processes and their applications[END_REF][START_REF] Jacobsen | Point process theory and applications: marked point and piecewise deterministic processes[END_REF], or from an image processing point of view [START_REF] Descombes | Effect of site management on leaf area, early biomass development, and stand growth efficiency of a Eucalyptus grandis plantation in South Africa[END_REF], we present here only the basics of marked point process modeling along with specificities for application to high spatial resolution images of Eucalyptus plantations (Fig. 3; step 2).

Marked point process modeling from zenithal images of forests

A high-resolution image of a forest canopy is a two dimensional lattice of spectral reflectance where tree crowns are distinguishable from their surrounding environment by intensity contrasts between sunlit and shadowed canopy portions or soil portions [START_REF] Proisy | Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images[END_REF]. As visual analysis naturally tends to individualize tree crowns, it is sensible and logical to model the image scene as a particular production of a point process having as marks discs that approximate the crown shapes (Fig. 4; left). Each MPP corresponds to a random configuration of discs belonging to the configuration space M P X u , where P is the space containing the geographic positions of the crown centers, and M is the space containing the marks. In our study, the mark was regarded as the ground-projected crown half-diameter (or radius). Space X was thus defined as:

] [ ] ; 0 [ ] ; 0 [ max min; r r Y X M P X M M u u u (1)
where X M is the width and Y M the length of the image scene, and r min is the minimum value and r max the maximum value considered for the radius of the discs. The PPM model was defined through a probability density with respect to the Poisson measure [START_REF] Van Lieshout | Markov point processes and their applications[END_REF] and an energy function was stated using a Gibbs formulation [START_REF] Descombes | Object Extraction Using a Stochastic Birth-and-Death Dynamics in Continuum[END_REF][START_REF] Descombes | Effect of site management on leaf area, early biomass development, and stand growth efficiency of a Eucalyptus grandis plantation in South Africa[END_REF]. Finding the best configuration for the model then became an energy minimization exercise through the space of all possible random configurations featuring a finite number of disc-objects. The total energy U(x) of a given configuration was defined as the sum of two terms, namely:

) ( ) ( ) ( x U x U x U d p (2)
where U p (x) is a term introducing constraints between two adjacent objects (prior energy) and U d (x) is a term that expressed how well the system made by the objects and their neighborhood fitted image reflectance (data energy). The first term U p (x) cumulates the elementary energy computed for each adjacent pair of objects while the second This local data energy U d (x) was computed through a function of the radiometric contrast between the pixels inside a given disc-object and the pixels composing its surrounding annulus. The higher this contrast, the more likely was the existence of a sunlit crown over a shadowed background. Let us note as µ i , ı i , n i and µ 0 , ı 0 , n 0 the mean intensity, standard deviation and number of pixels inside a given disc-object and inside its surrounding ring, respectively. For every disc brighter than its corona i.e. µ o <µ i , a contrast function derived from a Student's mean-comparison test is computed as follows:

0 0 2 0 0 2 0 1 1 1 2 ) ( n n n n n n x S i i i i V V P P (3)
Then, data energy is derived from the comparison of S(x) with a contrast threshold T (T>0 and T defined and used as part of model parameterization:

°°® ! otherwise , 1 ) ( T - 1 S(x) , ) ( ) ( x S if T x S T x U d (4)
when S(x) < 1 -T, U d (x) < 0, the disc fits the data. Setting the contrast threshold to a low value tends to increase the number of discs accepted around the same location, while high values of T may increase omissions in the event of low contrasts.

The prior energy term U p (x) is defined here to penalize any overlap between two objects and is a function of the normalized extent of the discs' intersection (Fig. 4; left). This function is parameterized using an overlapping threshold value, noted T, that must be adjusted in relation to both image characteristics and forest type. Finally, overall prior

energy is described as ) ( ) ( x x U p JK , where ) (x

K

is the number of disc pairs with nonnull intersection in the current configuration and J is a parameter that weights prior energy comparatively to data energy U d (x).

Multiple Births and Deaths optimization algorithm

In order to find the configuration that minimized total energy U(x), we used a 'multiple births and deaths' process. The algorithm, embedding an annealing technique, was first proposed in [START_REF] Descombes | Object Extraction Using a Stochastic Birth-and-Death Dynamics in Continuum[END_REF] and consists of the following steps:

1) Initialize both a temperature E to 1 and a birth rate 3) Sort the discs of the current configuration with respect to their data energy in decreasing order (i.e. from the worst to the best cases). 4) Remove each disc with a probability depending on both the temperature parameter and the total energy difference between the configurations with and without this object. 5) If the process has not converged (i.e. if total energy changes significantly during the last iteration), both temperature and birth rate are decreased.

3.3.

Identifying the best sets of parameters respectively. Various values for J, 7 and T were also used (see 4.1). For each parameters combination, all detected disc-objects (center coordinates plus radius) were drawn and overlaid on the analogue patterns provided by the visual interpretation of each expert (Fig. 3; step 3). A search was conducted for all detected discs located within a capture distance of less than the maximum radius of the tree crown expected at a given observation date. Capture distance was taken from field measurements to be the average tree crown half diameter at the considered age plus one standard deviation of the related crown radius distribution.

When only one disc was captured in the neighborhood of a known tree, the detection was considered as validated (situation 1). Detection process limitations generated three other possible situations: no discs captured (situation 2), multiple captures (situation 3), and false detection i.e. no reference position could be associated with a given detected disc (situation 4). We therefore defined four variables for characterizing the quality of detection for any given set of parameters. Situation 1 (one tree, one disc) was described through the number of validated discs (noted nbv). Situation 2 informs on omission error as quantified by the number of omitted discs (nbo) whereas situations 3 and 4 contribute to the commission error through the number of multiple detections (nbm) and the number of false detections (nbf). The number of trees visually detected was noted nbr.

Associated rates were then computed from ratios of nbv, nbo, nbl or nbv on nbr.

Finally, and as done in several other studies (e.g. [START_REF] Pitkänen | Individual tree detection in digital aerial images by combining locally adaptive binarization and local maxima methods[END_REF]Pouliot et al., 2005) we used an overall accuracy index, noted AI, to measure the performance of any given set of parameters:

100 (%) u nbf nbm nbo nbv nbv AI (5)
This index corrects the validated rate by taking account both of omission and commission errors. It should be noted that computation of this overall index was strictly dependent upon the availability of reference data for tree locations.

Local density validation and large-scale mapping

We validated the potential of the MPP to detect young Eucalyptus trees by comparing local density values estimated by MPP with fully independent data coming from the 112 field plots measuring ~590 m² as previously presented. Only the best set of parameters that minimized unavoidable omission and commission errors was used (Fig. 3; step 4).

Our last aim was to compute a MPP-derived local density map for the entire plantation area from the number of MPP-detected discs (Fig. 3; step 5). In order to cover the whole Eucalyptus plantation imaged by Worldview, the local density index had to be computed over approximately 800 image windows (each of 50x50m), i.e. about 170 hectares of plantation. We also indicated an error index for MPP-detection by averaging the number of omitted, redundant and false discs, noted nboe, nbme and nbfe, respectively.

For information, it takes around 20 minutes to process a 50x50m image (500x500 pixels)

extracted using a Matlab® routine on a 3.5GHz 64 bit PC running with a conventional RAM configuration of a few gigabytes. Processing time increases with the specified range of radii.

4.

Results

The quality of MPP detections was analyzed in two steps. The first aimed to identify both the parameters (J,7,T) and observation dates that maximized the overall accuracy index.

It was carried out using broad ranges of crown radius values. The second step evaluated how the range of crown values, i.e. [r min , r max ] pairs affected the final results.

Selecting tree detection parameters

MPP parameters

The May, August and January images were analyzed using various pairs of minimum and maximum radius values, i.e. [5, 8], [8, 10] and [10, 12] pixels, respectively. Detection scores were systematically computed in reference to the two visual expert appraisals.

For the sake of conciseness, only a subset of the results and figures obtained is presented in this paper.

Detection scores varied only by a few percentage points when data energy threshold value 7 was varied from 0.05 to 0.8 (not shown). This parameter was therefore set to 0.1 in subsequent analyses, without any great loss of generality. The effect of the J parameter that weights prior energy relative to data energy was also tested by considering values ranging from 0.5 to 6. For each value of J, we simultaneously examined the effect of the overlapping threshold by taking 0 T . The results obtained for all observation dates and all T values showed that setting J below 1 generated far too many multiple detections while J values above 1 reduced multiple detections to less than 15% for the May image and to around 10% for both the August and January images (result not shown). Also, overall detection quality was stable for all values of J t 1. Thus, for the remainder of the study, J was set to 2.

Variations in detection scores were analyzed by considering a range of overlapping threshold T values from 0 (no overlap) to very high values of up to 0.5 for May and August (50% overlapping) and 0.7 for January. We observed that the overall accuracy index for the May image (Fig 5; top) was relatively stable at about 60%, with T less than 0.2. Beyond this value, the overall accuracy index rapidly deteriorated due to a simultaneous decrease in validated detections and an increase in multiple detections.

Image acquisition date

As shown in figure 5, the accuracy index AI for MPP-based detections was clearly higher when using the August image than when using the May and January images, with mean values of 82%, 51% and 30%, respectively.

Considering the May acquisition, index values were about 70% as long as the overlapping threshold T remained below 0. crown enlargement (Fig. 2) is a likely explanation for this poor performance of the MPP analysis at this plantation age.

However, it should be noted that false detections (disc-objects that cannot be associated with a visually detected tree) were zero for all tested parameters and all observation dates, indicating that the marked point process had been appropriately parameterized.

Range of disc radii using the August image

Here we focused on analyzing the August image. The main parameters were set to the previously established values of J, 7 and T in order to determine how the choice of minimal/maximal radius values impacted omission and commission errors and overall detection quality (Table 2) for the four zones with different plantation row orientations (Fig. 6).

Since no significant difference was found between the tree density values provided by the two independent visual expert appraisals (see 3.1), densities were considered as a sound basis to compute overall detection scores. Interestingly, the validated detections rate exceeded 82% for all four zones, showing that row orientation had little or no impact on this result. However, the overall accuracy index indicated that median density values were about 15% lower in A and B than those observed for zones C and D. Rates of no detection and multiple detection ranged from 5 to 8% in zones A and B while in zones C

and D these values did not exceed 3%.

A detailed analysis using various [r min , r max ] pairs (Table 2) indicated that the MPP analysis performed least well taking values in the upper part of the realistic tree crown radius interval, i.e. with r min SL[HOV P ZKHUHDV XVLQJ >r min , r max ] pairs close to the lower boundary i.e. [7,11], [8,11] or [9,11] boosted both the overall index and the validated rate to values above 90% by minimizing omission and commissions errors to a few percent. The potential of the MPP analysis is illustrated in Figure 7 which shows examples of the detection results obtained with [r min , r max ] = [9,11].

Additionally, a comparison of visually-vs. MPP-detected trees pointed to MPP detection providing a 4% underestimation, but despite this difference the two assessments appeared to be linearly correlated throughout the range of Eucalyptus tree densities.

This indicates that the discrepancy between the visual-and MPP-based density estimates is not density-dependent when considering Worldview images acquired for 9month old Eucalyptus plantations.

Validation of MPP-derived local density

We computed the local density on the 112 validation plots using the best parameter combinations i.e. 7=0.1, J=2 and T=0.1 and using one of the best radius ranges, i.e. [7,11] pixels. The comparison of the MPP-derived and ground truth collected local density values correlated well (r²=0.88) with a root mean square error of 31 trees/ha (Fig. 8). The slope coefficient of the regression between both datasets did not significantly depart from 1 with a 95% confidence interval. There was no mean to quantify omissions or commission errors, since ground truth plots did not feature XY locations of trees.

However, this result validated the ability of MPP detection to assess tree density using the parameter values previously identified by visual detection. It also confirmed the robustness of MPP density estimates across gradients of tree densities, tree structural characteristics and row orientations.

Large-scale mapping of local density

Lastly, we generated a local tree density map for the entire plantation (~800 image windows of 50x50m) using one of the best radius ranges, i.e. [8,11] pixels. This parameterization yielded a number of trees ranging from 333 to 420 trees to compare with the 415 theoretically planted trees for a 50x50m planted area. The associated values for nboe, nbme, nbfe corresponded to 17, 14 and 0 trees, respectively. The resulting map showed local density variations ranging from 316 trees mainly in the northwestern part (zone A) to maximum density values in zone D (Fig. 9). Local density across the plantation was clearly not uniformly distributed.

Discussion

Potentialities and further improvements of MPP analysis

High spatial resolution satellite images can supply useful quantitative information for forest management as they can be used to analyze the forest scene on an individual tree scale. For instance, pixelwise classification algorithms based on optical image radiometry are generally ill-suited to detect young individual trees and measure their characteristics. The problem of tree crown delineation in fine spatial resolution images calls for methods that combine both an analysis of radiometric contrast and a priori knowledge of tree size and physiognomy. In response to this, and once it had been verified that tree crowns can be represented by circular discs, the marked point process modeling of forest image properties has an undeniable potential for monitoring Eucalyptus forest plantations. One particular advantage of using this MPP method is its ability to integrate a priori and ideally field-derived data on forest characteristics when setting the ranges of parameter values that determine point process configurations.

Values for disc radii must match realistic ranges of tree crown half-diameters on the observation date in order to avoid false detections and improve MPP performance. Data on tree crown overlapping area are crucial for MPP analysis and can be deduced from ground-truth references in one or more plots. Yet, very little crown overlap is observed in this kind of even-aged Eucalyptus plantation, and, over time, overlap is of very limited importance. The overlap parameter, T, is therefore more intended to express possible crown deformation that may occur from tree-tree interactions along the plantation rows.

Crown modeling via circles cannot account for such deformation plasticity, but allowances made for some overlapping provide the necessary plasticity. However, in our study, field data showed that little crown deformation occurred, and interpretation of the optimal overlap values is probably related to sun and/or sensor position distorting apparent tree position or crown shape.

Some difficulties may arise when transforming radiometric variations expressed as square pixels into disc-objects of adequate size. Here, we found in our preliminary tests that simply oversampling the original image up to 10 cm pixels greatly improved detection as the marked point process functions better when disc-objects are described by tens of pixels rather than a few pixels. More accurate disc positioning through to a finer evaluation of the radiometric distance between a crown and its periphery may improve data energy minimization. This of course relies on choosing an appropriate value for threshold parameter 7. In our case of Eucalyptus crown detection in Worldview panchromatic images and for the 10 cm resolution used, a sensitivity analysis led us to the conclusion that substantial variations in 7 do not significantly affect detection.

Dealing with other images such as infrared channels or Lidar-derived canopy-height models would require additional sensitivity testing. Additionally, computation of the radiometric distances needed to evaluate the data energy term may be improved, especially if spatial resolution and pixel numbers are high.

Prospects and limitations

Although the set of parameters used to implement the method is theoretically large, we have shown here that, for a specific class of application (i.e., a certain type of tree stand and a given type of image), most parameters have a limited effect on the result and can be safely assigned to predetermined values. In our study of a Eucalyptus plantation, we noted that we could take 7=0.1 and J=2 as the energy threshold and weighting of data energy, respectively. We also propose that the overlapping parameter may be taken in IXUWKHU VWXGLHV DV T in the absence of more specific information. The only remaining parameters concern the range of disc radii, which has been shown to have a substantial effect on the results, but which can easily be adapted to each specific case by field measurements of crown radii.

This work highlights that the number of trees and possibly tree locations can be mapped from metric spaceborne images during early growth stages before tree crowns join.

Considering Eucalyptus plantations, mapping can be used to illustrate that planting success is heterogeneous in space and pinpoint where the problem of tree mortality or insufficient growth is the most acute. The potential for stand density mapping from individual tree detection must now be tested in diverse types of even-aged, open-canopy stands. For instance, the monitoring of ecosystems marked by frequent primary successions, such as Amazonian mangroves, could benefit from the mapping of natural colonization over extensive mud banks by fast-spreading species such as Avicennia germinans [START_REF] Proisy | Mud bank colonization by opportunistic mangroves: A case study from French Guiana using lidar data[END_REF]. The method could also be used to detect emergent tree crowns sparsely distributed over a number of hectares of closed canopy forest, and this could be tested by setting the overlapping threshold and disc radii at values distinct from those observed for canopy trees. However, as suggested by the analysis of the January image, corresponding to trees aged 13 months, tree canopy closeness makes detection impossible with the available native image resolution. In comparison, [START_REF] Vauhkonen | Combining tree height samples produced by airborne laser scanning and stand management records to estimate plot volume in Eucalyptus plantations[END_REF] concluded that airborne scanning laser data allowed tree detection during the whole rotation time in 4x3m Eucalyptus plantations. Moreover, when considering natural tropical terra firme forests, tree crown size cannot be restricted to a given interval of a few decimeters. In addition, tree crowns often overlap or closely adjust one to another, resulting in non-circular shapes. Hence, any attempt to determine the number of trees per hectare (except emergent/dominant trees) by MPP analysis on images with a spatial resolution of 50 cm to 1m will probably fail in natural forests. A texture analysis providing information on modal crown size on a forest hectare scale [START_REF] Couteron | Predicting tropical forest stand structure parameters from Fourier transform of very high-resolution remotely sensed canopy figures[END_REF][START_REF] Proisy | Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images[END_REF] may be the best way to extract information from images of closed canopies at these resolutions.

Prospects for improving the management of Eucalyptus plantations

Large-scale automatic mapping of local tree density within commercial eucalypt plantations is possible by MPP-based single tree detection. Generated maps contribute to our overall knowledge of plantations, and could influence decision-making by locating site-specific silvicultural failures that reduce productivity. An early detection of high mortality areas could lead to modifications in current fertilization and weeding practices to adapt them more to local low stand density (precision forestry). Local densities are a good indicator of the potential volume at harvest, particularly when combined with information about average individual tree volume obtained from adequate field measurements. The advantage of using imagery compared to forest inventories is that individual-scale and thus local information can be obtained over large areas. A lower than expected tree density can inform the manager of possible problems occurring in that area, and this can either be dealt with directly or at the next rotation. Since the best detection results are obtained when trees are aged 6 to 9 months, this is too late to replace missing trees, but the manager may decide to replant the local area entirely if economically profitable, or change the clone, treatments or other operations at the next rotation. The price of high spatial resolution satellite images will very likely decrease in the future while availability will increase. Moreover, the application of such MPP-based analyses at each successive rotation may provide a basis for building temporal archives of plantation performance that may be useful for adjusting stocking densities and management practices. Finally, the potential of the MPP approach will probably increase with its application to fin-scale Lidar-derived canopy heights models that already proved usefulness for monitoring Eucalyptus plantations (e.g. [START_REF] Packalén | ALS-based estimation of plot volume and site index in a eucalyptus plantation with a nonlinear mixed-effect model that accounts for the clone effect[END_REF].

Concluding remarks

A better understanding of how the most influential MPP parameters need to be set for different forest types, tree ages, crown forms and sun and viewing conditions may be obtained by simulating reflectance images using 3D stand mockups and a radiative transfer model, like for instance the DART model [START_REF] Gastellu-Etchegorry | 3D modeling of satellite spectral images, radiation budget and energy budget of urban landscapes[END_REF]. Here, as explained by [START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF], realistic canopy images can be simulated by applying the DART model to 3D forest mockups consisting of "lollipop" trees (ellipsoids on a stick) built from field inventories and allometric relationships. This is a very convenient method for developing and testing the robustness of remote sensing methods dedicated to the application of high spatial resolution imagery to forests.

Simulating images of Eucalyptus plantations may help to both adjust the set of MPP parameters for each age (and related levels of total height and crown development) and determine optimal acquisition parameters in terms of sun and viewing angles (le [START_REF] Le Maire | Very high resolution satellite images for parameterization of tree-scale forest process-based model[END_REF]. This approach may pave the way to providing key attributes for forest management, such as trees location maps, an abacus of crown diameter growth during early stages, and areas of poor growth or high mortality. Information such as this is crucial when investigating the effect of planting density on yield [START_REF] Vanclay | Modelling Forest Growth and Yield: Applications to Mixed Tropical Forests[END_REF] or when

simulating competition for light [START_REF] Binkley | Explaining growth of individual trees: Light interception and efficiency of light use by Eucalyptus at four sites in Brazil[END_REF]le Maire et al., 2012a;[START_REF] Binkley | Light absorption and use efficiency in forests: Why patterns differ for trees and stands[END_REF] and when using a simple representation of tree crowns with geometrical volumes, such as ellipsoids or cones [START_REF] Medlyn | A MAESTRO Retrospective[END_REF]. The inset shows the position of the illustration within the embedding system of rows. 23 24

  October and March. Annual temperatures range from 13.3 °C in the coldest months (June to August) to 27.2 o C in the warmest season (December to February), with 19°C being the average. In addition, temperatures in the cold season rarely drop below 5 o C. The study area is almost entirely flat except in the western and northwestern part of the stand

  U d (x) cumulates the local energy and includes variations in pixel values within and around each object in the configuration.

  image pixels randomly from those not yet associated with a disc. Then, select discs with radius between r min and r max randomly from the newly born disc radius. a uniform distribution.

  MPP modeling was applied to the 18 image extracts (each of 50x50m). The range of disc radius variations was adjusted based on field measurements of the radii of a sample of crowns at the three plantation ages considered (see 2.1). The [r min , r max ] pairs values used to analyze the January image were taken according to the October measurements, assuming minor changes in canopy dimensions afterwards due to canopy closure. The largest interval range for [r min , r max ] pairs was thus set to[5, 10], [7, 12] and [9, 15] pixels in the analysis of the May, August and January 10 cm oversampled images,

  2, but then fell to 40% (Fig 5; top). The rate of missed detections exceeded 15% while that of validated detections (single detections) was less than 60%. Conversely, the August image clearly provided good overall scores exceeding 80% when using T (Fig 5; middle). Very poor detection quality was obtained from the January image, with the overall accuracy index peaking at no more than 35% for T = 0.4 (Fig 5; bottom). Canopy closure along plantation rows due to

Figure 1 .

 1 FIGURE CAPTIONS 1 2 Figure 1. Extract from the Worldview image acquired in August 2010 of the Eucalyptus 3 plantation area and location of the study site near São Paulo, Brazil (box map). The four 4 zones of different plantation row azimuths are noted A to D. The four filled frames inside 5 zones C and D indicate the areas where extensions of tree crowns were recorded on a 6 total of 4x83 trees. The empty square frames symbolize the 18 training plots (each of 7 50x50m) from which the visual analysis was completed. The filled white square frames 8 represent the 112 ground truth plots where tree counting was obtained from a field 9 inventory undertaken in May 2010. 10 11 Figure 2. Appearance of the four plantation row orientations (from left to right) from 12 selected 50x50m Worldview images in May, August and January (from top to bottom). 13 Crown size variability and variation of local tree density of trees are also apparent. 14 15 Figure 3. Flowchart of processing steps. 16 17 Figure 4. Illustration based on a particular configuration of three disc-objects of the main 18 aspects of marked point process fitting. Left: Illustration of overlapping area ǹ between 19 two disc-objects and principles defining the prior energy term in relation to the 20 overlapping threshold T. Right: Example of circular disc-objects as positioned on a 10 cm 21 pixels image to maximize reflectance contrasts between discs and surrounding rings. 22The inset shows the position of the illustration within the embedding system of rows. 23 24

Figure 5 .Figure 9 .

 59 Figure 5. Effect of overlapping threshold T on detection quality for May (6 months after 25 planting, top), August (9 months, middle) and January (13 months, bottom) images. 26 Circular dots indicate accuracy index AI (see text). Square, star and triangle markers 27 correspond to rates of validated, multiple and omitted detections, respectively. All values 28 are averaged over three plots from those inside zones A, C or D (one per zone) using 29 three values of data energy threshold 7 values of 0.05, 0.1 and 0.2. The standard error 30 in overall accuracy index is represented by dashed lines. 31 32 Figure 6. Results from the MPP-based modeling of the image acquired in August with 33 J=2, 7=0.2 and T=0.2 (see text). Detection rates and accuracy index are computed for 34 the four zones of the plantation and for each expert visual interpretation as median 35 values over the [r min , r max ] intervals listed in table 3. Circular dots indicate the accuracy 36 index AI. Square, star and triangle markers correspond to rates of validated, multiple and 37 omitted detections, respectively. 38 39 Figure 7. Example of MPP detection results obtained over the 4 plots presented in figure 40 2 using the best parameters combination established for the August image analysis i.e. 41 J=2, 7=0.2 and T=0.2 and [r min , r max ]=[8, 11]. Yellow discs with '+' indicate single 42 detections. Blue square dots and green triangles with 'x' indicate multiple detection and 43 no detections, respectively. 44 45 Figure 8. Comparison between local density values obtained from a field survey and 46 from MPP analysis. Two plots of strongly contrasting tree density are shown. The solid 47 black line gives the regression line of the form X a Y u (intercept coerced to zero). 48 49

Table 1

 1 Main acquisition parameters for Worldview images.

	Satellite	Date of acquisition	Acquisition time (GMT)	Spatial resolution (m)	Sensor azimuth angle (°)	Sensor elevation angle (°)	Sun Azimuth angle (°)	Sun Elevation angle (°)
	Worldview2 11-May-10	13:29	0,5	53,5	70,7	33,5	68,1
	Worldview2 01-Aug-10	13:43	0,5	297,7	83,2	32,5	42,7
	Worldview1 29-Jan-11	13:50	0,5	337,6	81,6	82,3	66,7

Table 2

 2 Median scores obtained using various [r min , r max ] pairs and in comparison with the two visual expert appraisals. Minimum and maximum scores are given in brackets.

			Number	AI (%)	Omitted detections Single detections Multiple detections
	r min	r max	of
			couples	Expert #1 Expert #2 Expert #1 Expert #2 Expert #1 Expert #2 Expert #1 Expert #2
				82.8	84.9	1.6	2.2	90.6	91.9	5.4	3.7
	7	8 -12	5	(68.4-	(69.1-	(0.0-	(0.0-	(81.2-	(81.7-	(0.7-	(0.9-
				95.2)	94.8)	10.7)	12.4)	97.5)	97.3)	15.2)	13.3)
				83.2	86.4	3.9	3.8	90.8	92.7	3.3	2.6
	8	9 -12	4	(70.1-93.3)	(72.1-94.8)	(0.7-10.7)	(0.2-13.1)	(82.4-96.5)	(83.8-97.3)	(0.7-11.6)	(0.7-7.9)
				81.2	84.8	5.7	4.4	89.6	91.4	2.9	2.3
	9	10 -12	3	(70.6-95.7)	(70.9-94.7)	(1.6-13.5)	(0.8-14.2)	(82.7-97.8)	(83.0-97.1)	(0.5-9.3) (0.7-7.3)
				75.9	76.4	11.4	11.3	86.3	86.6	1.9	0.8
	10	11 -12	2	(65.2-91.1)	(66.2-92.7)	(3.8-19.0)	(3.5-19.1)	(78.9-95.4)	(79.6-96.2)	(0.5-5.4) (0.0-2.8)

Acknowledgements

This work was supported by the ORFEO program, a Centre National d'Etudes Spatiales (CNES) Accompaniment Program for the use of PLEIADES images. We especially wish to thank Claire Tinel for her help in image sorting and are grateful to the Eucflux project (http://www.ipef.br/eucflux/) for providing inventories of the Eucalyptus plantation plots.

We thank Peter Biggins for the revision of the English.


Evguenia Kopylova
email: evguenia.kopylova@lifl.fr

Laurent No

H Él Ène Touzet

SortMeRNA: Fast and accurate filtering of ribosomal RNAs in metatranscriptomic data

Motivation:

The application of Next-Generation Sequencing (NGS) technologies to RNAs directly extracted from a community of organisms yields a mixture of fragments characterizing both coding and non-coding types of RNAs. The tasks to distinguish among these and to further categorize the families of messenger RNAs and ribosomal RNAs is an important step for examining gene expression patterns of an interactive environment and the phylogenetic classification of the constituting species. Results: We present SortMeRNA, a new software designed to rapidly filter ribosomal RNA fragments from metatranscriptomic data. It is capable of handling large sets of reads and sorting out all fragments matching to the rRNA database with high sensitivity and low running time.

INTRODUCTION

The application of NGS technologies for metatranscriptomic profiling has been a successful venture in practice. Scientists may now gain access to the full set of coding and non-coding RNA in a community of organisms, which becomes particularly important for samples which cannot be cultivated outside their native environment [START_REF] Stewart | Metatranscriptomics analysis of sulfur oxidation genes in the endosymbiont of solemnya velum[END_REF][START_REF] Bomar | Directed culturing of microorganisms using metatranscriptomics[END_REF][START_REF] Shi | Metatranscriptomics reveals unique microbial small rnas in the ocean's water column[END_REF]. The initial challenge of metatranscriptomic sequenced data analysis is to sort apart the RNA fragments based on their biological significance. Messenger RNAs (mRNA) cast a universal glimpse on the gene expression patterns between interactive species. Likewise, the ribosomal RNAs (rRNA) disclose information on the community's structure, evolution and biodiversity, and prevail in classification and phylogenetic analyses. The rRNA can comprise up to 90% of total RNA. Various prior-to-sequencing procedures, such as mRNA amplification kits, can help to enrich the yield of mRNA [START_REF] Gilbert | Gene expression profiling: metatranscriptomics[END_REF]. However, these kits are not fully satisfactory since secondary steps may be required to verify if the resulting material is an accurate representative of the initial samples [START_REF] Nygaard | Limitations of mrna amplification from small-size cell samples[END_REF]. New software have been recently developed to address this issue, they identify and isolate rRNA fragments from a set of sequenced reads. The first set of programs, Meta-RNA 3 [START_REF] Huang | Identification of ribosomal rna genes in metagenomic fragments[END_REF], SSUALIGN [START_REF] Nawrocki | Infernal 1.0: inference of rna alignments[END_REF] and rRNASelector [START_REF] Lee | rrnaselector: a computer program for selecting ribosomal rna encoding sequences from metagenomic and metatranscriptomic shotgun libraries[END_REF] share a common algorithmic approach, to represent an rRNA family database using a probabilistic model. Both Meta-RNA and rRNASelector use prebuilt Hidden Markov Models (HMM) and consequently sort reads against the database with the HMMER3 package [START_REF] Eddy | Profile hidden markov models[END_REF] whereas SSU-ALIGN uses Convariance Models to support secondary structure information. An alternative algorithm outside the domain of probabilistic models is riboPicker [START_REF] Schmieder | Identification and removal of ribosomal rna sequences from metatranscriptomes[END_REF] which uses a modified version of the Burrows-Wheeler Aligner [START_REF] Li | Fast and accurate short read alignment with burrowswheeler transform[END_REF]. Lastly, BLASTN [START_REF] Altschul | Basic local alignment search tool[END_REF] is used in numerous home-made workflows for this problem. With BLASTN however, reads should be compared against all sequences of an rRNA database to achieve a good sensitivity level. In all cases, computational time is still an issue to handle large collections of reads.

In this paper we describe SortMeRNA, an efficient filter requiring only a representative set for an rRNA database and rapidly sorting through millions of reads. The underlying algorithm is analogous to the seeding strategy, focusing on finding many short regions of similarity between an rRNA database and a read. SortMeRNA also takes advantage of redundancy between homolog sequences, as HMMs do, and builds a compressed model of all rRNA sequences. The generated results adhere to the accuracy of the HMM-based programs and are computed in a fraction of the time.

SYSTEM AND METHODS

Algorithm overview

We assume having a collection of unassembled reads and a database of rRNA sequences, and we want to sort out reads that match to the database. The general principle behind our algorithm is to search for many short similarity regions between each read and the rRNA database. We scan each read with a sliding window, and the accepted reads are those which have more than a threshold number of windows present in the database. For a given read and a given window on the read, we authorize one error (substitution, insertion or deletion) between the window and the rRNA database.

To achieve this task in an efficient manner, the rRNA database is stored in a Burst trie coupled with a lookup table that speeds up the access to the Burst trie and takes advantage of conserved regions in the rRNA sequences. For a given read and a given window on Fig. 1. SortMeRNA Algorithm Overview. The set of representative rRNA sequences is preprocessed in the following way: (1a) The lookup table stores all of the s 2 -mers and their number of occurrences which exist in the rRNA database. (1b) The Burst trie is a data structure which stores the rRNA database. (2a) The algorithm takes as input a collection of reads provided by the user, and for each read, a sliding window w of even length s ∈ [14, 20] moves across the read. (2b) For each window w, the prefix w[1.. s 2 ] and suffix w[s.. s 2 + 1] are translated into a decimal value between 0 and 2 s -1.

(3) If the value exists in the lookup table with a high frequency (see Section 1.1 of the supplementary file), the remaining part of the window is searched in the Burst trie. This is done with a cyclic traversal between the universal Levenshtein automaton and the Burst trie that determines whether the subpattern is present in the rRNA database with at most 1 error. For every letter traversed in the Burst trie, a bitvector is passed to the universal Levenshtein automaton to verify if the number of encountered errors remains less or equal to 1. (4) After all windows have been traversed, if the number of accepted windows exceeds a certain threshold (see in Section 2.5) then the read is accepted and classified as rRNA. the read, we find the set of windows present in Burst trie using the universal Levenshtein automaton. This comparison is done by performing a parallel traversal between the Levenshtein automaton and the Burst trie.

0 #1 1 #1 2 #1 3 #1 4 #1 0 #0 1 #0 2 #0 3 #0 4 #0 0 1 0 0 Σ Σ Σ Σ Σ ǫ ǫ
Figure 1 globally illustrates this framework. The length s of the sliding window is a parameter of the algorithm, further discussed in Section 2.5. The acceptance of a read depends also on the ratio of matched windows. Let r be this parameter. This choice will also be discussed in Section 2.5.

The Burst trie to store an rRNA database

The Burst trie [START_REF] Heinz | Burst tries: A fast, efficient data structure for string keys[END_REF] is a fast and versatile data structure which effectively stores a large number of strings such as an rRNA database. Unlike the standard trie, the binary search tree or other variants, which often adopt an equal rate of memory access among the cache or main memory, the Burst trie can exploit the modern cache architecture by addressing memory closest to the CPU. It is capable of reducing the number of trie nodes by 80% while maintaining performance similar to a hash table [START_REF] Askitis | Redesigning the string hash table, burst trie, and bst to exploit cache[END_REF]. Given a sequence vz, the Burst trie can store the prefix v as a link of trie nodes and the suffix z as an array of characters appended to the last trie node. Normally, subtrees become more sparse in the depth of a trie and representing them as reduced 'buckets' of contiguous memory preserves space and boosts cacheefficiency. When the number of sequences sharing a common prefix v reaches a fixed threshold, the appended bucket of suffixes bursts to form a new trie node and smaller sub-buckets. To optimize memory access during subtree traversal, the threshold size of a bucket should be less than the lower level cache. A systematic use of this trie can be observed in the fastest sorting algorithm for large sets of strings, the Burstsort [START_REF] Sinha | Cache-conscious sorting of large sets of strings with dynamic tries[END_REF].

Fig. 2. Let s = 16, the Burst trie below is constructed on the first six 17mers of an rRNA sequence. The 'char f lag' describes whether a pointer is set to a trie node '1', a bucket '2' or neither '0'. Additional information on the origin of the 17-mer directly follows each element, as shown in the dashed bucket.

A 0 rRN Ai C 1 U 2 A 3 G 4 A 5 A 6 U 7 G 8 A 10 G 11 U 12 U 13 U 14 G 15 A 16 U 17 C 18 U 19 G 20 G 21 C 22 U 23 C 24 A 25 . . . A C U A G A A U G A G U U U G A U C U A G A A U G A G U U U G A U C U A G A A U G A G U U U G A U C U 17 -mer # 3 A 0 C 1 G 2 U 3 1 2 2 2 null null null null char f lag trienode* child node void* bucket uagaaugaguuugauc 2 agaaugaguuugaucu 2 aaugaguuugaucugg 2 A 0 C 1 G 2 U 3 2 2 2 0 null null null null null ugaguuugaucuggc uagaaugauuugau aaugaguuugaucug strand hashid
Following a similar method of an array-structured trie as described in [START_REF] Sinha | Cache-efficient string sorting using copying[END_REF], our Burst trie is assembled exactly on the nucleotide alphabet {a, c, g, u}. As illustrated in Figure 2, the trie stores every unique (s + 1)-mer substring in an rRNA database, since we look at windows of length s with at most 1 error between any two words. The information on whether the (s+1)-mer belongs to a forward strand, the reverse-complement or both (strand), and its origin (hashid) follows each entry in a bucket. When the exact location of the (s + 1)-mer needs to be found in an rRNA database, the hashid value serves as an index in a complementary table storing this information. Nearly one-quarter of the 16S rRNA positions are 99-100% conserved [START_REF] Cannone | The comparative rna web (crw) site: An online database of comparative sequence and structure information for ribosomal, intron, and other rnas[END_REF][START_REF] Mears | Modeling a minimal ribosome based on comparative sequence analysis[END_REF] and this moderates the size of the trie since many identical or closely similar substrings are shared between sequences.

We use an additional optimization to improve access into the Burst trie. Since we consider at most one error between the window and the database, we have this simple property: For every two words such that the edit distance between them is bounded by 1, there exists a common substring of length s 2 which is either a prefix or a suffix of the two words. We apply this property to construct a lookup table storing all s 2 -mers existing in the rRNA database. Note that for s in [14,20], transposing the nucleotide alphabet onto a binary equivalent, such that {a, c, g, u} = {00, 01, 10, 11}, we can represent each s 2 -mer in s bits which maps to a unique integer value. Upon completion of the forward and reverse Burst tries, a scan of each trie is performed to record the existence of all s 2 -mers and, if present, associated pointers to the trie node representing the immediate letter following the prefix. The precomputed lookup table quickly determines whether an exact match of the prefix or suffix exists in the Burst tries and furthermore it provides us with direct access to the remaining part of the word in the Burst trie.

The lookup table also allows us to take into account distribution of s 2 -mers in the rRNA database. A multiple sequence alignment of an rRNA database can clearly define areas of high nucleotide conservation and emphasize the evolutionary origins shared between organisms. In a similar manner, the lookup table defines highly conserved areas by keeping only frequent s 2 -mer occurrences in the rRNA database. Before a window is traversed in the Burst trie, its prefix or suffix must exist in the lookup table. This notion enforces that a read matches closely to one region in a database rather than multiple scattered ones leading to a false alignment (see Section 1.1 of the supplementary file).

The universal Levenshtein automaton

Fig. 3. The nondeterministic Levenshtein automaton for p = acgu and k = 1. The s #e notation for each state corresponds to s number of characters read in the pattern p and e number of errors recorded. The initial state is 0 #0 and the three final states are 3 #0 , 4 #0 and 4 #1 . Each non-final state has three outgoing arcs, one for each type of edit operation.

0 #1 1 #1 2 #1 3 #1 4 #1 0 #0 1 #0 2 #0 3 #0 4 #0 a c g u Σ Σ Σ Σ Σ ǫ ǫ ǫ ǫ deletion ǫ substitution insertion Σ match
The classical nondeterministic Levenshtein automaton for a pattern p and a number of errors k recognizes the set of strings which are at most edit distance k to p (see Figure 3). This automaton is not suitable for computation because of the presence of multiple active states and epsilon transitions. This may be overcome by transforming the automaton into an equivalent deterministic form. However, the resulting automaton may be exponential in the length of p and likewise dependent on it. In [START_REF] Schulz | Fast string correction with levenshtein automata[END_REF] and [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF], a universal Levenshtein automaton was characterized based upon insightful observations of the classical one. The term universal conveys its one-time construction and independency of p. The intuition arises from the symmetry of the nondeterministic automaton, which applies the same set of transition rules to every new input character and each new set of active states is a subset of a known bounded superset. A set of bitvectors symbolizing the homology of p and a candidate string serve as input to the automaton. In full generality, the size of the automaton is exponential in a function of k [START_REF] Mitankin | Universal Levenshtein Automata. Building and Properties[END_REF]. In our case, since k = 1, it remains sufficiently small. The set of bitvectors representing the similarity of two strings are precomputed using the following definition.

Definition 2.1. [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF] The characteristic vector χ(w, V ) of a symbol w ∈ in a word V = v1 . . . vn ∈ * is the bitvector of length n where the ith bit is set to 1 iff w = vi.

The technical details of n ≤ 2k + 2 and the prefix of k symbols of '$' appended to the pattern p can be found in the (2004) paper cited above.

Example 2.1. Let k = 1, the input word W = acaga and the pattern p = $acuaga, then χ1(a, $acu) = 0100, χ2(c, acua) = 0100, χ3(a, cuag) = 0010, χ4(g, uaga) = 0010, χ5(a, aga) = 101 are the computed characteristic bitvectors. It follows that {χ1, . . . , χ5} is the characteristic bitvector array carrying the similarity information of x and p .

Beginning from χ1 to χ |s| , the bitvectors are sequentially passed into the universal Levenshtein automaton. Each bitvector leads to a transition between states (in constant time) corresponding to the number of errors encountered thus far. If some χi reaches a failure state, greater than k errors exist between s and p, and the strings are rejected. The automaton only recognizes two strings if the input of the last bitvector χ |s| leads to a final state.

Match of a read with the dynamic bitvector table

At this point, matching a window w of length s on the read against the rRNA database amounts to first checking whether the prefix or the suffix of length s 2 of w is present in the lookup table, then determining if the universal Levenshtein automaton for w recognizes some word in the Burst trie. For the second step, we have to implement a rapid traversal of the Levenshtein automaton which relies upon the precomputation of bitvectors for w. At every depth of the Burst trie, we assume that the symbol q in χi(q, V ) appears as one of {a, c, g, u} with equal probability. Ultimately during traversal, the bitvector of the actual residing nucleotide is chosen. Figure 4 shows the precomputation of bitvectors for p = $acuaga in Example 2.1. If the string x = acaga existed in the trie, then the highlighted set of bitvectors {0100, 0100, 0010, 0010, 101} would form the bitvector array (see Section 1.2 of the supplementary file for a graphic example).

Fig. 4. The precomputed bitvector table for pattern p = $acuaga covering all possibilities of q for k = 1. The first bit in each entry of column i = 0 represents the $ symbol and is always set to '0'. When the window is shifted by one position, the subsequent pattern p changes simply by the removal of the first character in the prefix and the addition of a new character in the suffix. Hence, rather than recomputing the bitvector table for each new window, a series of bitwise operations are taken to modify it, as demonstrated in Figure 5.

Fig. 5. The modification of the bitvector table from pattern p 1 = $acuaga to p 2 = $cuagaa for k = 1. Columns 0-2 of p 2 are equal to columns 1-3 of p 1 , except for column 0, where the most significant bit (MSB) of every bitvector represents the symbol $ and is set to '0'. Column 3 of p 2 equals to column 4 of p 1 with an additional bit appended. The appended bit is set to '1' in the bitvector correspoding to the newly appended character, otherwise it is set to 0. Column 4 of p 2 is equal to column 3 of p 2 , although the MSB is not considered. The same rule applies to column 5 of p 2 , where the two MSBs of the column 3 bitvectors are not considered. 0100 1000 0000 0000 000 00 0000 0001 0010 0100 100 00 0010 0100 1000 0000 000 00

p 2 : $ i c 0 u 1 a 2 g 3 a 4
a 5

Following a pre-order path, the traversal of the Burst trie begins at the root node. Through knowledge of the nucleotide letter and the depth of the node being visited, the coinciding bitvector is accessed in the precomputed bitvector table, indifferent to whether the node is a trie node or a character in the bucket. Subsequently, the bitvector is passed to the universal Levenshtein automaton which decides whether to continue traversal of the current subtree or backtrack to the first branching point with a non-failure Levenshtein state and recommence traversal of a new substree. In this manner, a complete traversal of the Burst trie remains unlikely as backtracking occurs each time the edit distance between the pattern and a traversed branch exceeds k. To further speed up Burst trie traversal for every window, a 'backwards dictionary' approach as described in [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF] was implemented. The original algorithm builds two dictionaries, one for the forward strings and the second for their reverse equivalents. In this manner, the same window can be traversed quickly from both ends.

Parameter setting

The algorithm depends on two parameters: The length s of the sliding window, and the minimal proportion r of accepted windows in a read. To find a robust choice for s and r, we ran the algorithm for several values of s and r on several rRNA databases and for several sets of reads.

We purposely designed four databases with distinctive features: Small 16S and large 23S subunit, varying identity percentage and from distinct phylogeny tree subparts, Set 1: 16S, 80% identity (2262 rRNA) Set 2: 16S, 80% identity, truncated phylogeny tree (2187 rRNA) Set 3: 23S, 95% identity (1969 rRNA) Set 4: 23S, 95% identity, truncated phylogeny tree (1906 rRNA) These databases were constructed by applying the ARB package [START_REF] Ludwig | Arb: A software environment for sequence data[END_REF] and UCLUST [START_REF] Edgar | Search and clustering orders of magnitude faster than blast[END_REF] to sequences from SILVA [START_REF] Pruesse | Silva: a comprehensive online resource for quality checked and aligned ribosomal rna sequence data compatible with arb[END_REF] (see Section 2.1 in the supplementary file). Next, we constructed datasets of simulated rRNA and non-rRNA reads using the software MetaSim [START_REF] Richter | A sequencing simulator for genomics and metagenomics[END_REF]. We used two sequencing error models, Roche 454 and Illumina, because the errors for Roche 454 mainly originate as indels and for Illumina as substitutions. The length of the reads differs as well: ≥ 200nt for Roche 454 and 100nt for Illumina technology. To test the sensitivity on Set 1 and Set 3, we constructed 300,000 Roche 454 reads and 1,000,000 Illumina reads on the entire SILVA database minus the sequences used for the representative rRNA database. To measure the sensitivity for discovering new species with Set 2 and Set 4, the same number of reads was simulated only on the truncated sections of the Bacteria phylogeny tree. To test the selectivity, the non-rRNA reads were simulated using the NCBI bacterial genomes library with rRNAs masked (see Section 2.2 in the supplementary file).

The parameter values were varied as: s ∈ [14, 20] and r ∈ (0, 1). The main conclusion is that s = 18, r = 0.15 for Roche 454 reads and s = 18, r = 0.25 for Illumina reads give best sensitivy/selectivity balance for all rRNA databases. Moreover, varying r within short ranges does not significantly affect the results (see Section 2.3 of the supplementary file). We use these values as default settings in all subsequent analyses.

IMPLEMENTATION

SortMeRNA is implemented in C++ and freely distributed under the GNU general public license (GPL). It can be downloaded from http://bioinfo.lifl.fr/RNA/sortmerna. The software uses OpenMP functions to parallelize filtering of the reads.

The input criteria are a fasta/fastq file of letter space reads produced by Roche 454 or Illumina technologies, and a fasta file of rRNA sequences. There are eight rRNA databases included in the software package covering the small (16S/18S), large (23S/28S) and 5/5.8S ribosomal subunit rRNAs, which were all derived from the SILVA and RFAM databases. Additionally, the user can work with their own RNA databases.

EXPERIMENTAL EVALUATION

The performance of SortMeRNA was measured in terms of sensitivity, selectivity and real-data analysis compared to the software SSU-ALIGN [START_REF] Nawrocki | Infernal 1.0: inference of rna alignments[END_REF]), Meta-RNA [START_REF] Huang | Identification of ribosomal rna genes in metagenomic fragments[END_REF]), rRNASelector [START_REF] Lee | rrnaselector: a computer program for selecting ribosomal rna encoding sequences from metagenomic and metatranscriptomic shotgun libraries[END_REF]), riboPicker [START_REF] Schmieder | Identification and removal of ribosomal rna sequences from metatranscriptomes[END_REF]) and BLASTN [START_REF] Altschul | Basic local alignment search tool[END_REF]). All tests were performed on an Intel(R) Xeon(R) CPU W3520 2.67GHz machine with 8GB RAM, L1 cache size of 32 KB, L2 cache size of 256 KB and L3 cache size of 8192 KB. Since riboPicker and SSU-ALIGN do not provide a direct option for multi-threading, all tests were carried out using one thread.

rRNA databases

We created two new representative databases: 16S rRNA with 85% identity (7,659 sequences) and 23S rRNA with 98% identity (2,811 sequences) (see Section 3.2 of the supplementary file). The 16S rRNA database was used by SortMeRNA, riboPicker, BLASTN and SSU-ALIGN, and the 23S rRNA database by SortMeRNA, riboPicker and BLASTN. SSU-ALIGN was written for aligning small ribosomal subunits and does not provide models for 23S rRNA. riboPicker was also tested with a more comprehensive database made available from their website: All 16S and 23S rRNA sequences taken from SILVA, RDP-II, Greengenes, NCBI archaeal and bacterial genomes, and HMP (3,232,371 16S and 19,602 23S unique sequences). The results for this larger database are indicated by riboPicker * in the subsequent tables. For Meta-RNA and rRNASelector, we used the HMMs provided with the software.

Simulated reads

Sensitivity for 16S rRNA. 300,000 Roche 454 and 1,000,000 Illumina 16S rRNA reads were simulated in the same manner as described in Section 2.5. The performance results can be viewed in Table 1. All software programs except riboPicker and SSU-ALIGN have a sensitivity level higher than 97%, and even higher than 99% for BLASTN and SortMeRNA. The sensitivity for riboPicker is very low (56%) because BWA-SW works well with error rates 2%-3% for 100-200nt reads, and loses sensitivity for new species. As expected, the sensitivity increases with a larger database (indicated riboPicker * ). Considering the computation time, SortMeRNA runs in less than 2 minutes, or 72x faster than the next fastest tool with proportionate sensitivity (Meta-RNA). Note also that BLASTN executes at a very slow speed (several hours), because reads should be compared against all of sequences in the representative database.

Selectivity for 16S rRNA. 1,000,000 Roche 454 and 1,000,000 Illumina non-16S rRNA reads were simulated in the same manner as described in Section 2.5. The performance results can be viewed in Table 2. All programs have a selectivity level higher than 99.98%.

The number of false positives for the HMM-based programs remains comparable to SortMeRNA for both Illumina and Roche 454 reads. The difference in the simulated data results between Meta-RNA and rRNASelector can be attributed to the number of bacteria vs. archaea rRNA sequences used in the construction of the HMMs, as well as additional parameter settings in rRNASelector. riboPicker * and BLASTN show the lowest selectivity. Concerning the running time, the order of the fastest programs is rRNASelector, Meta-RNA and SortMeRNA. Both rRNASelector and Meta-RNA use the HMMER3 package, which applies a pre-filter to quickly reject sequences which would score very low in the HMM. This acceleration heuristic gives these programs a competitive advantage on the artificial dataset for selectivity where all of the sequences are negative.

Results for 23S rRNAs are analogous in terms of accuracy and running time. They can be found in Table A

Real data

The metatranscriptomic datasets SRR106861 of a photosynthetic microbial community and SRR013513 of a tidal salt marsh creek from 454 sequencing were downloaded from the NCBI Sequence Read Archive. The results for 16S and 23S can be viewed in Table 3 and Table 4 respectively, and the overlap of the results between tools in Figure 6 and Figure 7. The results obtained with SortMeRNA are very close to the ones obtained with HMM-based methods. riboPicker finds only a fraction of all potential rRNAs, which confirms its low sensitivity for small databases. The majority of 16S reads found only by riboPicker * (1,298) map to mRNA. For 23S analysis in Table 4 and Figure 7, approximately 99% of the excess reads of Meta-RNA (12,112) and rRNASelector likewise map to 28S, along with 83% of the (624) reads found only by BLASTN and Meta-RNA. The (537) reads found only by BLASTN map to 16S rRNA, ncRNA and mRNA.

DISCUSSION

SortMeRNA has shown to be a rapid and efficient filter which can sort a large set of metatranscriptomic reads with high accuracy comparable to the HMM-based programs. SortMeRNA implements seeds with errors (substitution and indel) and this important characteristic renders the algorithm robust to errors of different types of sequencers while providing the ability to discover new rRNA sequences from unknown species.

The method used by the algorithm is universal and flexible. The database can be constructed on any family of sequences provided by the user. Moreover, the algorithm does not require a multiple sequence alignment file to build the database, as HMM-based programs do, and this is an advantage when sequences are hard to align or only partial sequences are available. Another advantage of SortMeRNA is the small number of parameter settings required by the program (see Section 2 of the supplementary data).

  g u . . . c a g u c g u . . . c a g c sliding window w of length s 2b. for each window

  and Table B under Section 3.3 of the supplementary file.

Table 3 .

 3 Runtime for the SRR106861 metatranscriptome of 105,873 reads against a 16S rRNA database of 7,659 sequences. Fig. 6. Venn diagram for reads classified as 16S rRNA by BLASTN, Meta-RNA, SortMeRNA and riboPicker * in the SRR106861 metatranscriptome.

		rRNA	run time latency memory (%)
	SortM eRN A	27046 34s	1x	4.8
	riboP icker	11389 4m10s	7x	2.3
	riboP icker *	27195 39m3s	69x	30.8
	BLAST N	27061 1h29m 157x	0.6
	M eta-RN A	27111 10m33s 18x	1.8
	rRN ASelector 27085 10m40s 18x	0.8
					Meta-
	BLASTN
	27, 061 reads	17		26	RNA 27, 111 reads
		1	46	0
	6	1, 104	0	1, 298
		38	25, 893	0
	SortMeRNA	4	0		riboPicker *
	27, 046 reads		0		27, 195 reads

Table 4 .

 4 Runtime for the SRR013513 metatranscriptome of 207,368 reads against a 23S rRNA database of 2,811 sequences. Fig. 7. Venn diagram for reads classified as 23S rRNA by BLASTN, Meta-RNA, SortMeRNA and riboPicker * in the SRR013513 metatranscriptome.

		rRNA	run time latency memory (%)
	SortM eRN A	94395	51s	1x	3.8
	riboP icker	71937	10m2s	12x	3.9
	riboP icker *	84152	36m27s 43x	5.5
	BLAST N	94439	3h42m 261x	0.9
	M eta-RN A	106698 1h33m 109x	4.8
	rRN ASelector 107900 1h36m 113x	3
					Meta-
	BLASTN
	94, 439 reads	537	12, 112	RNA 106, 698 reads
		290	624	0
	3	8, 851	0	14
		1, 151	83, 960	38
	SortMeRNA	0	139		riboPicker *
	94, 395 reads		1		84, 152 reads

Table 1 .

 1 SENSITIVITY. 1,000,000 of MetaSim simulated Illumina (100nt) and 300,000 Roche 454 (≥ 200nt) rRNA reads against a representative 16S rRNA database of 7,659 sequences.

				Illumina				Roche 454
		rRNA	run time latency memory (%) sensitivity (%) rRNA	run time latency memory (%) sensitivity (%)
	SortM eRN A	998615 1m39s	1x	8.5	99.861	299979 1m43s	1x	6.3	99.993
	riboP icker	558607 18m45s 11x	6.8	55.860	123024 18m36s 11x	5.6	41.008
	riboP icker *	999941 6h33m	238x	35.3	99.994	299999 9h	314x	34	99.999
	BLAST N	995322 23h52m 868x	3.0	99.532	299978 18h35m 649x	1.4	99.992
	M eta-RN A	983332 2h	72x	33.3	98.333	299980 1h57m	68x	12.9	99.993
	rRN ASelector 974118 1h47m	64x	17.4	97.411	299976 2h	70x	7	99.992
	SSU -ALIGN	971221 6h49m	248x	0.1	97.122	299902 5h50m	204x	0.1	99.967

Table 2 .

 2 SELECTIVITY. 1,000,000 of MetaSim simulated Illumina (100nt) and 1,000,000 Roche 454 (≥ 200nt) non-rRNA reads against a representative 16S rRNA database of 7,659 sequences.

				Illumina				Roche 454
		rRNA run time latency memory (%) selectivity (%) rRNA run time latency memory (%) selectivity (%)
	SortM eRN A	17	2m9s	2x	7.6	99.9983	13	3m42s	1x	10.2	99.9987
	riboP icker	7	10m22s 8x	6.7	99.9993	3	29m45s 9x	16.8	99.9997
	riboP icker *	158	56m37s 42x	35.1	99.9842	53	2h43m	49x	45.2	99.9947
	BLAST N	33	14m22s 11x	0.3	99.9967	33	16m12s 5x	0.3	99.9967
	M eta-RN A	11	1m33s	1x	0.1	99.9989	11	3m41s	1x	0.2	99.9989
	rRN ASelector 10	1m20s	1x	0.1	99.9990	11	3m21s	1x	0.2	99.9989
	SSU -ALIGN	8	3h51m 173x	0.1	99.9992	11	10h30m 188x	0.1	99.9989

ACKNOWLEDGEMENTS

This research was supported by ANR project MAPPI (ANR-2010-COSI-004), LIFL (UMR CNRS 8022 Université Lille 1) and Inria Lille Nord-Europe. Project MAPPI is associated with the Tara Oceans expedition (oceans.taraexpeditions.org), where the principal tasks involve the development of new software for mapping and assembling metagenomic and metatranscriptomic data.