M Dahoumane

D Dauvergne

J Krimmer

H Mathez

C Ray

E Testa

A H Walenta

Y Zoccarato

A Low Noise and High Dynamic Charge Sensitive Amplifier-Shaper associated with Silicon Strip Detector for Compton Camera in hadrontherapy

come

Abstract-A 8 channel Front End Electronics (FEE) circuit has been designed and fabricated in 0.35 µm CMOS process from Austria Micro System to be coupled with the Silicon Strip Detector (SSD) of the Compton Camera for quality control of hadrontherapy. Each channel includes a Charge Sensitive Amplifier (CSA) followed by two parallel CR-RC shapers. Slow and fast shapers, with 1 μs and 15 ns shaping time, are used to measure the energy and to time stamp all events respectively. The two sides of the SSD are read thanks to a configurable system for holes and electrons. The CSA presents an open loop gain of 67 dB and 90 degrees phase margin assuring a high stability. The circuit has been successfully tested. The test results are in good agreement with analytic and simulation calculations. Here, we describe the principles and present measured performances of the prototype. A high linearity over the range of 3E3 to 3E6 electrons is reached with a conversion gain of 3.6 mV/fC. The circuit achieves an ENC (Equivalent Noise Charge) of 412 electrons rms. 75% of the total noise is generated by the small value of the feedback resistor chosen to avoid pile up phenomenon due to the 1E5 hits/s occupancy rate. A cross-talk of 2 % was measured, 99% of which is due to the power supply disturbances. The power supply dissipation is 21 mW/channel for 3.3 V supply voltage. The area of this design is 2871×1881 μm 2 including pads.

I. INTRODUCTION

adrontherapy is a radiotherapy technique which consists in irradiating tumorous cells with protons or ions. During irradiation, part of the incident projectiles undergoes nuclear fragmentation and prompt-gamma rays can be used to control the ion range, because their emission profile is found to be correlated with the range, for both protons and carbon ions [START_REF] Min | Prompt gamma measurements for locating the dose falloff region in the proton therapy[END_REF], [START_REF] Testa | Monitoring the Bragg peak location of 73 MeV⁄u carbon ions by means of prompt γ-ray measurements[END_REF]. Since excited fragments emit γ-rays almost instantaneously (below the ns range). Prompt-γ can be used to monitor online and in real-time the ion range. A Compton Camera detecting system has been proposed [START_REF] Richard | Design Guidelines for a Double Scattering Compton Camera for Prompt-gamma imaging During Ion Beam Therapy: A Monte Carlo Simulation Study[END_REF], [START_REF] Roellinghoff | Design of a Compton camera for 3D prompt-gamma imaging during ion beam therapy[END_REF] and is expected to be a very efficient and fast technique for the online dose control and high resolved 3D image reconstruction. In the proposed system shown in Fig. 1, the scatter detector is composed of a stack of 2 mm thick doublesided Silicon Strip Detectors (SSD) providing the spatial and temporal coordinates (x, y, z, t) and the energy deposited by the recoil electron during the Compton interaction of γ-rays. The scattered photon is further absorbed in a scintillator absorber. The SSD detector will be coupled with a multi-channel integrated readout electronics fulfilling the following particular requirements imposed by the application: dynamic range of 1000: 0,48 fC to 480 fC, two signal polarities readout: P and N sides of the SSD, counting rate: 1E5 hits/s, noise: ENC = 300 electrons rms, shaping Time (slow) = 1 µs for energy measurement, shaping Time (fast) =15 ns for time stamp and detector capacitance C d = 15 pF. For this aim, a 8-channel Front End readout Electronics (FEE) prototype is designed and fabricated using the AMS CMOS 0.35 µm process. This paper presents the design details and test results of such a prototype circuit. II. DESIGN ARCHITECTURE Fig. 2 illustrates a general diagram of the circuit AC-coupled with a double sided SSD detector. Each channel includes a CSA as a first stage followed by two parallel shapers. The first one, with long shaping time (1 µs), measures the input charge and optimizes the signal to noise ratio. The second one, with small shaping time (15 ns), is dedicated for timing measurement. The detector, with a capacitance C d , produces charge pulses that are integrated on the CSA feedback capacitor C f [START_REF] Paul | Prospects for charge sensitive amplifiers in scaled CMOS[END_REF], [START_REF] Hofmann | Development of readout electronics for monolithic integration with diode strip detectors[END_REF]. The circuit has been designed to readout H both P and N polarities by an appropriate feedback resistor R f or a T resistor network. The voltage references of the shapers are selectable according to the feedback network of the CSA. Classical cascode architecture [START_REF] Radeka | Design of a charge sensitive preamplifier on high resistivity silicon[END_REF] has been chosen to implement the CSA stage of the readout chain as shown in Fig. 3. The geometry of each component is optimized for the foreseen performances. The first order transfer function of the CSA in response to a pulse charge stimulus is given by:

A. The Charge Sensitive Amplifier (CSA) description

f d in out C G C j Zf I V 0 0 1 ) ( ω ω ω + - = , (1)
where Z f is the equivalent feedback impedance composed by R f in parallel with C f . and Cd the detector capacitance. The value of C f (250 fC) is fixed by the maximum signal of the dynamic range and the power supply voltage value which is 3.3 V in 0.35 µm CMOS process we used. The value of the feedback resistance is chosen to discharge the feedback capacitor C f by a time constant

f f Cf R C = τ
and so to avoid the pile up phenomenon which could be caused by the high counting rate of the Compton Camera (1E5 hits/s). Thus, the small value R f = 4 MΩ needed has a very harmful effect on the noise as it will be detailed in the section III. In Fig. 3, the NMOS input transistor M1 has been chosen to get a better power supply rejection ratio and higher speed than the PMOS one. The transistor P1 is used as an active load which is mandatory for wide dynamic range such as that of our application (1000). It sets the bias current of the input transistor M1 at 5 mA from 3.3 V power supply voltage. In order to minimize the ballistic deficit, G 0 C f >>C d , where G 0 is the open loop gain of the CSA, C f the feedback capacitance and C d the detector capacitance. To ensure up to 99 % of signal charge transfer, G 0 should be around 67 dB according to the following formula which gives the charge transfer rate α:

f d C G C 0 1 1 + = α . ( 2
)
The rise time of the CSA output for a pulse charge stimulus is determined by the amplifier time constant τ CSA as:

in d f CSA R C C G f Cd = = 0 0 2π τ , (3)
where

0 0 2 G f c π ω = (4)
is the Gain Bandwidth Product. A pole-zero R and C feed forward compensation technique is used to improve stability without affecting the speed of the CSA [START_REF] Sansen | Feedforward Compensation techniques for High Frequency CMOS amplifier[END_REF]. The Bode diagram simulation of the CSA gives an open loop gain of 67 dB and a phase margin of 90°. The loop gain and Cadence STB tool simulation were also performed. The obtained results ensure high stability and very weak loss of the transferred signal charge. The signal collected on the P strip side of the SSD has a negative shape on the CSA output. The authorized output dynamic range maintaining all transistors in the saturation mode is given by the approximate following formula:

V fF fC C Q V f out 92 . 1 250 480 max 0 - = - = - = , ( 5
)
where Q 0 is the injected charge. The total dynamic range corresponding to both polarities equals ±1.92 V. So, a DC level shifting is needed in the case of P polarity. For this purpose, a T resistor network is added in parallel with Rf and it is selected to read P side signals (Cf. Fig. 2). Its equivalent resistor value is given by:

t network f R R R R R R / 2 1 2 1 ) ( + + = . ( 6
)
Due to high output swing, the CSA does not include a source follower needed to increase the driving capability. This is not an issue for the complete readout chain (SCA + Shaper). However, the output signal of the CSA tested separately is slowed down and the measured gain is drastically reduced by the effect of the output load of the lines and pads.

B. The shaper implementation

A CR-RC signal shaping has been chosen to implement the two shapers. The slow shaper is directly connected to the output of the CSA while the fast one is connected through a buffering stage (Cf. Fig. 2). Active shaping is performed by using a two stage compensated Miller OTA (Operational Transconductance Amplifier). The shaping stages perform a pulse shaping primarily to optimize the signal-to-noise ratio (SNR) of the system and also to measure the energy and the time with the slow shaper around 1 µs shaping time and the fast one with a shaping time of 15 ns, respectively. Fig. 4 shows the CR-RC shaper and emphasizes the OTA architecture. The non-filtered noise of the shaper can be an issue in some cases as detailed in section III. CR-RC shaper transfer function in the frequency domain is given by:

( ) 2 1 ) ( s s s H τ τ + = , ( 7
)
where, τ is the shaping time. Here, we chose the derivation and the integration time constants to be equal:

2 2 1 1 C R R C shaping = = τ . ( 8
)
In time domain, the gain is given by:

e R R G SH 1 2 = . ( 9
)

III. NOISE STUDY

The most critical issue in this design is the noise. The desired energy resolution of the Compton camera which gives an acceptable spatial resolution of the whole hadrontherapy system must be below 1 keV. This corresponds to an input ENC (Equivalent Noise Charge) of 300 electrons rms. Fig. 5 illustrates the dominant noise sources. Generally, there are two main sources of noise in a CSA: a parallel and a series noise. The parallel noise is generated by the thermal noise of the feedback resistor and the detector leakage current fluctuations. The input referred power spectral density (PSD) of the parallel noise is given by:

2 4 2 ) ( n f G i i R kT qI S = + = ω , (10)
where, q is the electron elementary charge, I G is the detector leakage current, k is the Boltzmann constant and T is the temperature.

The noise input voltage spectrum (series noise PSD) of a MOS transistor in saturation is given by:

f A e WLf C K g kT S n ox f m v + = + = 2 3 8 ) (ω . ( 11
)
The first term is the channel thermal noise and the second term is the 1/f or Flicker noise [START_REF] Sansen | Limits of low noise performance of detector readout front ends in CMOS technology[END_REF]. Optimal noise matching is found by varying input transistor dimensions. Total (series and parallel) noise output voltage spectrum is given by:

2 2 2 2 2 1 ) ( f n n f d out C i f A e C C S ω ω + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = . ( 12
)
This formula is obtained through the following approximations: infinite open loop gain of CSA and considering an ideal integrator (i.e. no feedback resistance).

The output noise voltage of the complete chain is performed by integrating the product of CSA's output PSD by the shaper's Transfer function as in:

π ω ω ω 2 ) ( ) ( 2 0 2 d j H S V out out ∫ ∞ = . ( 13
)
Total input referred noise after shaping is given by: π τ τ where C t is the total input capacitance composed mainly by C d and grid-source C gs capacitance of the CSA input transistor, τ is the filtering time constant (shaping time). In this design, calculations were done to find the optimal value of W/L of the input transistor leading to the best trade-off between minimizing series and 1/f noise and maintaining low value of Ct. According to the formula (14), the dominant noise type depends on τ, C t and the value R f imposed by the counting rate of the physics experience. TABLE I shows a noise simulation result of the readout chain. One can see clearly that the noise due to the feedback resistor (here R f is split into 4 resistors for layout matching) represents 75 % of the total noise. The noise of the shaper can also be a key issue in case of low gain of the shaper given by the R2/R1 ratio (Fig. 4). As we have a wide dynamic range at the output of the CSA, this ratio is fixed to 4. The global theoretical gain of the is:

47 . 1 1 2 = = e R R G SH . ( 15
)
The current noise in the channel of transistor P4 (Fig. 4) also generates noise in the gate through the gate-channel capacitance C ox WL. The gate noise is due to the capacitive coupling frequency depending (f 2 ). The gate noise is about

2 2 5 4 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = T m g f f kTg i , ( 16
) gs m T C g f π 2 = , ( 17
)
where f T is the cut off frequency of the MOSFET. This noise is important at high frequencies and is not filtered by the shaper.

IV. LAYOUT DESIGN

The multi-channel nature of the design which will fit the form of the strips of the SSD detector (crosstalk issue), low noise and high speed are parameters that require a special care when drawing the layout of the input stage (CSA). Therefore, the R in must be minimized by splitting the input transistor into a matrix of small transistors perfectly identical. The size of each channel has exactly the width of a pad (100 µm). This makes the channel multiplication easier and minimizes the mismatch between channels. Fig. 6 shows the layout of the ASIC, an emphasis of the input transistor is carried out.

V. EXPERIMENTAL RESULTS

A prototype has been fabricated in a CMOS 0.35μm process from Austria Micro System. It includes 8 channels of the described FEE readout of the Compton camera. A photography of the chip is shown in Fig. 7. The circuit dissipates a power of 21 mW per channel from 3.3 V supply voltage and occupies an area of 2871×1881 μm 2 , including pads.

A. Test setup

A test board was designed to characterize the prototype experimentally. Fig. 8 shows a photography of the test board and a synoptic of its modules. A charge injection circuit is used to simulate the detector charge which will be generated by the stopping of a Compton recoil electron in the SSD. This circuit is composed of a fast pulse generator terminated by a 50 Ω resistor, a 1 pF injection capacitor and a15 pF capacitor to simulate the detector capacitance. The response of the ASIC to an input pulse is sent to a 50 Ω buffering stage to drive different signal processing devices (e.g. scope, external amplifiers…).

B. Test results and circuit performances

The circuit has been successfully tested. The test results are in good agreement with analytic and simulation calculations. All functionalities have been validated. Fig. 9 shows outputs waveforms of the circuit for both negative and positive signal polarities. The measured shaping time of the slow shaper is 0.7 µs. However, that of the fast one (35 ns) is greater than the calculated one (15 ns). This is due to the convolution between the proper shaping time and the peaking time of the CSA (30 ns). The circuit achieves a high linearity even over the needed dynamic range, mostly less than 1 % of INL. A conversion gain of 3.6 mV/fC, very close to post layout simulation (3.9 mV/fC), is obtained up to 600 fC. The noise evaluation has been performed according to the measurement standards. Noise and signal were measured in the same conditions (temperature and electromagnetic environment). The input referred noise ENC is obtained by:

out inj rms qV Q V ENC = , (18)
where Q inj is the injected charge, q is the electron charge, V out is the output voltage of the circuit corresponding to Q inj and V rms is the Standard Deviation which is given by:

) ( ) ( 2 2 setup asic setup Vr ms σ σ σ - + = = . ( 19
)
An ENC of 412 electrons rms was measured on the channel number 1 for electron configuration. 75% of this amount comes from the combination of the feedback resistor R f (4 MΩ) chosen to avoid the pile-up and the value of the shaping time (1µs). An ENC of 657 electrons rms is measured for the hole working mode. The performances obtained for hole configuration are slightly lower than those achieved for the electron one. This is due to the feedback network equivalent resistor used to shift the DC level in case of P polarity (Fig. 2). No significant dispersion of the ENC between channels is relevant. The ENC is plotted of each channel in Fig. 11. One can see the difference between the two configurations (electron and hole).

The noise which depends on the detector capacitance was also quantified. Fig. 12 plots the calculated, schematic and postlayout simulation as well as tested curves of the ENC versus detector capacitance. The slope of the theoretical curve is very small as no additional component has been taken into account in the calculations. The slope of the tested curve is slightly higher than that of simulation. It is due to additional external noisy components such as bias resistors. The parallel configuration and the proximity of channels give a cross talk of ~2 % (electron mode), mostly (≥99%) coming from the supply voltage disturbance. In Fig. 13, we injected a charge (in the left side) alternatively at the input of each single channel, and each time we measured the crosstalk on the seven other channels (front side of the figure) from the nearest to the farthest channel. The dispersion noticed on the crosstalk depends on which channel we injected the signal. It can be caused by mismatch between injection capacitors. The same crosstalk was measured for positive polarity (hole mode).

After electrical characterization, the ASIC was coupled to a 8 strip double-sided SSD detector to measure the energy spectrum of a barium radioactive source ( 133 Ba), providing photon lines at 30.9, 53, 80.9, 276, 303, 356 and 383 keV. Fig. 14 plots the energy spectrum measured on the N strips and P strips. Note that all P strips are connected together. The measurements were performed using the ASIC on one side and using a classical electronics on the other side of the detector. An energy resolution of 8 keV FWHM has been measured for both systems. This preliminary result does not reflect the real noise of the readout electronics and the detector because of the non-optimal test environment. More precisely this is due to the length of the cables connecting strips to FEE electronics. This issue will be solved by mounting the detector and the ASIC on the same test board and thereby minimizing the length of connections and their capacitances.

C. CSA + fast shaper chain performances

The fast shaper's branch of the chain has also been tested successfully. The test of this branch is more specific. Indeed, the charge collection time (arrival time of the signal charge at the CSA input) varies from 30 ns to 130 ns. Therefore, contrary to the slow shaper's branch which amplitude does not depend on the charge collection time, the amplitude of the fast shaper's branch depends on the collection time. This means that with the same charge, we can have different conversion gains, so making noise measurements more complicated. The area delimited by the fast shaper's output curve and the base line should be the same for a given charge, independently on the arrival time of the signal. This criterion can be used to calculate the charge at the input of the chain (CSA) after calibration. Such a measurement will be done in the next version of this circuit. A particular care must be taken of the slew rate when designing the OTA of the fast shaper. In Fig. 15 TABLE II gives a summary of the circuit performances. The measured parameters are close to schematic and post layout simulation ones. The exception of the CSA block is due the absence of a source follower stage on its output. Therefore, the gain is lowered in the test by the output load and consequently the measured input referred noise is increased. An indirect measurement of the CSA parameters using a complete chain and the gain of the slow shaper is carried out and gives results close to those of simulations. In order to read out the signals delivered by the Silicon Strip Detectors constituting the Compton camera for promptgamma monitoring of hadrontherapy, a 8 channel front end readout electronics has been designed and fabricated in 0.35 µm CMOS process from AMS. Each channel is built of a Charge Sensitive Amplifier (CSA) connected to two parallel shapers slow and fast. Exhaustive tests were performed and all the functionalities of the above-described circuit have been validated. The series noise of this design was minimized to ~5% of the total noise thanks to high optimization of the input transistor. In the future version of this circuit, we will replace the feedback resistor by a MOS switch and thus eliminate nearly 75% of the total noise. Also, further functionalities will be integrated such as a comparator providing a logic level for time stamp, multiple shaping and multiple gain.

A

  Low Noise and High Dynamic Charge Sensitive Amplifier-Shaper associated with Silicon Strip Detector for Compton Camera in hadrontherapy M. Dahoumane, D. Dauvergne, J. Krimmer, H. Mathez, C. Ray, E. Testa, A.H. Walenta, Y. Zoccarato

Fig. 1 .

 1 Fig. 1. Monitoring configuration system: the prompt γ-ray emission points are reconstructed by intersecting the ion trajectory, given by the hodoscope and the Compton cone, reconstructed with the Compton camera [4].

Fig. 2 .

 2 Fig. 2. General overview of the readout electronics coupled with a double sided SSD detector. The 8 strips of the SSD are AC-coupled with the 8 parallel channels of the circuit. The P and N strips are read with two different selectable configurations.

Fig. 3 .

 3 Fig. 3. The CSA architecture including RC feedback network.

Fig. 4 .

 4 Fig. 4. Shaper scheme on the top and Miller architecture of the OTA designed to perform the shaping function in the bottom of the figure.

Fig. 5 .

 5 Fig. 5. Main Noise sources of a CSA.

Fig. 6 .

 6 Fig. 6. Layout of the full chip fabricated in CMOS 0.35 µm process of AMS: a zoom on the CSA stage and its input transistor is shown.

Fig. 7 .

 7 Fig. 7. Photography of the chip showing its various blocks: 8 parallel complete channels, bias blocks and individual testing blocks.

Fig. 8 .

 8 Fig. 8. Test set up and photo of the test board supporting the ASIC (device under test). On the top, is shown a synoptic of the test board.

Fig. 9 .

 9 Fig. 9. Screenshot of the response of a complete chain to P and N polarities. In (a) is shown the CSA+slow shaper and in (b) CSA+fast shaper. Linearity measurement was also performed. The top of Fig. 10 shows the measured output and the fitted curves of the readout chain (with slow shaper) as a function of the input charge (in fC). The bottom of this figure shows the INL (Integral Non Linearity) for the two working configurations (electron on the left and hole on the right side).

Fig. 10 .

 10 Fig. 10. Output vs. injected charge (Qinj) response of a complete chain for P and N polarities on the top side. INL vs. Qinj in the bottom sideof the figure.

Fig. 11 .

 11 Fig. 11. ENC vs. input charge measured on all channels for electron and hole working modes.

Fig. 12 .

 12 Fig. 12. Calculation, scheme and layout simulation and test ENC vs. detector capacitance.

Fig. 13 .

 13 Fig. 13. Crosstalk between the CSA+Slow shaper channels measured for the electron configuration.

Fig. 14 .

 14 Fig. 14. Energy spectrum of the 133 Ba source measured with discrete electronics and the described ASIC connected to a 2mm thick 8-strip SSD (11x11 mm2 active area). P strips are connected together.

  , we highlight this phenomenon. One can see the slew rate issue before enhancing the slew rate (left side of the figure); this issue disappears by enhancing the slew rate (right side of the figure). Indeed, a high slew rate is needed to drive the output load for fast shape and high dynamic signals (~35 ns shaping time in case of the fastest signal). This issue is present only in case of falling signals which turn the output transistor of the shaper's OTA out of saturation region.

Fig. 15 .

 15 Fig. 15. Slew rate issue in case of negative, high amplitude and fast signal.

TABLE I .

 I NOISE SIMULATION SUMMARY: CONTRIBUTION OF VARIOUS COMPONENTS IN %

TABLE II .

 II PERFORMANCE SUMMARY VI. CONCLUSION AND UPCOMING WORK

ACKNOWLEDGMENT First, we would thank Mr. E. Bechetoille for his very useful technical help. We would also like to thank Mr. D. Delaunay and Mr. R. Della Negra for their contribution on the test of this circuit.


/+'%&'!-"++'3)' 9 :!4.)3)44'! 8+*%'

We report advances in understanding how femtosecond laser ripples formation in metal works. The surface topology features are discussed in terms of periodicity and contrast of the pattern formation, related to chronological laser-induced events. Resonant excitation of Surface Plasmons (SP) in metallic gratings show that SP wave excited during the femtosecond laser pulse can initiate the observed patterning and play a key role. Metals behavior under nonequilibrium conditions on the picosecond timescale is then investigated to correlate the amount of material experiencing solid-to-liquid transitions and the subsequent ripples contrast. With the de-rived observation, a calculation of the transient nonequilibrium thermodynamic characteristics of excited nickel is also performed allowing defining characteristic timescales of thermocapillary proc-ess which may occur under multipulse irradiation.

1-Introduction

Application potential for ultrashort laser pulses relies on their capacity to confine and localize energy on the smallest scales. Despite diffraction limit, in particular conditions, ultrafast laser excitation can structure well below the incident wavelength. In particular, it can trigger nanoscale self-arrangement of matter on irradiated surfaces [1]- [7]. Appearing for fluences under the ablation threshold and developing with varying periodicity and amplitudes, the manifestation of periodic surface patterns with subwavelength dimensions, usually known as Laser-Induced Periodic Surface Structure (LIPSS) or « ripples », is a not well-understood but inherent phenomenon. Two distinct types of LIPSS are typically observed when solid materials are irradiated upon ultrashort multiple pulses: Low Spatial Frequency LIPSS (LSFL) and High Spatial Frequency LIPSS (HSFL) [8]- [12]. The microscopic mechanism of ripples formation is still under discussion and various kinds of surface wave excitation methods have been proposed so far besides the conventional interference pattern of the incident or refracted light and scattered wave [13]- [17]. Whereas ultrafast laser excitation stage is the most investi-gated mechanisms in recent literature, an emphasis has to be put on the surface dynamics in this particular irradiation regime to explain the subsequent spontaneous formation of nanostructured patterns. Therefore, we describe here several approaches to shine light on different aspects. A comprehension effort related to laser excitation mechanisms will be developed, able to assist and validate experiments performed on metals. In this context, we investigate hypotheses assuming that LSFL, having a spatial period close to the irradiation wavelength, are formed due to optical interference of the incident laser radiation with a surface plasmon. HSFL formation mecha-nism is still quite controversial for reasons related to its polarization dependence and it will not be discussed in this paper. If the most part of experiments are performed on a practically smooth surface, we will present irradiation re-sults obtained on pre-structured surfaces. Grating-coupled surface plasmon resonance on metallic surface has been exploited to investigate the correlation between ripples formation under ultrashort laser exposure and SP genera-tion conditions. For Transverse Magnetic (TM) irradiation conditions, a well-defined period of grating was a seed for low spatial frequency LIPSS formation [18]. To design matter transformation patterns with predictable properties for a wide range of materials, the contrast of the nanostruc-turing process is another key parameter [5], [19]- [20]. Fol-lowing localized photoexcitation, the sub-surface region undergoes transient melting and ultra-rapid quenching at the nanoscale. A thorough understanding of the involved physical mechanisms during the matter transformation will therefore help to control LIPSS morphology. Measuring the contrast of structures obtained on nickel and copper sur-faces which are supposed to have different behavior under nonequilibrium conditions, the effects of energy coupling and relaxation strength have been correlated to LIPSS for-mation [20]. The amount of material experiencing solid-to-liquid transitions appears to be correlated to the ripples height. Finally, the characteristic development time of LIPSS formation will be discussed in the context of matter transformation phases using an hydrodynamic approach. Simulations of laser-matter interaction coupling hydrody-namics and heat transfer will be presented to estimate laser energy deposition, electron-ion nonequilibrium stage, ma-terial heating and growth. The subsequent phase transfor-mation and the lifetime of molten region can have implica-tions on thermocapillarity phenomena. This will be effec-tively estimated by applying characteristic timescales fram-ing the positive feedback of ripple growth under multipulse irradiation. Finally, the depth of the formed pattern will also affect the grating coupling varying ripples periodicity via electromagnetic feedback effects.

Plasmon wavelength defining LIPSS periodicity

As LIPSS performed on absorbing materials such as metals or excited semiconductors have been mostly found to grow along the direction perpendicular to the laser polarization, LSFL has been early attributed to the excitation of surface plasmon polaritons [8], [21]. A photon cannot transform into a SP without a coupling system (i.e a prism or a grating). Since a rough surface corresponds to a superposition of a number of different gratings, the mechanism underlying the coupling of SP on a rough interface is essentially the same as that via a single corrugation [22]. It is well known that in the specific case of a resonance harmonic grating, the incident laser wave impinging on the grating at a certain angle transforms into a diffracted wave propagating along the grating surface as an SP wave.. The proposed method allows correlating the excitation of plasmons in the optical domain to the formation of ripples. We have investigated the formation of LIPSS on pre-structured nickel surfaces having periodic gratings with a prechosen range of periodicities upon irradiation by single linearly polarized Ti:sapphire laser pulses (800 nm, 150 fs) in air as a function of the grating period and orientation of the polarization [18]. Considering only the first diffraction order, n = 1, the condition of phase-matching simplifies to a SP wave with a wavelength equal to the grating spacing. For the experimental part, a Mach-Zehnder interference scheme has been used to create a number of corrugation gratings on nickel plates coated with a thin photoresist deposited on a Ni substrate by spin coating. Different grating periods G ranging from 440 nm to 800 nm, with 10 nm increment, have been fabricated. The struc-tures surfaces were irradiated by a single fs laser pulse. Fig. 1a) shows an example of periodic nanorelief written on a Ni sample before irradiation.

From the complex optical index of Ni, 2.48+i4.38 given in Ref. [START_REF]Handbook of optical constants of solids: ED Palik Academic[END_REF], we estimate that the maximal coupling should occur at G = 790 nm. To show that the grating of 10 nm thickness does not modify the expected resonant period, we have calculated the 0th-order reflection efficiency. This calculation is performed for an incident TM wave at 800 nm irradiating a Ni surface grating of 10 nm depth by Rigorous Coupled Wave Analysis (RCWA) [START_REF] Lyndin | [END_REF]. A rectangular profile of the grooves has been assumed, according to AFM measurements of the gratings. Fig. 1b) illustrates the 0th-order reflection efficiency response; it can be seen that the maximal resonant coupling is close to the expected value of 790 nm, for this relatively low grating depth. Consequently, it can be assumed that the resonant position is not affected by this relatively small grating depth as it will be discussed in more details in section 4.

Figure 1: Grating-coupled surface plasmon enhanced excitation. Example of a representative pre-structured surface at 750 nm (a). 0th order reflection efficiency calculated by RCWA show-ing a resonance expected for nm (b). SEM image of the LIPSS formation observed for the specific grating period nm (c) and corresponding 2D-FT image (d). Radially averaged autocorrelation function R(r) processed from the selection on the SEM image (e).

The irradiation of each grating by a single pulse at a fluence of 1.42 J/cm2 for the two polarizations at normal incidence reveals that ripples formation only occurs for TM irradiation and a well defined grating period of G = 750 nm. Fig. 1c) shows the LIPSS structure generated for this grating periodicity. Single pulse irradiation allows achieving only the early stage of ripple formation and we observe that the pattern is not homogeneously distributed in the irradiated spot. 2D Fourier transform (2D-FT) on the selected part of the acquired image is displayed in Fig. 1d), showing a vertical sickle-shaped feature confirming that ripples are perpendicular to the laser polarization as ex-pected for LSFL on metals. The SEM image and the corre-sponding 2D-FT reveal that the period of the formed ripples is 750 nm. This result indicates that the specific grating period for which SP mode is efficiently excited is correlated to the ripples structure. The observed resonance is however different from the expected one around 790 nm. Hence, it is supposed that the ultrafast material excitation aects the dielectric properties during the laser pulse, modifies the plasmonic dispersion relation and thus results in a dierent grating periodicity for which resonance occurs. A radially averaged autocorrelation function R(r) has also been calculated in the direct-space to take into account the k vectors concentrated along the arcs in a (kx, ky) plane. R(r) is calculated from the autocorrelation of pixels of the binary SEM image as a function of the radial distance r [START_REF]Radially Averaged Autocorrelation Function has been processed with ImageJ, a public domain java-based image processing software[END_REF]. The calculated function is displayed on Fig. 1e). The oscillating behavior of R(r) is characteristic of a long-range order with periodicity 750 nm, which is more difficult to observe in the SEM image in Fig. 1(c). The decay of the oscillation amplitude is mainly due to features such as ir-regular orientation, bifurcation and limited periodic structures. However, the correlation of the LIPSS appears to extend over several periods. The observed spatial periods of the LIPSS, somewhat smaller than the expected one, can be explained by transient changes of the dielectric proper-ties of the laser-excited material [8], [18], [22]. To illustrate this, we have calculated SP wavelength for a range of optical index. The required optical index allowing to reach SP = G 750 nm is ñ=0.9+i2. The optical index corresponding to an unexcited material is ñ=2.48+i4.38. Several optical index are possible to reproduce this SP period.

Hydrodynamic stages

The plasmonic approach presented in the first section establishes preliminary steps defining orientation and energy modulation guiding the mechanisms of LIPSS formation. Once the process is triggered, further steps towards nanostructuration follow and several processes are responsible for the final morphology. They require extended approach defining the role of material response from a thermodynamic and hydrodynamic perspective. In particular, material behaviour is strongly dependent on the electronic response upon ultrafast excitation. The effect of the laser irradiation with the response of the material is commonly described by the two-temperature model (TTM).

The energy balance between the electron cloud and the lattice is described by the temporal evolution of the electron and the lattice temperatures through a coupling parameter related to the electron-phonon scattering rate and by thermal diffusion processes occurring within each subsystems. The presence of a transient nonequilibrium stage has strong implications on the energy distribution and thus on the final morphology. As energy relaxation dynamics are material dependent, we generated LIPSS on representative materials with different electronic configuration in order to investigate the influence of and thermal diffusion efficiency on the specific contrast of LIPSS. Metals with a large range of coupling strengths were exposed to ultrafast laser excitation and the patterning surface has been investigated exsitu by atomic force microscopy, and scanning electron microscopy [20].

The experimental results were obtained for Ni and Cu irradiated by 50 pulses of 80 J (150 fs) at 800 nm. The focal spot obtained in the case of nickel appears much wider than for copper due to different absorption coefficient. Nevertheless, absorbed fluence allows to generate LIPSS at the periphery of the focal spot.The relative strength and the electron temperature (Te) dependence of the electron-phonon coupling in the metals were identified as key factors affecting the initial energy redistribution and the dimensional size of the region undergoing transient melting and resolidification. It appears that growth of the contrast and increase of the ripples height are correlated to the evolution of the constant of coupling and to the amount of material experiencing solid-to-liquid transitions.

The amplitudes varies as a function of the material since Ni has an averaged peak-to-peak amplitude H200 nm while Cu presents a nearly two fold lower contrast H100 nm. These two materials differ mainly by their electron-We have seen that LIPSS are readily formed upon ul-trafast excitation on metals by laser pulses with a fluence sufficient to induce phase transformation. Following sin-gle-pulse exposure, a quasi-amorphous resolidified layer of several tens of nm in thickness has been experimentally observed [29]- [30]. These results suggest that the freezing structures morphology can depend on the thermodynamic properties of the irradiated materials. To investigate the transient thermodynamic stages of the metal surface, we have performed the simulation of the nonequilibrium heating of a nickel sample under laser irradiation using a one dimensional two temperature hydrodynamic code (Esther) [31]- [32]. The approach solves, according to a Lagrangian scheme, the fluid equations for the conservation of mass, momentum and energy for electronic and ionic species. The interaction between the laser field and the target is calculated by solving the Helmholtz wave equation in the inhomogeneous media. During relaxation, the material thermodynamic properties are described by the Bushman-Lomonosov-Fortov multiphase equation of state covering a large range of densities and temperatures between con-densed states and hot plasmas [33]. This calculation en-ables us to monitor thermodynamic transformations of the surface, the energy profile, the onset of phase transition and hydrodynamic movement in the longitudinal direction (cor-responding to laser propagation in normal incidence) to correlate them with the expected evolution of the surface topology. As a result, the evolution of thermodynamic con-ditions experienced by the surface is determined and de-picted on Fig. 4. The figure indicates time-resolved tem-perature [Fig. 4(a)] and density [Fig. 4(b)] profiles of the heated surface under a 150 fs pulse laser irradiation at 800 nm for a fluence regime allowing exceeding melting tem-perature at the surface but preventing an ablation process (0.35 J/cm2 in the simulation). As a fact, ripples are sup-posed to appear only when a sufficient energy is deposited at the surface to induce liquid phase transformation and material flow, whereas for higher fluence vaporization and ablation of the surface will loosen the modulated energy deposition then inhibiting LIPSS formation. The simula-tions show also the kinetics of melting and solidification under ultrafast excitation. A time of a few picoseconds is required to reach the melting temperature, with the rise slope being related to the electronic temperature depend-ence of the electron-phonon coupling. Heat conduction to the bulk of the sample leads to cooling below the melting temperature, causing resolidification on several hundreds of picoseconds timescale at the surface. The averaged ve-locity of solidification can be estimated at around 80 m.s-1, close to the one observed in experiments performed on Ag thin films heated by a fs excitation [34]. Depth of melting evolves as a function of time and reaches a maximal value of 85 nm at 17 ps. This time corresponds to the characteris-tic value for acoustic propagation at the sound velocity. Therefore, the depth of the molten zone decreases after shock relaxation but electronic thermal conduction contin-ues to favorize the increase the thickness of the melted layer in the following tens of picoseconds, up to around one hundred of picoseconds. We note that the size of the liquid layer is in agreement with the estimated electron diffusion depth in the solid Ni at the melting threshold since electronic thermal diffusion defines the heated depth in a nonequilibrium system [20]. The time scales of melting ( ps) and resolidification ( ps) as well as the melting depth ( nm) estimated in this simulation, are consistent with other calculations [35] and with most experimental observations [35]- [39].

Feedback process

The above calculation corresponds to single pulse con-ditions. However in multipulse exposure conditions, the morphology of the final irradiation differs strongly from the observed single-shot one. A feedback is therefore invoked. The feedback that causes coherent structures to grow can originate from different mechanisms. Ripples growth requires a positive feedback and we discuss in this last section possible feedback processes related to surface tension effects and plasmonic resonance alterations which can lead to non-trivial structure development.

1D-hydrodynamic calculation restricts to study in-depth effects whereas due to the large surface-to-volume ratio at nanoscale, the surface tension is supposed to play an im-portant role in lateral organization [40]. A possible scenario of surface motion in a hydrodynamic process is that strong gradient of surface tension can cause a thermocapillary force (also referred to Marangoni effect or Bénard-Marangoni convection) since decreases as the equilibrium temperature (referred as T) increases [41]. A gradient of surface tension occurs being driven by a temperature gradi-ent T which will force the liquid to flow away from regions of low surface tension. Therefore, thermal gradients can drive the liquid flow from hotter to colder regions leading to a modulated material transport on mesoscopic scales. The time required for one convective cycle should be smaller than the total laser dwell time in terms of melt life-time [40], [START_REF] Bäuerle | Laser Processing and Chemistry[END_REF]. At least several pulses (N) can be required to induce a significant contrast of the induced structure It results from previous considerations that plasmonic stage does not necessarily select the dominant modulation of the final pattern since thermocapillary waves can also contribute to a nonlinear growth of specific periodicities. Never-theless, feedback process induced by several laser shots can involve both plasmonic and capillary processes since SP are surface modes strongly localized at the surface of the metal, and therefore they are sensitive to any changes that occur in the vicinity of the air/metal interface. Indeed, the approach described in section 2 is based on a plasmonic response serving as the starting inhomogeneous absorption of the surface. This mechanism is supposed to determine both the orientation of the structures, due to a polarization dependency, and an initial periodicity which should evolve in time since LIPSS contrast increases as the number of pulses N increases. Hence, surface modification results in an increasing grating depth which, in turn, changes the surface wave dispersion. This surface modification shifts the position of the resonant grating period, where the amplitude of the shift depends on the LIPSS contrast. Figure 2 illustrates the zero-th order reflexion efficiency evolution from a grating with varying depth modulation from 10 nm to 100 nm. The calculation has been performed for a TM wave at 800 nm incident on a rectangular surface grating of nickel using RCWA [START_REF] Lyndin | [END_REF]. As it has been already shown [46], the SP resonance shifts toward a smaller periodicity for deeper gratings and a broadening of the resonance is observed. At this point, we note that a grating depth of about 75 nm is required to re-produce the resonance observed in our experiment at 750 nm and that the shift discussed in section 2 cannot be at-tributed to an undesired pertubative effect. For stuctures deeper than 30 nm, a large laser absorption is expected, accelerating the structuring process and promoting positive cumulative effects for a specific wavelength. Coherent structures are then supposed to develop with a dominant period related to the most enhancing or the most rapidly growing period given by the optimal multi-pulse electro-magnetic and hydrodynamic feedback.

Conclusion

In conclusion we have indicated some fundamental as-pects deriving from the plasmonic coupling, the electronphonon nonequilibrium stage, the subsequent sequence of phase transformation and finally the surface dynamics via surface tension gradients. Grating-coupled Surface Plas-mon Polaritons on metallic surface has been exploited to investigate the correlation between ripples formation under ultrashort laser exposure and SP g e n e r a t i o n conditions. The main finding of this experiment was an experimental dem-onstration of the resonant coupling of the fs-laser pulse to the grating structure, based on SP excitation, and its impact on the early stage of LIPSS formation. The understanding formation dynamics of regular self-organized nanoscale structures under laser exposure of materials still required rigorous concepts related to modulated energy deposition, laser-induced nonequilibrium phenomena and the related relaxation paths involving mass transport. Nevertheless we have shown that in the particular case of strongly absorbing materials which can undergo capillary action as metals, hydrodynamic modeling approaches should be involved to explain the observed patterns. Extended understanding ef-fort has to be made to check if optically-induced modula-tion can be preserved when the system undergoes phase transitions and how feedback mechanisms can cause the growth of grating periodicities by coupling to the diffrac-tive laser light, involving plasmonic/hydrodynamic effects.

Figure 2 :

 2 Figure 2 : Calculation of 0th-order reflection efficiency as a func-tion of grating period for different grating depth.


Rayan Chikhi

Guillaume Rizk

Space-efficient and exact de Bruijn graph representation based on a Bloom filter

The de Bruijn graph data structure is widely used in nextgeneration sequencing (NGS). Many programs, e.g. de novo assemblers, rely on in-memory representation of this graph. However, current techniques for representing the de Bruijn graph of a human genome require a large amount of memory (≥ 30 GB). We propose a new encoding of the de Bruijn graph, which occupies an order of magnitude less space than current representations. The encoding is based on a Bloom filter, with an additional structure to remove critical false positives. An assembly software implementing this structure, Minia, performed a complete de novo assembly of human genome short reads using 5.7 GB of memory in 23 hours.

Introduction

The de Bruijn graph of a set of DNA or RNA sequences is a data structure which plays an increasingly important role in next-generation sequencing applications. It was first introduced to perform de novo assembly of DNA sequences [START_REF] Idury | A new algorithm for DNA sequence assembly[END_REF]. It has recently been used in a wider set of applications: de novo mRNA [START_REF] Grabherr | Full-length transcriptome assembly from RNA-Seq data without a reference genome[END_REF] and metagenome [START_REF] Peng | Meta-IDBA: a de novo assembler for metagenomic data[END_REF] assembly, genomic variants detection [START_REF] Peterlongo | Identifying SNPs without a reference genome by comparing raw reads[END_REF][START_REF] Iqbal | De novo assembly and genotyping of variants using colored de bruijn graphs[END_REF] and de novo alternative splicing calling [START_REF] Sacomoto | KISSPLICE: de-novo calling alternative splicing events from RNA-seq data[END_REF]. However, an important practical issue of this structure is its high memory footprint for large organisms. For instance, the straightforward encoding of the de Bruijn graph for the human genome (n ≈ 2.4 • 10 9 , k-mer size k = 27) requires 15 GB (n • k/4 bytes) of memory to store the nodes sequences alone. Graphs for much larger genomes and metagenomes cannot be constructed on a typical lab cluster, because of the prohibitive memory usage.

Recent research on de Bruijn graphs has been targeted on designing more lightweight data structures. Li et al. pioneered minimum-information de Bruijn graphs, by not recording read locations and paired-end information [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF]. Simpson et al. implemented a distributed de Bruijn graph to reduce the memory usage per node [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF]. Conway and Bromage applied sparse bit array structures to store an implicit, immutable graph representation [START_REF] Conway | Succinct data structures for assembling large genomes[END_REF]. Targeted methods compute local assemblies around sequences of interest, using negligible memory, with greedy extensions [START_REF] Warren | Targeted assembly of short sequence reads[END_REF] or portions of the de Bruijn graph [START_REF] Peterlongo | Mapsembler, targeted and micro assembly of large NGS datasets on a desktop computer[END_REF]. Ye et al. recently showed that a graph roughly equivalent to the de Bruijn graph can be obtained by storing only one out of g nodes (10 ≤ g ≤ 25) [START_REF] Ye | Exploiting sparseness in de novo genome assembly[END_REF]. Conway and Bromage observed that the self-information of the edges is a lower bound for exactly encoding the de Bruijn graph [START_REF] Conway | Succinct data structures for assembling large genomes[END_REF]:

log 2 ( 4 k+1 |E| ) bits,
where k + 1 is the length of the sequence that uniquely defines an edge, and |E| is the number of edges. In this article, we will consider for simplicity that a de Bruijn graph is fully defined by its nodes. A similar lower bound can then be derived from the self-information of the nodes. For a human genome graph, the self-information of |N | ≈ 2.4 • 10 9 nodes is log 2 ( 4 k |N | ) ≈ 6.8 GB for k = 27, i.e. ≈ 24 bits per node.

A very recent article [START_REF] Pell | Scaling metagenome sequence assembly with probabilistic de bruijn graphs[END_REF] from Pell et al. introduced the probabilistic de Bruijn graph, which is a de Bruijn graph stored as a Bloom filter (described in the next section). It is shown that the graph can be encoded with as little as 4 bits per node. An important drawback of this representation is that the Bloom filter introduces false nodes and false branching. However, they observe that the global structure of the graph is approximately preserved, up to a certain false positive rate. Pell et al. did not perform assembly directly by traversing the probabilistic graph. Instead, they use the graph to partition the set of reads into smaller sets, which are then assembled in turns using a classical assembler. In the arXiv version of [START_REF] Pell | Scaling metagenome sequence assembly with probabilistic de bruijn graphs[END_REF] (Dec 2011), it is unclear how much memory is required by the partitioning algorithm.

In this article, we focus on encoding an exact representation of the de Bruijn graph that efficiently implements the following operations:

1. For any node, enumerate its neighbors 2. Sequentially enumerate all the nodes The first operation requires random access, hence is supported by a structure stored in memory. Specifically, we show in this article that a probabilistic de Bruijn graph can be used to perform the first operation exactly, by recording a set of troublesome false positives. The second operation can be done with sequential access to the list of nodes stored on disk. One highlight of our scheme is that the resulting memory usage is approximated by 1.44 log 2 ( 16k 2.08 ) + 2.08 bits/k-mer. For the human genome example above and k = 27, the size of the structure is 3.7 GB, i.e. 13.2 bits per node. This is effectively below the self-information of the nodes. While this may appear surprising, this structure does not store the precise set of nodes in memory. In fact, compared to a classical de Bruijn graph, the membership of an arbitrary node cannot be efficiently answered by this representation. However, for the purpose of many applications (e.g. assembly), these membership queries are not needed. We apply this representation to perform de novo assembly by traversing the graph. In our context, we refer by traversal to any algorithm which visits all the nodes of the graph exactly once (e.g. a depth-first search algorithm). Thus, a mechanism is needed to mark which nodes have already been visited. However, nodes of a probabilistic de Bruijn graph cannot store additional information. We show that recording only the visited complex nodes (those with in-degree or outdegree different than one) is a space-efficient solution. The combination of (i) the probabilistic de Bruijn graph along with the set of critical false positives, and (ii) the marking scheme, enables to perform very low-memory de novo assembly.

In the first Section, the notions of de Bruijn graphs and Bloom filters are formally defined. Section 3 describes our scheme for exactly encoding the de Bruijn graph using a Bloom filter. Section 4 presents a solution for traversing our representation of the de Bruijn graph. Section 6 presents two experimental results: (i) an evaluation of the usefulness of removing false positives and (ii) an assembly of a real human dataset using an implementation of the structure. A comparison is made with other recent assemblers based on de Bruijn graphs.

de Bruijn graphs and Bloom filters

The de Bruijn graph [START_REF] Idury | A new algorithm for DNA sequence assembly[END_REF], for a set of strings S, is a directed graph. For simplicity, we adopt a node-centric definition. The nodes are all the k-length substrings (also called k-mers) of each string in S. An edge s 1 → s 2 is present if the (k -1)-length suffix of s 1 is also a prefix of s 2 . Throughout this article, we will indifferently refer to a node and its k-mer sequence as the same object.

A more popular, edge-centric definition of de Bruijn graphs requires that edges reflect consecutive nodes. For k ′ -mer nodes, an edge s 1 → s 2 is present if there exists a (k ′ + 1)-mer in a string of S containing s 1 as a prefix and s 2 as a suffix. The node-centric and edge-centric definitions are essentially equivalent when k ′ = k -1 (although in the former, nodes have length k, and k -1 in the latter).

The Bloom filter [START_REF] Kirsch | Less hashing, same performance: Building a better bloom filter[END_REF] is a space efficient hash-based data structure, designed to test whether an element is in a set. It consists of a bit array of m bits, initialized with zeros, and h hash functions. To insert or test the membership of an element, h hash values are computed, yielding h array positions. The insert operation corresponds to setting all these positions to 1. The membership operation returns yes if and only if all of the bits at these positions are 1. A no answer means the element is definitely not in the set. A yes answer indicates that the element may or may not be in the set. Hence, the Bloom filter has one-sided errors. The probability of false positives increases with the number of elements inserted in the Bloom filter. When considering hash functions that yield equally likely positions in the bit array, and for large enough array size m and number of inserted elements n, the false positive rate F is [START_REF] Kirsch | Less hashing, same performance: Building a better bloom filter[END_REF]:

F ≈ 1 -e -hn/m h = 1 -e -h/r h ( 1
)
where r = m/n is the number of bits per element. For a fixed ratio r, minimizing Equation 1 yields the optimal number of hash functions h ≈ 0.7r, for which F is approximately 0.6185 r . Solving Equation 1 for m, assuming that h is the optimal number of hash function, yields m ≈ 1.44 log 2 ( 1 F )n. The probabilistic de Bruijn graph is obtained by inserting all the nodes of a de Bruijn graph (i.e all k-mers) in a Bloom filter [START_REF] Pell | Scaling metagenome sequence assembly with probabilistic de bruijn graphs[END_REF]. Edges are implicitly deduced by querying the Bloom filter for the membership of all possible extensions of a k-mer. Specifically, an extension of a k-mer v is the concatenation of either (i) the k -1 suffix of v with one of the four possible nucleotides, or (ii) one of the four nucleotides with the k -1 prefix of v. The probabilistic de Bruijn graph holds an over-approximation of the original de Bruijn graph. Querying the Bloom filter for the existence of an arbitrary node may return a false positive answer (but never a false negative). This introduces false branching between original and false positive nodes.

3 Removing critical false positives

The cF P structure

Our contribution is a mechanism that avoids false branching. Specifically, we propose to detect and store false positive elements which are responsible for false branching, in a separate structure. To this end, we introduce the cF P structure of critical False Positives k-mers, implemented with a standard set allowing fast membership test. Each query to the Bloom filter is modified such that the yes answer is returned if and only if the Bloom filter answers yes and the element is not in cF P .

Naturally, if cF P contained all the false positives elements, the benefits of using a Bloom filter for memory efficiency would be lost. The key observation is that the k-mers which will be queried when traversing the graph are not all possible k-mers. Let S be the set of true positive nodes, and E be the set of extensions of nodes from S. Assuming we only traverse the graph by starting from a node in S, false positives that do not belong to E will never be queried. Therefore, the set cF P will be a subset of E. Let P be the set of all elements of E for which the Bloom filter answers yes. The set of critical false positives cF P is then formally defined as cF P = P \ S.

Figure 1 shows a simple graph with the set S of correct nodes in regular circles and cF P in dashed rectangles. The exact representation of the graph is therefore made of two data structures: the Bloom filter, and the set cF P of critical false positives. The set cF P can be constructed using an algorithm that limits its memory usage, e.g. to the size of the Bloom filter. The set P is created on disk, from which cF P is then gradually constructed by iteratively filtering P with partitions of S loaded in a hash-table.

Dimensioning the Bloom filter for minimal memory usage

The set cF P grows with the number of false positives. To optimize memory usage, a trade-off between the sizes of the Bloom filter and cF P is studied here. Using the same notations as in the definition of the Bloom filter, given that n = |S|, the size of the filter m and the false positive rate F are related through Equation 1. The expected size of cF P is 8n • F, since each node only has eight possible extensions, which might be false positives. In the encoding of cF P , each k-mer occupies 2 • k bits. Recall that for a given false positive rate F, the expected optimal Bloom filter size is 1.44n log 2 ( 1 F ). The total structure size is thus expected to be

1.44n log 2 1 F Bloom filter + (16 • Fnk) cF P bits ( 2
)
The size is minimal for F ≈ (16k/2.08) -1 . Thus, the minimal number of bits required to store the Bloom filter and the set cF P is approximately

n • (1.44 log 2 ( 16k 2.08 ) + 2.08). (3)
For illustration, Figure 2-(a) shows the size of the structure for various Bloom filter sizes and k = 27. For this value of k, the optimal size of the Bloom filter is 11.1 bits per k-mer, and the total structure occupies 13.2 bits per k-mer. Figure 2-(b) shows that k has only a modest influence on the optimal structure size. Note that the size of the cF P structure is in fact independent of k.

In comparison, a Bloom filter with virtually no critical false positives would require F • 8n < 1, i.e. r > 1.44 log 2 (8n). For a human genome (n = 2.4 • 10 9 ), r would be greater than 49.2, yielding a Bloom filter of size 13.7 GB.

Additional marking structure for graph traversal

Many NGS applications, e.g. de novo assembly of genomes [START_REF] Miller | Assembly algorithms for next-generation sequencing data[END_REF] and transcriptomes [START_REF] Grabherr | Full-length transcriptome assembly from RNA-Seq data without a reference genome[END_REF], and de novo variant detection [START_REF] Sacomoto | KISSPLICE: de-novo calling alternative splicing events from RNA-seq data[END_REF], rely on (i) simplifying and (ii) traversing the de Bruijn graph. However, the graph as represented in the previous section neither supports (i) simplifications (as it is immutable) nor (ii) traversals (as the Bloom filter cannot store an additional visited bit per node). To address the former issue, we argue that the simplification step can be avoided by designing a slightly more complex traversal procedure [START_REF] Chikhi | Localized genome assembly from reads to scaffolds: practical traversal of the paired string graph[END_REF].

We introduce a novel, lightweight mechanism to record which portions of the graph have already been visited. The idea behind this mechanism is that not every node needs to be marked. Specifically, nodes that are inside simple paths (i.e nodes having an in-degree of 1 and an out-degree of 1) will either be all marked or all unmarked. We will refer to nodes having their in-degree or out-degree different to 1 as complex nodes. We propose to store marking information of complex nodes, by explicitly storing complex nodes in a separate hash table. In de Bruijn graphs of genomes, the complete set of nodes dwarfs the set of complex nodes, however the ratio depends on the genome complexity [START_REF] Kingsford | Assembly complexity of prokaryotic genomes using short reads[END_REF]. The memory usage of the marking structure is n c C, where n c is the number of complex nodes in the graph and C is the memory usage of each entry in the hash table (C ≈ 2k + 8).

Implementation

The de Bruijn graph structure described in this article was implemented in a new de novo assembly software: Minia3 . An important preliminary step is to retrieve the list of distinct k-mers that appear in the reads, i.e. true graph nodes. To discard likely sequencing errors, only the k-mers which appear at least d times are kept (solid k-mers). We experimentally set d to 3. Classical methods that retrieve solid k-mers are based on hash tables [START_REF] Marais | A fast, lock-free approach for efficient parallel counting of occurrences of k-mers[END_REF], and their memory usage scale linearly with the number of distinct k-mers. To avoid using more memory than the whole structure, we implemented a constant-memory k-mer counting procedure (manuscript in preparation). To deal with reverse-complementation, k-mers are identified to their reverse-complements.

We implemented in Minia a graph traversal algorithm that constructs a set of contigs (gap-less sequences). The Bloom filter and the cF P structure are used to determine neighbors of each node. The marking structure records already traversed nodes. A bounded-depth, bounded-breadth BFS algorithm (following Property 2 in [START_REF] Chikhi | Localized genome assembly from reads to scaffolds: practical traversal of the paired string graph[END_REF]) is performed to traverse short, locally complex regions. Specifically, the traversal ignores tips (dead-end paths) shorter than 2k + 1 nodes. It chooses a single path (consistently but arbitrarily), among all possible paths that traverse graph regions of breadth ≤ 20, provided these regions end with a single node of depth ≤ 500. These regions are assumed to be sequencing errors, short variants or short repetitions of length ≤ 500 bp. The breadth limit prevents combinatorial blowup. Note that paired-end reads information is not taken into account in this traversal. In a typical assembly pipeline (e.g. [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF]), a separate program (scaffolder ) can be used to link contigs using pairing information.

Results

Throughout the Results section, we will refer to the N50 metric of an assembly as the longest contig size, such that half the assembly is contained in contigs longer than this size.

On the usefulness of removing critical false positives

To test whether the combination of the Bloom filter and the cF P structure offers an advantage over a plain probabilistic de Bruijn graph, we compared both structures in terms of memory usage and assembly consistency. We retrieved 20 million E. coli short reads from the Short Read Archive (SRX000429), and discarded pairing information. Using this dataset, we constructed the probabilistic de Bruijn graph, the cF P structure, and marking structure, for various Bloom filter sizes (ranging from 5 to 19 bits per k-mer) and k = 23 (yielding 4.7 M solid k-mers).

We measured the memory usage of both structures. For each, we performed an assembly using Minia with exactly the same traversal procedure. The assemblies were compared to a reference assembly (using MUMmer), made with an exact graph. The percentage of nucleotides in contigs which aligned to the reference assembly was recorded.

Figure 3 shows that both the probabilistic de Bruijn graph and our structure have the same optimal Bloom filter size (11 bits per k-mer, total structure size of 13.82 bits and 13.62 per k-mer respectively). In the case of the probabilistic de Bruijn graph, the marking structure is prominent. This is because the graph has a significant amount of complex k-mers, most of them are linked to false positive nodes. For the graph equipped with the cF P structure, the marking structure only records the actual complex nodes; it occupies consistently 0.49 bits per k-mer. Both structures have comparable memory usage.

However, Figure 3 shows that the probabilistic de Bruijn graph produces assemblies which strongly depend on the Bloom filter size. Even for large sizes, the probabilistic graph assemblies differ by more than 3 Kbp to the reference assembly. We observed that the majority of these differences were due to missing regions in the probabilistic graph assemblies. This is likely caused by extra branching, which shortens the lengths of some contigs (contigs shorter than 100 bp are discarded).

Below ≈ 9 bits per k-mer, probabilistic graph assemblies significantly deteriorate. This is consistent with another article [START_REF] Pell | Scaling metagenome sequence assembly with probabilistic de bruijn graphs[END_REF], which observed that when

G G G G G G G G G G G G G G G 3 100 4527
Bloom filter size (bits/kmer)

G G G G G G G G G G G G G G G 0
Bloom filter size (bits/kmer) Differences with exact assembly (Kbp)

Fig. 3: Whole structures size (Bloom filter, marking structure, and cF P if applicable) of the probabilistic de Bruijn graph with (top right) and without the cF P structure (top left), for an actual dataset (E. coli, k = 23). All plots are in function of the number of bits per k-mer allocated to the Bloom filter. Additionally, the difference is shown (bottom left and bottom right) between a reference assembly made using an exact de Bruijn graph, and an assembly made with each structure.

the false positive rate is over 18% (i.e., the Bloom filter occupies ≤ 4 bits per k-mer), distant nodes in the original graph become connected in the probabilistic de Bruijn graph. To sum up, assemblies produced by the probabilistic de Bruijn graph are prone to randomness, while those produced by our structure are exact.

de novo assembly

We assembled a complete human genome (NA18507, SRA:SRX016231, 142.3 Gbp of unfiltered reads of length ≈ 100 bp, representing 47x coverage) using Minia. After k-mer counting, 2,712,827,800 solid k-mers (d = 3) were inserted in a Bloom filter dimensioned to 11.1 bits per solid k-mer. The cF P structure contained 78,762,871 k-mers, which were stored as a sorted list of 64 bits integers, representing 1.86 bits per solid k-mer. A total of 166,649,498 complex k-mers (6% of the solid k-mers) were stored in the marking structure using 4.42 bits per solid k-mer (implementation uses 8⌈ k 32 ⌉ bytes per k-mer). Table 1 shows the time and memory usage required for each step in Minia.

We compared our results with assemblies reported by the authors of ABySS [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF], SOAPdenovo [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF], and the prototype assembler from Conway and Bromage [START_REF] Conway | Succinct data structures for assembling large genomes[END_REF]. Table 2 shows the results for four classical assembly quality metrics, and the time and peak memory usage of the compared programs. We note that Minia has the lowest memory usage (5.7 GB), seconded by the assembler from Conway and Bromage (32 GB). The wall-clock execution time of Minia (23 h) is comparable to the other assemblers; note that it is the only single-threaded assembler.

The N50 metric of our assembly (1.2 Kbp) is slightly above that of the other assemblies (seconded by SOAPdenovo, 0.9 Kbp). All the programs except one assembled 2.1 Gbp of sequences.

We furthermore assessed the accuracy of our assembly by aligning the contigs produced by Minia to the GRCh37 human reference using GASSST [START_REF] Rizk | GASSST: global alignment short sequence search tool[END_REF]. Out of the 2,090,828,207 nucleotides assembled, 1,978,520,767 nucleotides (94.6%) were contained in contigs having a full-length alignment to the reference, with at least 98% sequence identity. For comparison, 94.2% of the contigs assembled by ABySS aligned full-length to the reference with 95% identity [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF].

To test another recent assembler, SparseAssembler [START_REF] Ye | Exploiting sparseness in de novo genome assembly[END_REF], the authors assembled another dataset (NA12878), using much larger effective k values. SparseAssembler stores an approximation of the de Bruijn graph, which can be compared to a classical graph for k ′ = k + g, where g is the sparseness factor. The reported assembly of the NA12878 individual by SparseAssembler (k + g = 56) has a N50 value of 2.1 Kbp and was assembled using 26 GB of memory, in a day. As an attempt to perform a fair comparison, we increased the value of k from 27 to 51 for the assembly done in Table 2 (k = 56 showed worse contiguity). The N50 obtained by Minia (2.0 Kbp) was computed with respect to the size of SparseAssembler assembly. Minia assembled this dataset using 6.1 GB of memory in 27 h, a 4.2× memory improvement compared to SparseAssembler.

Step [START_REF] Conway | Succinct data structures for assembling large genomes[END_REF], ABySS [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF], and SOAPdenovo [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF]. Contigs shorter than 100 bp were discarded. Assemblies were made without any pairing information.

To the best of our knowledge, Minia is the first method that can create contigs for a complete human genome on a desktop computer. Our method improves the memory usage of de Bruijn graphs by two orders of magnitude compared to ABySS and SOAPdenovo, and by roughly one order of magnitude compared to succinct and sparse de Bruijn graph constructions. Furthermore, the current implementation completes the assembly in 1 day using a single thread.

De Bruijn graphs have more NGS applications than just de novo assembly. We plan to port our structure to replace the more expensive graph representations in two pipelines for reference-free alternative splicing detection, and SNP detection [START_REF] Peterlongo | Identifying SNPs without a reference genome by comparing raw reads[END_REF][START_REF] Sacomoto | KISSPLICE: de-novo calling alternative splicing events from RNA-seq data[END_REF]. We wish to highlight three directions for improvement. First, some steps of Minia could be implemented in parallel, e.g. graph traversal. Second, a more succinct structure can be used to mark complex k-mers. Two candidates are Bloomier filters [START_REF] Chazelle | The bloomier filter: an efficient data structure for static support lookup tables[END_REF] and minimal perfect hashing.

Third, the set of critical false positives could be reduced, by exploiting the nature of the traversal algorithm used in Minia. The traversal ignores short tips, and in general, graph regions that are eventually unconnected. One could then define n-th order critical false positives (n-cF P ) as follows. An extension of a true positive graph node is a n-cF P if and only if a breadth-first search from the true positive node, in the direction of the extension, has at least one node of depth n + 1. In other words, false positive neighbors of the original graph which are part of tips, and generally local dead-end graph structures, will not be flagged as critical false positives. This is an extension of the method presented in this article which, in this notation, only detects 0-th order critical false positives.

Fig. 1 :

 1 Fig. 1: A complete example of removing false positives in the probabilistic de Bruijn graph. (a) shows S, an example de Bruijn graph (the 7 non-dashed nodes), and B, its probabilistic representation from a Bloom filter (taking the union of all nodes). Dashed rectangular nodes (in red in the electronic version) are immediate neighbors of S in B. These nodes are the critical false positives. Dashed circular nodes (in green) are all the other nodes of B; (b) shows a sample of the hash values associates to the nodes of S (a toy hash function is used); (c) shows the complete Bloom filter associated to S; incidentally, the nodes of B are exactly those to which the Bloom filter answers positively; (d) describes the lower bound for exactly encoding the nodes of S (selfinformation) and the space required to encode our structure (Bloom filter, 10 bits, and 3 critical false positives, 6 bits per 3-mer).

Fig. 2 :

 2 Fig. 2: (a) Structure size (Bloom filter, critical false positives) in function of the number of bits per k-mer allocated to the Bloom filter (also called ratio r) for k = 32. The trade-off that optimizes the total size is shown in dashed lines. (b) Optimal size of the structure for different values of k.

Table 1 :

 1 Details of steps implemented in Minia, with wall-clock time and memory usage for the human genome assembly. For constant-memory steps, memory usage was automatically set to an estimation of the final memory size. In all steps, only one CPU core was used.This article introduces a new, space-efficient representation of the de Bruijn graph. The graph is implicitly encoded as a Bloom filter. A subset of false positives, those which introduce false branching from true positive nodes, are recorded in a separate structure. A new marking structure is introduced, in order for any traversal algorithm to mark which nodes have already been visited. The marking structure is also space-efficient, as it only stores information for a subset of k-mers. Combining the Bloom filter, the critical false positives structure and the marking structure, we implemented a new memory-efficient method for de novo assembly (Minia).

		Time (h) Memory (Gb)
	k-mer counting	11.1	Constant (set to 4.0)
	Enumerating positive extensions 2.8	3.6 (Bloom filter)
	Constructing cF P	2.9	Constant (set to 4.0)
	Assembly	6.4	5.7 (Bloom f.+ cF P + mark. struct.)
	Overall	23.2	5.7
	7 Discussion

Table 2 :

 2 de novo human genome (NA18507) assemblies reported by our assembler (Minia), Conway and Bromage assembler

Source code available at http://minia.genouest.org/

Acknowledgments

The authors are grateful to Dominique Lavenier for helpful discussions and advice, and Aurélien Rizk for proof-reading the manuscript. This work benefited from the ANR grant associated with the MAPPI project (2010-2014).