Shin-Ichiro Sato
email: sato.shinichiro2@qst.go.jp

Takahiro Makino

Takeshi Ohshima

Tomihiro Kamiya

Wataru Kada

Osamu Hanaizumi

Veljko Grilj

Natko Skukan

Ivan Sudić

Milko Jakšić

Michal Pomorski

Gyorgy Vizkelethy

Transient Current Induced in Thin Film Diamonds by Swift Heavy Ions

Keywords: Diamond crystal, radiation induced effects, detectors, membranes

Single crystal diamond is a suitable material for the next generation particle detectors because of the superior electrical properties and the high radiation tolerance. In order to clarify fundamental characteristics of diamond particle detectors, transient currents generated in diamonds by single swift heavy ions (26 MeV O 5+ and 45 MeV Si 7+ ) were investigated. Two dimensional mappings of transient currents by single ion hits were also measured. In the case of 50 µm-thick diamond, both the signal height and the collected charge were reduced by the subsequent of ion hits and the charge collection time was extended. Strong signals appeared occasionally when most of the signals were fully degraded. These results were thought to be attributable to the polarization effect in diamond and it appeared only at the positive bias. In the case of 6 µm-thick diamond membrane, an "island" structure was found in the 2D mapping of transient currents.

Signals in the islands

showed different applied bias dependence from signals in other regions, indicating different crystal quality. The instability of signals resulting from the polarization effect was also observed in the 6 µm-thick diamond membrane. Simulation study of transient currents based on the Shockley-Ramo theorem clarifies that the polarization appears only in the irradiated region and both the electric field and the hole lifetime are reduced by the polarization effect.

Introduction

Diamond is an attractive material for radiation detectors to be used in harsh environments because of its superior electronic properties, such as high carrier mobility, wide bandgap, high radiation tolerance, and high breakdown voltage. Although the study for diamond radiation detectors started in the 1940s [1,2], many problems about the sample preparation, such as availability, crystal quality, and control of impurity, made the development difficult because natural diamonds were used in the earliest years. Since significant improvement has been recently achieved in the crystal growth technology of diamond, research and development for single-crystal diamond radiation detectors has been extensively carried out in recent years and their superior characteristics has been demonstrated [3]. For instance, it has been recently reported that single-crystal diamonds grown by chemical vapor deposition (CVD) showed almost 100 % charge collection efficiency and high energy resolution (less than 1 %) in the detection of 5.486 MeV alpha particles from 241 Am [4].

With the development of research on diamond particle detectors, however, some issues has also come out. One of the critical issues on diamond particle detectors is the polarization effect: temporal degradation of signal amplitude during ion irradiation. Being generally known to occur in insulator crystals [5], the polarization effect has been one of the important topics of interest in diamond detectors [6][7][8]. The polarization effect disappears by changing the polarity of applied bias and occasionally mitigated during ion irradiation [9]. This is more pronounced for short range ions, which produce high ionization in a small volume [10]. When charge carriers generated by radiation are trapped in deep levels in diamond, they cannot be released in a short time and give rise to a space charge (internal electric field). The external electric field is screened by the space charge and as a result, the collected charge signal amplitude decreases. This is a brief explanation of the polarization effect, and the detail mechanism is still less well understood. The origin of these deep levels is completely unidentified. It is also still unclear how the charge collection process of generated carriers are affected by the polarization effect.

In this study, transient currents in diamond generated by swift heavy ions were measured using Time Resolved Ion-Beam Induced Current (TRIBIC) technique [11] in order to investigate the charge collection process in diamonds. The TRIBIC technique is a powerful tool for investigating the collection process of generated charge by incident ions.

Two dimensional distribution of transient currents was also investigated.

Experimental

Two samples were used in this study: a single-crystal CVD (sc-CVD) diamond with the thickness of 50 µm and a thin sc-CVD diamond membrane with the thickness of 6 µm.

The fabrication method of diamond membrane was reported elsewhere [12]. The detectors were of electrical grade (substitutional nitrogen content [N] < 5 ppb) and optical grade ([N] < 5 ppm), respectively, which were commercially available from Element Six Ltd. Metal electrodes were formed on both sides of the samples to collect charges generated by incident ions (26 MeV O 5+ and 45 MeV Si 7+ ). Current-voltage (IV) and capacitance-voltage (CV) characteristics were investigated prior to TRIBIC measurements.

The experimental setup of TRIBIC measurement is illustrated in Fig. 1. Comparing to conventional Ion Beam Induced Charge (IBIC) technique [13], the TRIBIC measurement gives information on time variation of charge collection, i.e. transient current, in addition to generated charge by a single ion incidence. The detail of TRIBIC measurement system was reported elsewhere [11]. The samples were irradiated at normal incidence with a focused microbeam of 26 MeV O 5+ and 45 MeV Si 7+ ions at the Sandia National Laboratories (SNL), USA. Ion-induced transient currents generated by single ion hits were measured from the top electrode of sample by the 20 GHz oscilloscope (Tektronix, DPO 72004) via a linear amplifier (Picosecond Pulse Labs, 5840B) and a bias tee. The transient currents and their positional information were measured, and the two dimensional (2D) mapping of the transient currents was recorded. The spot size of the microbeam was about 1 µm in diameter. The microbeam was scanned from left to right and from bottom to top, and the scan area was 50×50 µm 2 . A bias voltage ranging from -50 V to +50 V was applied to the top electrode of sample (the same side as ion incidence) and the bottom electrode was grounded. The variation of transient currents due to applied bias voltages was also investigated.

Energy depositions and creation of number of electron-hole pairs by incident ions were calculated by SRIM [14]. The calculation results are shown in Fig. 2. The projected ranges of 26 MeV-O and 45 MeV-Si in diamond were calculated to be 7.7 µm and 8.7 µm, respectively. Therefore, 45 MeV-Si ions completely pass through the 6 µm-thick diamond membrane and 38.9 MeV (86.4%) of the energy was estimated to be deposited by electronic excitation. The diamond density of 3.515 g/cm 3 and the average ionization (electron-hole pair creation) energy of 13.2 eV [15] were used in these calculations. and charge collection efficiency (CCE) was also obtained from the collected charge, the deposited ion energy, and the average ionization energy. The charge generated in 50 µmthick diamond by 26 MeV Oxygen was estimated to be 0.318 pC. Figure 5 shows variations of the peak heights and CCEs of TRIBIC signals as a function of applied electric field. Both the peak height and the CCE increased with increasing absolute value of electric field. The value of CCE was higher at positive bias than at negative bias; 86.2 % at +10 kV/cm and 60.4 % at -10 kV/cm. This trend corresponds to the results reported elsewhere [4,10]. According to the Shockley-Ramo theorem [16], in the case of charge generation near one electrode, CCE as a function of applied electric field can be explained by a single-carrier Hecht equation:

Results and Discussion

TRIBIC signals of 50 µm-thick diamond

CCE = $% & ' 1 -exp - ' $% & (1).
Where W, E, and µτ are the thickness, the electric field, and the mobility-lifetime product, respectively. A uniform electric field in the diamond is assumed in Eq. ( 1). The fitting result is shown as the red line in Fig. 5. The values of µt product were (1.8±0.2)×10 -6 cm 2 /V for holes and (6.4±0.8)×10 -7 cm 2 /V for electrons. These values are substantially lower than the values reported elsewhere [4,10,17,18] which can be explained as follows.

Firstly, since density of electron-hole pairs generated by swift heavy ions is much higher than that of light ions like protons and alpha particles, some carriers recombine promptly after generation, the value of CCE for swift heavy ions are generally lower than that for light ions [19]. As a result, the lower value of µτ product was obtained from the CCE for heavy ions than the reported values. Also, since we performed irradiation experiments using this sample prior to this study, significant radiation degradation might have been already caused even though diamond is super radiation tolerant. The applied electric fields were +10 kV/cm (black) and -10 kV/cm (red). The linear amplifier was not used for this measurement. The microbeam was scanned from left to right, and from bottom to top, and the transient current was recorded only when the signal height was higher than a trigger level of the oscilloscope. Therefore, the signal amplitude was reduced gradually from the bottom left to the top right in the 2D maps and eventually the signal height fell below the trigger level.

As a result, no signal was observed even if the ions hit, and the blank region appeared at the upper region of 2D maps. Since a saturation tendency of decrease in TRIBIC signal which was observed after ~1000 hits was just due to the trigger level, the actual signal amplitude would continue to decrease after ~1000 hits. However, interestingly, strong These results are obviously attributable to the polarization effect. When ion induced carriers (electrons and holes) are trapped in deep levels in diamond, the trapped carriers interfere with the external electric field and the charge collection is interrupted. The detail of polarization effect is discussed in section 3.3. We therefore conclude that the inhomogeneous 2D maps in Figs. 6(a) and (b) do not indicate the inhomogeneity of crystal quality, but are simply caused by enhancing the polarization effect with increasing number of ion hits. The results of Fig. 6 also indicate that the polarization effect induced by previous ion hits are widely spread bidimensionally and the degradation of signal amplitude occurs even if the position of subsequent ion hits is different from that of precedent ion hits. The spatial distribution of the polarization effect could be estimated to be over the range of microns, although it is difficult to estimate the exact value in this study.

The strong signal F occasionally observed after 1200 hits were comparable to the signal before the polarization occurred (the signal A), and the weak signals (the signal E) appeared for a while after a strong signal appeared once. It is speculated that the polarization effect was temporarily cleared out by excess accumulation of carriers due to the repetition of ion hits [9]. On the other hand, the polarization effect is polaritydependent and the inhomogeneous mapping was never observed when the negative bias was applied. Also, the signal height at -10 kV/cm, which is shown as the signal A' in Fig. 6(d), was about 60% compared to the signal at +10 kV/cm.

TRIBIC signals of 6µm-thick diamond membrane

The TRIBIC signals of the 6 µm-thick diamond membrane by 45 MeV Si 7+ ion hits are summarized in Fig. 7. Figures 7(a At the negative bias condition (Fig. 7(c)), the signals at the islands (the signal G) were stronger than the signals at the other region (the signal H). Meanwhile, at the positive bias condition (Fig. 7(d)), the signals at the other region (the signal K) were stronger than that at the islands (the signal I). In addition, the weaker signals were found at the coast of islands (the signal J). The signals at the other region appeared sparsely, indicating most of signals generated by ion hits at the other region were below the trigger level (~0.2 mA) and could not be detected by the oscilloscope. The signal height at the other region at -5 kV/cm (the signal H) was about 60% compared to +5 kV/cm (the signal K).

The trend observed at the other region is very similar to the result of Fig. 6, indicating the polarization effect could also appear in the 6 µm-thick diamond membrane. Although the polarization effect has been often observed in diamond whose thickness is considerably longer than projected range of ions, swift heavy ions which have high electronic energy deposition can induce the polarization effect even if diamond is thin enough for ions to penetrate. On the other hand, the signals at the islands showed a different trend from the other region, suggesting the existence of dislocation, inclusion, or grain boundary. The difference of bias polarity dependence is thought to be due to the difference of crystal quality in the islands. It is unclear at the present stage whether the island structure was formed during crystal growth or during the RIE etching for the fabrication of membrane structure.

Analysis of the polarization effect based on the Shockley-Ramo theorem

In order to analyze the signal degradation due to the polarization effect shown in Fig. 6(d), a simulation of transient currents based on the Shockley-Ramo theorem developed by Li et al. [20] was performed. A one-dimensional model of charge collection in semiconductors can be obtained from the Shockley-Ramo theorem and this is suitable to clarify how the polarization effect affects the charge collection process. Figure 8(a) shows the conceptual diagram of the model proposed in this study. It was assumed that an electric field was formed uniformly in diamond and the polarization occurred only in the irradiated region. In this case, the transient current, i(t), is represented by the following equations:

𝑖 𝑡 = 𝑖 / 𝑡 + 𝑖 1 𝑡 (2)
𝑖 / 𝑡 = 2 3' exp - 5 % 6
Λ 𝑥, 𝑡 • 𝑣 / 𝑥 ' < 𝑑𝑥

𝑖 1 𝑡 = 2 3' Λ 𝑥, 𝑡 • 𝑣 1 𝑥 • exp - 5 % >? @ < 𝑑𝑥 + Λ 𝑥, 𝑡 • 𝑣 1 𝑥 • exp - 5 % >A ' @ 𝑑𝑥 (3)
Where i e (t), i h (t), and t denote the electron current, the hole current, and the time, respectively. q, Λ, G, R, and W are the electron charge, the electronic energy deposition by incident ions, the average ionization energy (13.2 eV [15]), projected range of incident ions, and the sample thickness. The data of Λ(x, t = 0) was obtained from Fig. 2. The electron lifetime, τ e was assumed to be constant for the polarization effect. On the other hand, it was assumed that the hole lifetime in the irradiated region (0 ≤ 𝑥 ≤ 𝑅), τ h1 was reduced by the polarization effect and the hole lifetime in the other region (R< 𝑥 ≤ 𝑊), τ h2 was unchanged. The drift velocity, v is represented by the Canali model [21]:

𝑣 H = $ I & J KL M I N O P QI 𝑖 = 𝑒, ℎ (5).
Where µ, v s , and E are the drift mobility, the saturation velocity, and the electric field. The subscripts, e and h denote electron and hole, respectively Constant values were used for both the electron and the hole mobilities, and the drift velocity was changed only by the change in electric field. Values of 1.5×10 7 cm/s and 1.05×10 7 cm/s were used for electron and hole saturation velocity [2]. Finally, the calculated current was convolved with a

Gaussian that included all terms of the system resolution, such as stray capacitance, etc.

The typical response time was 6.5×10 -10 s for the 50 µm-thick diamond and 3.0×10 -10 s for the 6 µm-thick diamond membrane.

Initial values of mobility and lifetime were determined from fittings of the signals A and A' in Fig. 6(d), which were the transient currents without polarization. The obtained values are summarized in Table 1 and the fitting results are shown in Fig. 8(b). The products of mobility and lifetime for electrons and holes in Table 1 were slightly smaller than the µτ products which were obtained in Fig. 5. This is attributed to the difference of incident ions [19]. This model also explains well the transient currents in the 6-µm thick diamond membrane.

45 MeV-Si ions completely penetrate the 6-µm diamond membrane and induce the polarization, i.e. the decrease in both electric field and hole lifetime, throughout the depth direction. In this case, however, the electron current is comparable with the hole current and we were unable to determine individual parameters for electrons and holes. Therefore the same parameters were used for electrons and holes to simulate the signal K. The results are shown as blue lines in Fig. 8(c). The calculated transient current was in agreement with the signal J when the electric field and hole lifetime were reduced from +5.0 kV/cm and 5.0 ns to +2.0 kV/cm and 3.0 ns, respectively. Further discussion should be conducted for the change in electron lifetime in the irradiated region, although no variation of electron lifetime was assumed in this study. Since the polarization effect is observed at the condition of positive bias only, it is difficult to discuss the existence of electron lifetime killer at the present stage.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Experimental setup of TRIBIC measurement.

Figure 3

 3 Figure3shows the IV characteristics of the 50 µm-thick diamond. Although the IV

Figure 4

 4 Figure 4 shows typical transient currents (TRIBIC signals) at the biases of ±50 V (±10

Fig. 3 .Fig. 4 .

 34 Fig. 3. Current-voltage characteristics of the 50 µm-thick diamond. Negative values in current were converted to the positive values in order to show in the semi-logarithmic scale.

Fig. 5 .

 5 Fig. 5. Va r i a t i o n s o f p e a k height of TRIBIC signal (closed squares) and charge collection efficiency (CCE, red open circuits) of the 50 µm-thick diamond due to 26 MeV O 5+ ion hits as a function of electric field. The line in red are the fitting curve of CCE by Eq. (1).

Fig. 6 .

 6 Fig. 6. (a) Two dimensional maps of peak height of TRIBIC signals and (b) collected charge induced by 45 MeV Si 7+ ion hits. The applied bias was +50 V (+10 kV/cm) and the irradiation area was 50×50

  signals appeared sparsely at the upper region of 2D mapping (after ~1200 hits, see the shaded zone in (c)). The typical weak and strong TRIBIC signals were shown as the signals E and F in (d). The shape of signals E and F was unchanged from signals D and A, respectively.

  ) and (b) show the 2D maps of peak height of TRIBIC signals at the biases of -3 V (-5 kV/cm) and +3 V (+5 kV/cm), respectively. The 2D maps of collected charges are not shown here since they were qualitatively no different from that of peak heights. The ion irradiation at a strong positive bias (+50 kV/cm) was performed just before the TRIBIC measurements and "island" structure appeared in the 2D map. Examples of TRIBIC signals at ±5 kV/cm are shown in (c) and (d), respectively. These signals are identified by the letters of G to K and their locations are listed in the table in Fig. 7. Figures 7(d) and (e) show the histograms of peak height and collected charge, respectively.

Fig. 7 .

 7 Fig. 7. Summary of TRIBIC data in the 6 µm-thick diamond membrane due to 45 MeV Si 7+ ion hits at ±5 kV/cm. (a) and (b) show the 2D maps of the peak heights at -5 kV/cm and +5 kV/cm, (c) and (d) show typical transient currents at -5 kV/cm and +5 kV/cm, respectively. (e) and (f) are histograms of the peak height and the collected charge, respectively.

  Although the polarization effect simply induces the reduction of electric field in diamond according to the general understanding, it was shown from the simulation study that the decrease in hole lifetime should also be taken into account to fit degraded transient currents after the polarization. Assuming both the electric field and the hole lifetime are reduced by the polarization effect only in the irradiated region (0 ≤ 𝑥 ≤ 𝑅), transient currents were calculated and the results for the 50-µm thick diamond are illustrated in Fig. 8(b). In Fig. 8(b), the experimental transient currents as well as the calculated results are shown. The signal D, which is the transient current after the polarization, could be simulated by considering both the decrease in hole lifetime and electric field in the irradiated region. The results showed the hole lifetime and the electric field were reduced from 5.45 ns and +10 kV/cm and to 1.1 ns and +0.50 kV/cm, respectively.

Fig. 8 .

 8 Fig. 8. (a) One-dimensional model of the polarization effect based on the Shockley-Ramo theorem. (b) and (c) are calculation results of the 50 µm-thick diamond (red lines) and the 6 µm-thick diamond membrane (blue lines), respectively. Electron and hole currents are drawn as the dashed and dasheddotted lines, respectively.

  Transient currents (TRIBIC signals) of diamonds induced by single swift heavy ion hits and their 2D mappings were investigated in this study. Variations of TRIBIC signals and CCEs due to different diamond thicknesses, applied biases, and bias polarities were clarified. Change in the TRIBIC signals by subsequent of ion hits was also investigated.From the results of TRIBIC measurements in the 50 µm-thick diamond by irradiated by 45 MeV-Si ions, the signal degradation and instability originating from the polarization effect were observed. At positive bias, signals were drastically degraded by the subsequent ion hits and became unstable with increasing number of ion hits. An increase in the charge collection time in addition to the decrease in signal height were observed.After 1200 ion hits, non-degraded signals appeared occasionally with fully degraded signals. The spatial distribution of the polarization effect could be estimated to be in the range of microns. At negative bias, on the other hand, neither signal degradation nor instability of TRIBIC signals appeared with the subsequent of ion hits. The signal height at negative bias was about 60% compared to the signal height at positive bias.The similar trend was found in the TRIBIC results of the 6 µm-thick diamond membrane by 45 MeV-Si ion hits. This fact strongly indicates that the polarization effect could occur in diamond membranes by swift heavy ions. Island structure was found in the 2D maps and signals in the islands showed different trend from the other regions. The inhomogeneity is thought to be caused by the presence of dislocation, inclusion or grain boundary in diamond.

Table 1 .

 1 Parameters used for the simulation of transient currents in Fig. 6(d).

	Electron mobility, µ e (cm 2 /Vs)	2.2×10 2
	Hole mobility, µ h (cm 2 /Vs)	3.1×10 2

Acknowledgement

This work was supported by JSPS KAKENHI Grant Number 26249149. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

* Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia

Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

Also, simulation study based on the Shockley-Ramo theorem was performed to understand the observed polarization effect. The proposed model explains well transient currents in both thick and thin diamond. It was shown from the simulation study that the polarization occurred only in the irradiated region and induced significant reduction of hole lifetime in addition to electric field.


Sawsan Kanj

Bru Thomas

¨ls

Ste ´phane Gazut

Shared Nearest Neighbor Clustering in a Locality Sensitive Hashing Framework

Keywords: density-based methods, metagenomic data, sequence clustering, locality sensitive hashing

We present a new algorithm to cluster high-dimensional sequence data and its application to the field of metagenomics, which aims at reconstructing individual genomes from a mixture of genomes sampled from an environmental site, without any prior knowledge of reference data (genomes) or the shape of clusters. Such problems typically cannot be solved directly with classical approaches seeking to estimate the density of clusters, for example, using the shared nearest neighbors (SNN) rule, due to the prohibitive size of contemporary sequence datasets. We explore here a new approach based on combining the SNN rule with the concept of locality sensitive hashing (LSH). The proposed method, called LSH-SNN, works by randomly splitting the input data into smaller-sized subsets (buckets) and employing the SNN rule on each of these buckets. Links can be created among neighbors sharing a sufficient number of elements, hence allowing clusters to be grown from linked elements. LSH-SNN can scale up to larger datasets consisting of millions of sequences, while achieving high accuracy across a variety of sample sizes and complexities.

INTRODUCTION

C lustering is usually defined as the task of unsupervised learning, where the class labels of the it has been shown in [START_REF] Wu | The uniform effect of k-means clustering[END_REF] that K-means tends to produce clusters of relatively uniform sizes and globular shapes, even if the data structure is endowed with varying cluster sizes or different shapes. This bias is known as the uniform effect of the K-means. Moreover, the number of clusters K has to be specified a priori, which is not trivial when no prior knowledge is available. To address these problems, methods based on estimating the density and/or the similarity among instances have been introduced [START_REF] Jarvis | Clustering using a similarity measure based on shared near neighbors[END_REF][START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF].

In [START_REF] Erto ¨z | Finding Topics in Collections of Documents: A Shared Nearest Neighbor Approach[END_REF], the authors presented an effective clustering method based on two key notions: the similarity between neighboring elements and the density around instances. This method, shared nearest neighbors (SNN), is a density-based clustering method and incorporates a suitable similarity measure to cluster data. After finding the nearest neighbors of each element and computing the similarity between pairs of points, SNN identifies core points, eliminates noisy elements, and builds clusters around the core elements. This method can yield improved performance compared with other clustering approaches with data of varying densities, and it can automatically infer the number of output clusters. However, this method has complexity O(n 2 ), where n is the number of instances in the dataset, arising from the computation of the similarity matrix, which can be prohibitive when dealing with very large amounts of data.

One interesting concept to reduce the burden of computing the similarity matrix is locality sensitive hashing (LSH). This concept was initially introduced to find approximate near neighbor information in high-dimensional space [START_REF] Indyk | Approximate nearest neighbors: towards removing the curse of dimensionality[END_REF][START_REF] Har-Peled | Approximate nearest neighbor: Towards removing the curse of dimensionality[END_REF]Wang et al., 2014a). The key idea is to hash elements into different buckets; then for a query instance x, to use instances stored in buckets containing x as candidates for near neighbors. This approximation reduces the query time complexity to O( log n) instead of O(n) (O(n) is the complexity for searching nearest neighbors for one instance). Therefore, the similarity matrix computation time can be reduced to O(n log n).

We propose here to retain the basic principle of LSH by randomly splitting the dataset into a number of smaller-sized subsets, using a family of hashing functions, so that similar elements will be hashed together with high probability. We then look for nearest neighbors of each element in its bucket, and construct links among elements sharing a significant number of neighbors to generate clusters. The proposed method, called LSH-SNN, has the advantage of reducing the complexity for computing the similarity matrix, while maintaining the same level of clustering accuracy.

In this study, we have evaluated the performance of the LSH-SNN method on metagenomic datasets of various sizes and complexities. We have also compared the results with another density-based clustering algorithm and the K-means method implemented in a popular sequence clustering software called Me-taCluster [START_REF] Yang | Metacluster: Unsupervised binning of environmental genomic fragments and taxonomic annotation[END_REF]Wang et al., 2014b). Many computational tools have been proposed in the literature to analyze metagenomic sequences that are generated from microbial communities. These tools can be grouped into two main categories: (1) supervised and (2) unsupervised methods. Supervised methods, often relying on sequence similarities and alignments of DNA fragments to reference sequences of known taxonomic origins, include tools such as Kraken [START_REF] Wood | Kraken: Ultrafast metagenomic sequence classification using exact alignments[END_REF], Clark [START_REF] Ounit | Clark: Fast and accurate classification of metagenomic and genomic sequences using discriminative k-mers[END_REF], and Diamond [START_REF] Buchfink | Fast and sensitive protein alignment using diamond[END_REF]. Unsupervised methods group metagenomic fragments based on intrinsic features, such as the statistics of l-mer frequencies extracted from fragments. Unsupervised methods, such as MetaCluster [START_REF] Yang | Metacluster: Unsupervised binning of environmental genomic fragments and taxonomic annotation[END_REF], MaxBin [START_REF] Wu | Maxbin: An automated binning method to recover individual genomes from metagenomes using an expectation-maximization algorithm[END_REF], and MetaProb [START_REF] Girotto | Metaprob: Accurate metagenomic reads binning based on probabilistic sequence signatures[END_REF], became attractive due to the lack of reference genomes for the bulk of micro-organisms identified by culture-independent approaches.

The rest of this article is organized as follows: Section 2 surveys related work on clustering; Section 3 recalls some background on local sensitive hashing and the SNN methods; and Section 4 introduces our method based on the combination of local sensitive hashing and SNN. Experimental results on synthetic and real datasets are illustrated in Sections 5 and 6, respectively, whereas Section 7 concludes the article.

RELATED WORK

Clustering methods look for similarities within a set of instances without any need for prior data labeling. Numerous methods have been proposed in literature to deal with clustering tasks. Existing algorithms can be grouped into five categories as proposed in [START_REF] Berkhin | A survey of clustering data mining techniques[END_REF] and [START_REF] Kotsiantis | Recent advances in clustering: A brief survey[END_REF]: partitioning methods, hierarchical methods, density-based methods, grid-based methods, and model-based methods. Hereafter, we will describe the main characteristics of these methods.

Partitioning methods construct K partitions of the data by grouping instances around the gravity center of each cluster. They can be divided into two main groups: the centroid methods such as K-means [START_REF] Hartigan | Algorithm as 136: A k-means clustering algorithm[END_REF], and the medoids ones [START_REF] Kaufman | Clustering by Means of Medoids[END_REF] such as the K-modes [START_REF] Huang | Extensions to the k-means algorithm for clustering large data sets with categorical values[END_REF] and the K-prototypes algorithms [START_REF] Huang | Clustering large data sets with mixed numeric and categorical values[END_REF]. Partitioning methods are simple to implement; however, the number of clusters K should be specified.

Hierarchical methods build a tree hierarchy, known as dendrogram, to form clusters in two different manners: agglomerative (bottom-up) and divisive (top-down). The former starts with singleton clusters and recursively merges them in a bottom-up strategy, whereas the latter breaks the dataset into smaller clusters in a top-down strategy. They use various local criteria to join or split clusters. Hierarchical methods have the advantage of handling any form of similarity without requiring the number of clusters to be known in advance. However, to construct a dendrogram, they suffer from time and space complexities, which are quadratic with respect to the number of clusters. Hierarchical clustering includes methods such as: BIRCH [START_REF] Zhang | Birch: An efficient data clustering method for very large databases[END_REF], CURE [START_REF] Guha | Cure: An efficient clustering algorithm for large databases[END_REF], and CHAMELEON [START_REF] Karypis | Chameleon: Hierarchical clustering using dynamic modeling[END_REF].

Density-based methods generate clusters based on the density of instances in a region. These methods are related to different concepts defining a point's nearest neighbors, such as density, connectivity, and boundary. Density-based methods can find arbitrary-shaped clusters; however, they output border instances, which may be unclustered and considered as outliers. Existing methods include DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF], OPTICS [START_REF] Ankerst | Optics: Ordering points to identify the clustering structure[END_REF], DENCLUE [START_REF] Hinneburg | An efficient approach to clustering in large multimedia databases with noise[END_REF], Jarvis-Patrick ( JP) [START_REF] Jarvis | Clustering using a similarity measure based on shared near neighbors[END_REF][START_REF] Jarvis | Clustering using a similarity measure based on shared near neighbors[END_REF]SNN (Erto ¨z et al., 2003) algorithms. SNN will be described in further detail in Section 3.2.

Grid-based methods quantize the space into a finite number of cells that form a grid structure. Clustering is then performed on the grid cells, instead of the database itself [START_REF] Liao | A grid-based clustering algorithm using adaptive mesh refinement[END_REF]. The main advantage of these methods is their fast processing time; however, they output clusters with either vertical or horizontal boundaries (no diagonal boundary can be detected). This category includes STING [START_REF] Wang | Sting: A statistical information grid approach to spatial data mining[END_REF], WaveCluster [START_REF] Sheikholeslami | Wavecluster: A multi-resolution clustering approach for very large spatial databases[END_REF], and CLIQUE [START_REF] Agrawal | Automatic subspace clustering of high dimensional data for data mining applications[END_REF] algorithms, among others.

In model-based clustering, it is assumed that the data are generated from K probability distributions and the goal is to find the distribution parameters (Yeung et al., 2001). Model-based methods are characterized by a small number of parameters; however, the computational burden can become significant if the number of distributions is large. Moreover, it is difficult to estimate the number of clusters. Many model-based clustering methods are described in the literature, such as Expectation-Maximization [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], Self Organized Maps net [START_REF] Kohonen | Self-Organizing Maps[END_REF], and AutoClass (Cheeseman and [START_REF] Cheeseman | Bayesian classification (autoclass): Theory and results[END_REF].

None of the categories mentioned earlier can directly cluster large amounts of instances of arbitrary shapes and at the same time automatically detect the appropriate number of clusters. The SNN algorithm from the density-based clustering category can deal with local density variations and automatically find clusters of different shapes. However, adapting this technique with massive data requires extensive storage and time costs, especially for the computation of the similarity matrix.

BACKGROUND

In this section, we briefly review some background on locality sensitive hashing (Subsection 3.1) and SNN algorithms (Subsection 3.2).

Locality sensitive hashing

LSH was first introduced in Indyk and Motwani (1998) as a classical geometric lemma on random projections, to quickly find similar items in large datasets. One or many families of hash functions map similar inputs to the same hash code. This hashing technique produces a splitting of the input space into many subspaces, called bins or buckets, with a high probability that instances originally close in their input space will be in the same bin or in adjacent bins within the LSH framework.

To alleviate the curse of dimensionality, each hash function projects the data to a lowerdimensional space (h : R d ! Z). Different techniques have been presented in the literature to generate hash functions. These techniques can be categorized into two families: min-hash [START_REF] Broder | On the resemblance and containment of documents[END_REF] and random projections (sim-hash). In document classification, min-hash is typically used when looking for textually similar documents by processing items and generating integers from strings of characters [START_REF] Leskovec | Mining of Massive Datasets[END_REF]. Random projections, on the other hand, are obtained via simple probability distributions such as p-stable distribution [START_REF] Datar | Locality-sensitive hashing scheme based on p-stable distributions[END_REF] and sign-random-projection [START_REF] Charikar | Similarity estimation techniques from rounding algorithms[END_REF].

Dealing with large datasets, LSH is typically combined with nearest neighbors techniques [START_REF] Buhler | Efficient large-scale sequence comparison by locality-sensitive hashing[END_REF] and used for data clustering [START_REF] Broder | On the resemblance and containment of documents[END_REF][START_REF] Boydell | High-throughput continuous clustering of message streams[END_REF]. To perform k-nearest neighbors, buckets [and sometimes their adjacent buckets [START_REF] Lv | Multi-probe lsh: Efficient indexing for high-dimensional similarity search[END_REF]] containing the query element are checked and all the existing instances are ranked according to their distances to the query element. To cluster high-dimensional data [START_REF] Haveliwala | Scalable techniques for clustering the web[END_REF][START_REF] Koga | Fast hierarchical clustering algorithm using locality-sensitive hashing[END_REF], similar elements contained in the same bucket can be joined to output clusters in a hierarchical way [START_REF] Rasheed | Efficient clustering of metagenomic sequences using locality sensitive hashing[END_REF].

Shared nearest neighbor

SNN is a density-based clustering approach that was successfully applied for finding groups of documents sharing a strong, coherent topic or theme [START_REF] Gowda | Agglomerative clustering using the concept of mutual nearest neighbourhood[END_REF][START_REF] Erto ¨z | Finding Topics in Collections of Documents: A Shared Nearest Neighbor Approach[END_REF][START_REF] Steinbach | The challenges of clustering high dimensional data[END_REF][START_REF] Berkhin | A survey of clustering data mining techniques[END_REF][START_REF] Moe ¨llic | Image clustering based on a shared nearest neighbors approach for tagged collections[END_REF][START_REF] Patidar | Analysis of different similarity measure functions and their impacts on shared nearest neighbor clustering approach[END_REF]. SNN handles clusters of widely differing sizes, densities, shapes, and in the presence of large amounts of noise and outliers. To exploit space density of the data, SNN relies on the concept of similarity based on the SNN approach. The similarity matrix is sparsified by keeping only the k-most nearest neighbors (knn). The SNN graph is then constructed by creating links between pairs of instances having each other in their respective knn lists. The weight of the link can be calculated either as the number of shared neighbors between two knn lists or using the ordering of these shared neighbors.

The algorithm determines the type of each instance (core, border, or noisy) by calculating its connectivity; that is, the number of links coming out of this instance, which will be compared with noisy and topic thresholds. Noisy instances are discarded and will never be used in the clustering process. Core instances form final clusters with their connected elements. The algorithm's behavior is controlled by four parameters: the number of nearest neighbors (noted as knn hereafter), the topic threshold, and two other thresholds governing the addition of elements to clusters. Depending on the parameter settings, many of the border instances can remain unlabeled because they are unconnected to core elements.

OUR PROPOSAL: THE LSH-SNN ALGORITHM

We describe the key idea of our algorithm, called LSH-SNN, in Section 4.1. We then describe how to tune the different parameters in Section 4.2.

Method description

SNN is a relatively effective unsupervised method to automatically find clusters of different shapes and densities. For n data items and knn nearest neighbors, the computational complexity of SNN is O(n 2 ), whereas its space complexity is O(knn à n). Hence, SNN can face important scalability issues when the number of instances n becomes large.

Our goal is to adapt a suitable framework for clustering large number of instances by using the SNN principle. It is motivated by large metagenomic datasets incurring high computational costs. To help reduce these costs, we consider the rationale of locality sensitive hashing as a framework for the development of the SNN method. In this framework, the n data instances are randomly partitioned into a number of smallersized subsets called buckets, and for each data point, we locate its approximate nearest neighbors inside its bucket. This approach, called LSH-SNN, has an advantage over SNN since it restricts the calculation of distances for each single item inside its bucket, whereas SNN needs to calculate n distance measures for each point before selecting the list of nearest neighbors.

LSH-SNN begins with the extraction of features from the sequence fragments by computing the frequencies of all possible l-mers (substrings of length l) in each of them (Section 4.1.1). Nearest neighbors of all sequences are then computed by applying the LSH technique, which involves splitting sequences into different buckets in such a way that similar sequences end up in the same bucket with a higher probability (Section 4.1.2). Elements stored in buckets containing a given sequence x are retrieved and ranked according to their distances to x to compute its nearest neighbors list. Shared neighbors are linked according to the SNN rule, and connected sequences form output clusters (Section 4.1.3). Finally, in the case of unclustered fragments, a last step is performed to assign them to the cluster that is the most similar in terms of l-mer distribution (Section 4.1.4).

4.1.1. l-mer frequency calculation. Sequence similarities are typically identified by comparing occurrence patterns of relatively short DNA substrings of length l between the sequences [START_REF] Yang | Metacluster: Unsupervised binning of environmental genomic fragments and taxonomic annotation[END_REF][START_REF] Tanaseichuk | Separating metagenomic short reads into genomes via clustering[END_REF].

Two broad scenarios can be used to assess l-mer-based similarities: Abundance-based methods make use of relatively large l values (l ! 20) to ensure the uniqueness of most l-mers [START_REF] Tanaseichuk | Separating metagenomic short reads into genomes via clustering[END_REF], whereas composition-based methods rely on smaller l values. Since DNA is a combination of four different types of nucleotides (A,T,G,C), there are at most 4 l l-mer combinations forming the feature vector. The frequency of each l-mer combination is normalized by dividing the number of occurrences by the fragment length.

Because of nucleotide base complementarities, the size of the feature vector can be reduced by half, that is, for a DNA sequence x of length s, the feature vector is given by:

x = x 1 ' x 2 ' . . . ' x m ð Þ
, where x i = f i s , f i is the frequency occurrence of a given l-mer combination and m represents the size of the vector (or number of descriptive attributes), m = 4 l =2 if l is odd, and (4 l + 4 l=2 )=2 if even.

4.1.2. Locality sensitive hashing. For convenience, we briefly recall some notations. Let X be the collection of n sequences of m-dimensional features, and let x 2 X denote an input sequence. Let k be the number of projections. For each i 2 [k], h i (x) is given by: h i (x) = sign(x:v i ), where v i is a vector whose components are randomly generated from a Gaussian distribution, for example, N (0' 1). This scalar projection gives one hash value for x. The hash code for x is then obtained by a concatenation of the k hash values, g(x) = (h 1 (x)' h 2 (x)' . . . ' h k (x)). LSH prepares r copies of g(:) to improve the hashing discriminative power [START_REF] Dasgupta | Fast locality-sensitive hashing[END_REF], Wang et al., 2014a) (to avoid confusion, k [in lower case] is the number of sampled bits, whereas K [in upper case] is the number of output clusters).

The feature vectors are first normalized with zero mean and unit variance. Each input sequence x is then indexed by a hash code g(x) = (h 1 (x)' h 2 (x)' . . . ' h k (x)) defining its bucket identity, and the hash code of all sequences in X constitutes a hash table. This projection produces a new k-dimensional space (k << m). Since the number of elements per bucket is typically much smaller than n, we need to ensure that similar sequences share the same bucket with a higher probability while minimizing random effects. To achieve this, r hash functions g 1 ' g 2 ' . . . ' g r are sampled independently, each generating a distinct hash table. For each sequence x, we then identify g 1 (x)' g 2 (x)' . . . ' g r (x) indexing the r buckets where x is mapped in each projection.

Note that any hash function could be applied in this step [START_REF] Pauleve | Locality sensitive hashing: A comparison of hash function types and querying mechanisms[END_REF] when distances are measured as angles between pairs of points. In this work, we demonstrate results based on random hash functions that are generated from a Gaussian distribution.

The projection of the data items into different buckets can be summarized as follows:

Algorithm 1: Hashing Input: Set of DNA sequences X of size n, number of repetitions r Output: Set of matrices of hash code T corresponding to X:

1: X = (X -mean(X))=std(X); 2: k = log (n);
8k is the number of projections or axis 3: for i = 1 to r do 4:

Create an m-by-k matrix A i where rows are identical and generated from a Gaussian distribution (N (0' 1)); 8m is the number of features 5:

T i = X:A i ; 6: end for For hashing the dataset, the time complexity is O(m • k • r) per sequence, since each sequence of dimension m will be processed by k hash functions repeated r times. Therefore, for n sequences, the time complexity of LSH is O(n • m • k • r). For a fixed l-mer (and hence of m) value, LSH has a complexity of

O(n • log (n) • r) (k = log (n); see Section 4.2).
4.1.3. Shared nearest neighbors. Once the space has been partitioned k • r times, a simple k-nearest neighbor classifier may be considered to find the nearest neighbors of a sequence x inside its bucket (i.e., sharing the same hash code) for all partitions. The union of r subsets of nearest neighbors for a given sequence is treated as its neighborhood list. Since the nearest neighbor lists are generated from sparsely populated buckets, the computational cost and runtime can be greatly improved.

Once the sets of nearest neighbors have been defined, the SNN algorithm follows two steps: computing the strength of links and the labeling of sequences. A link is created between two sequences x 1 and x 2 if they have each other in their respective neighborhood lists, and it can be weighted according to the product of positions of shared instances between these two lists, namely:

link(x 1 ' x 2 ) = X (knn + 1 -p 1 ) Ã (knn + 1 -p 2 )' (1)
where p 1 and p 2 are the positions of a shared neighbor in the lists of x 1 and x 2 . The knn lists are then transformed into a graph where sequences (nodes) are connected via weighted links. For each sequence x in the graph, the sum of its total links (conn x ) is computed to enable the selection of a subset of representative sequences according to a connectivity-based criterion (conn > topic threshold).

Algorithm 2 inset summarizes the SNN method. To check the nearest neighbors of a sequence x, n 0 distances need to be evaluated, where n 0 is the number of elements sharing the same hash code as x. Since we have n sequences, the time complexity for computing the nearest neighbor elements is

O(n • n 0 • m + C(knn))
, where C(knn) is a relatively small factor enabling the selection of knn near neighbors for each sequence [START_REF] Jarvis | Clustering using a similarity measure based on shared near neighbors[END_REF]. To assess the link between two sequences having each other in their respective knn list, two columns of size knn are selected and evaluated. The cost of this process amounts to O(n • knn • knn).

Therefore, the total complexity of the algorithm becomes:

Complexity = O(n • k • r + n • m • n 0 + n • knn • knn) = O(n • k • r + n • n 0 + n • knn • knn) (for fixed l -mer value) = O(n • log (n) • r + n • n 0 + n • n 0 ) knn = ffiffiffiffi n 0 p = O(n • log (n) • r + 2 • n • n 0 ) n 0 % n=2 log (n) %O(n • log (n) • r)
This algorithm is able to handle clusters of different densities, but it can leave a large number of nonnoisy sequences unclustered. To alleviate this problem, we define an additional step to relabel unclustered sequences.

4.1.4. Relabeling. This relabeling step was developed to reduce the number of unclustered sequences. Briefly, it identifies a subset of frequencies that are characteristic of each cluster and contributing most to the classifier's accuracy, discarding less relevant features. Each unclustered sequence is then added to the cluster that is the most closely related with respect to this subset of frequencies.

Relabeling proceeds by computing the mean of each cluster and dividing it by the mean of the other clusters. Discriminant l-mer frequencies, that is, those that most differentiate this cluster from others, are selected (cf line 7 in Algorithm 3 inset). Then, for a given sequence, we compute its distance to the mean of each cluster by using the subset of discriminant frequencies and assign it to the nearest cluster. If two clusters are nearly equally close to a given sequence, we keep the latter unlabeled to avoid increasing the number of misclassified instances.

Algorithm 2: SNN computation Input: Set of DNA sequences X of size n, set of matrices of hash code T corresponding to r projections of X, the number of nearest neighbors knn, strong, topic, and merge thresholds Output: Number of clusters K, predicted set of labels Ŷ corresponding to X: 1: for i = 1 to n do 2: for j = 1 to r do 3:

Select Z j the bucket indexing by the same hash code of x i , T ij ; 4:

Compute the distances of x i to each element in Z j ; 5:

Select the knn nearest neighbors elements and add them to N xi ; 6: end for 7: end for 8: for i = 1 to n do 9:

for j = 1 to r • knn do 10:

Select an element The implementation of this part of the algorithm is presented as a pseudo code in Algorithm 3 inset. Relabeling requires n • K operations, hence a complexity of O(n • K). The overall complexity of the algorithm depends on the complexity of the hashing functions and the SNN classifier. The total complexity is % O(n • log (n) • r).

Parameters in LSH-SNN

This section discusses the configuration of the LSH and SNN parameters, which impact the method's performance in terms of both runtime and clustering quality. Parameters were determined by grid search and focused on optimizing the V-measure (see Section 5.3).

LSH has two parameters, k and r: (1) The number k of bits sampled determines the number of instances inside the buckets, which, on average, is expected to be equal to n 0 = n 2 k . The total number of buckets is limited to max (n' 2 k ). If k takes a small value with respect to the number of sequences, then we would end up with a large number of sequences per bucket (n 0 ) and the time consumption of the SNN phase will be very high. On the other extreme, if k = n, we would get on average one sequence per bucket and there will be no knn lists to be constructed. In this study, we set k to log (n). ( 2) The number of projections r is the second parameter. For r = 1, we take the risk that two close elements end up in distinct buckets because of the random nature of the hashing. By increasing the number of projections, we increase the probability that these two elements are mapped to the same bucket. The number of projections r should, thus, be increased Algorithm 3: Relabeling of unclustered sequences Input: Set of DNA sequences X of size n, number of projections r, number of clusters K, predicted set of labels Ŷ corresponding to X Output: Predicted set of labels Ŷ:

1: for i = 1 to K do 2:
Compute the mean of cluster i (m i ); 3:

Compute the standard deviation of cluster i (sd i ); 4:

for j = 1 to K do 5:

Compute the mean of cluster j (m j ); 6:

for ind = 1 to m do 7:

then if (((m i [ind]=m j [ind] ! K) _ (m i [ind]=m j [ind] 1=K)) ^(sd i [ind] is small)) 8:
ind is a significant frequency; 9:

end if 10:

end for 11:

end for 12: end for 13: for i = 1 to n do 14:

if ŷi is not labeled then 15:

Compute the distance of x i to each cluster by using only the significant frequencies; 16:

Label x i with the cluster having the nearest distance to x i ; 17:

end if 18: end for to stabilize the results. On the other hand, for large values of r, distant elements may be mapped to the same bucket, provoking the r • knn lists to grow and ultimately leading to the same drawbacks as the initial SNN method. In this study, r was empirically set between 300 and 1200, depending on the size of the dataset.

Four parameters determine the outcome of the SNN procedure: These are the knn, topic, merge, and strong thresholds. (1) The size of the nearest neighbor list knn ultimately depends on the number of elements per bucket. To construct knn lists containing a sufficiently large number of closely related sequences that are consistent with the shared neighbors criterion, we need to choose a relatively large number of nearest neighbors. However, increasing knn may add more distant sequences to the same cluster, thus increasing the computational cost. On the other hand, decreasing knn will result in many smaller-sized clusters. A simple and pragmatic approach is to set knn = ffiffiffiffi n 0 p or knn = log (n 0 ), where n 0 = n 2 k . In our experiments, we fixed knn to ffiffi ( p n 0 ). (2) The topic threshold determines the proportion of highest connectivity links to be selected as representatives. This threshold ranged from 0:04 to 0:06 in our experiments. (3) The merge threshold represents the fraction of links to be used in the cluster merging process, and it was fixed to 0:02 in our experiments. (4) Finally, the strong threshold aims at reducing the number of unlabeled sequences (singleton clusters) and at picking representative elements. In our experiments, this parameter was set to 0:1.

EXPERIMENTS ON SEMI-SYNTHETIC DATASETS

In this section, we detail the experimental setup with the semi-synthetic datasets. We first describe the datasets (Section 5.1) and provide a short reminder about the methods that we compare our algorithm with (Section 5.2). We then detail the metrics used to evaluate the performance of our method (Section 5.3) and present and discuss the results in Section 5.4.

Datasets

We have used synthetic datasets of increasing sizes and complexities, composed of 600 base-pair length reads (DNA fragments) sampled from sequenced bacterial genomes that were assembled to create synthetic communities of known composition and abundance [START_REF] Gkanogiannis | A scalable assembly-free variable selection algorithm for biomarker discovery from metagenomes[END_REF], which were then used as ground truth to measure the clustering accuracy of the different algorithms evaluated. The mean coverage of the datasets was fixed to 1X, which means that, on average, a given position in the genome is covered by 1 read. The number of reads derived from each organism is equal to 5000.

To evaluate the performance of the various clustering modules on the benchmark datasets, we compared class memberships of elements (reads) in each dataset with the memberships induced by the clustering. Class memberships correspond to the genomes that the reads were sampled from, with the cardinality of the class set matching the samples' richness. For all the synthetic datasets, the read generation process was performed by using the Mason software [START_REF] Holtgrewe | Mason-A Read Simulator for Second Generation Sequencing Data[END_REF] with default error model parameters for Illumina reads (Mason can insert position-specific sequence modifications according to empirically calibrated and sequencing platform-dependent error models).

Benchmark methods

The LSH-SNN method is compared with MetaCluster [START_REF] Yang | Metacluster: Unsupervised binning of environmental genomic fragments and taxonomic annotation[END_REF], a popular compositional binning software based on the K-means algorithm with the Spearman footrule distance, which operates on relative rankings of the l-mer frequencies [START_REF] Yang | Metacluster: Unsupervised binning of environmental genomic fragments and taxonomic annotation[END_REF].

The main advantage of this approach is its simplicity, which underlies its ability to handle datasets harboring a relatively large number of genomes. However, its behavior can be sensitive to the random choice of initial cluster centers, and it may fail to output clusters when the data present nonglobular shapes. Moreover, the number of clusters should be predetermined, which is not trivial in the absence of prior knowledge [START_REF] Rokach | Clustering methods. Data Mining and Knowledge Discovery Handbook[END_REF].

The complexity of the K-means method is O(n

• m • K • t),
where n is the number of instances, m is the dimension of data, K is the number of clusters, and t is the number of iterations for convergence. As l-mer is fixed for all experiments, the overall complexity of MetaCluster becomes O(n • K • t).

We also compared our algorithm with the JP method [START_REF] Jarvis | Clustering using a similarity measure based on shared near neighbors[END_REF] coupled with the LSH indexing in a way similar to the LSH-SNN combination. JP also relies on the nearest neighbors similarity concept, but it simply merges elements in the same cluster if they share a sufficient number of neighbors, that is, the number of shared elements between two neighbors is greater than a user-predefined threshold kt, which we fixed to (r à knn

) 2 =2. The complexity of LSH-JP is O(n • k • r + n • n 0 + n • knn • knn) % O(n • log (n) • r).

Performance evaluation

We considered the following metrics to evaluate the performance of our method: Homogeneity, Completeness, V-measure, F-measure, and the Adjusted Rand Index.

Homogeneity evaluates the class distribution within each cluster; it is highest when each cluster contains only elements from a single class. Completeness, on the other hand, measures the distribution of cluster assignments within each class; it is highest when elements of a given class are assigned to a single cluster. Let C be the number of distinct organisms in a dataset of n sequences and K be the number of output clusters. These two measures are defined as follows [START_REF] Rosenberg | V-measure: A conditional entropy-based external cluster evaluation measure[END_REF]:

Homogeneity = 1 - P K i = 1 P C j = 1 n i' j n log n i' j P C j = 1 n i' j P C j = 1 P K i = 1 n i' j n log P K i = 1 n i' j n ' (2) Completeness = 1 - P C j = 1 P K i = 1 n i' j n log n i' j P K i = 1 n i' j P K i = 1 P C j = 1 n i' j n log P C j = 1 n i' j n ' (3)
The V-measure is defined as the harmonic mean of homogeneity and completeness. It is calculated as:

V -measure = 2 à homogeneity à completeness homogeneity + completeness : (4)
The clustering accuracy, F-measure, is defined as:

F -measure = X C j = 1 n j n max i F(i' j) (5)
2% and 4% for MetaCluster. These figures underlie the lower Adjusted Rand Index and F-measure values achieved by LSH-JP. On the other hand, the relabeling step in the LSH-SNN algorithm specifically aimed at reducing the number of unclustered sequences.

Table 1 indicates that LSH-SNN outperforms LSH-JP and MetaCluster in terms of homogeneity and Vmeasure, with the latter being the harmonic mean between homogeneity and completeness. LSH-SNN also slightly outperforms the other methods on the F-measure in four out of seven datasets, and it consistently yields the best performance on all the datasets in terms of the Adjusted Rand Index metric.

The LSH-SNN algorithm is implemented in the C++ programming language, and it uses the OpenMP application programming interface to support multiprocessing. Experiments were conducted on a Linux x86 64 server endowed with multi-core CPUs and 2 TB of RAM. The LSH-SNN and LSH-JP computations were parallelized on 48 cores.

Overall, these results demonstrate that LSH-SNN achieves accurate binning for DNA sequences as short as 600 bp, as compared with LSH-JP and MetaCluster's K-means implementation, and despite the latter being aware of the correct number of clusters (genomes) as an input parameter. The two highest peaks in Figure 1 correspond to the first clusters (shaded) of Table 2; these clusters contain the bulk of the pathogen genome.

EXPERIMENTS ON REAL DATASETS

To assess the behavior of the LSH-SNN method on real metagenomic data, we ran our algorithm on a set of 142,723 contigs resulting from a cross-assembly of 53 human microbiomes [START_REF] Alneberg | Binning metagenomic contigs by coverage and composition[END_REF], 43 of which harbored the Escherichia coli strain responsible for the 2011 Shiga toxin-producing E. coli (STEC) outbreak in Germany. We observed that (1) the number of output clusters was consistent with species diversity estimates derived from 16S ribosomal small subunit RNA analyses of the same samples [START_REF] Alneberg | Binning metagenomic contigs by coverage and composition[END_REF], (2) the pathologic E. coli strain genome was distributed among a very limited number of homogeneous clusters (Table 2), and (3) these pathogen-bearing clusters were the most discriminant variables to separate the individuals infected with the STEC strain from the remaining ones (Figure 1).

These observations demonstrate that the LSH-SNN algorithm performs rather decently on real-world datasets of biomedical relevance.

CONCLUSION

We have proposed an unsupervised composition-based method for binning large volumes of metagenomic sequences, without any prior knowledge of their reference genomes or the number of distinct genotypes present in the sample. LSH-SNN is based on two key steps: the hashing/indexing of the data space and the linking of sequences to build clusters. After computing the l-mer distribution of each sequence, LSH partitions the input space into buckets containing smaller subsets of sequences whose connectivity is evaluated via the SNN rule. A third step was added to reduce the number of singletons or unclustered sequences.

The LSH-SNN algorithm can scale to datasets containing millions of sequences, and it does not require the number of output clusters to be predetermined. Although the presented algorithm makes use of the SNN rule, we envision that the LSH concept could be combined with other sequence processing (including clustering) methods dealing with large data volumes, or it could be used on its own as exemplified in [START_REF] Berlin | Assembling large genomes with single-molecule sequencing and localitysensitive hashing[END_REF], where a MinHash LSH scheme was used to compute similarities between long noisy reads generated with a new single-molecule real-time sequencing technology.

The LSH-SNN algorithm was evaluated on seven semi-synthetic metagenomic datasets of different sizes and complexities (i.e., harboring different numbers of organisms/genotypes). We observed that LSH-SNN performs comparably or better on these datasets than the two other clustering algorithms tested (LSH-JP and K-means/MetaCluster). We should note, however, that, even though LSH-SNN significantly increases the size of the datasets that can be handled as compared with what can be achieved with the SNN method alone, its complexity does not compare favorably with Lloyd's heuristic underlying most K-means clustering engines. Therefore, the latter is probably more suited to the analysis of very large datasets containing billions of sequences, which are already generated nowadays from complex metagenomics samples (e.g., from soil) using high-throughput sequencing platforms. Alternatively, the LSH-SNN approach could be used to cluster contigs (sets of overlapping sequences) resulting from a preliminary (meta)genome assembly step, instead of being applied to raw (unassembled) reads as in this study.

FIG. 1 .

 1 FIG.1. The output clusters are arbitrary labeled with their identifier on the x axis. The y axis displays mutual information of the final clusters, which reflects the correlation of the clusters (treated as variables) with metadata labels (i.e., health status) of the individuals from the cohort.

Table 2 .

 2 Performance on Real Metagenomic Data

		No. of pathologic	Total no.
	Cluster_id	strain contigs	of contigs
	30	819	1082
	406	244	324
	33	9	24
	1	8	7527
	581	1	7

ACKNOWLEDGMENTS

This work was supported by the French Alternative Energies and Atomic Energy Commission (Commissariat l'Energie Atomique et aux Energies Alternatives), through its transversal program Technologies pour la Sante ´(MetaTarget project). The authors would like to thank Anestis Gkanogiannis for generating and sharing the semi-synthetic datasets, Dr. Christopher Quince for sharing the STEC outbreak contigs, and Marcel Salanoubat for overall support.

With

where n j is the cardinality of cluster C j .

Finally, the Adjusted Rand Index [START_REF] Hubert | Comparing partitions[END_REF][START_REF] Yeung | Details of the adjusted rand index and clustering algorithms, supplement to the paper: an empirical study on principal component analysis for clustering gene expression data[END_REF] computes a similarity measure between the observed and ideal clusterings as:

These measures all take values between 0 and 1, with higher values corresponding to better clustering accuracies.

Discussion

After tuning the parameters described in Section 4.2, we evaluated the performance of LSH-SNN as well as of two other clustering algorithms, LSH-JP and K-means (the MetaCluster implementation), on the different datasets (note that MetaCluster requires the number of component genomes to be provided as input for each dataset). Clustering accuracy was quantified by using the previously described metrics. The results are shown in Table 1 together with the rank (between parentheses) of each method, and they highlight (in bold) the best value for each evaluation criterion. The fraction of unclustered sequences (singleton elements) varied between 20% and 30% for LSH-SNN, 60% and 90% for LSH-JP, and between

AUTHOR DISCLOSURE STATEMENT

No competing financial interests exist.


Afifa Imtiaz
email: a.imtiaz@brgm.fr

Cecile Cornou
email: cecile.cornou@univ-grenoble-alpes.fr

Eleni Koufoudi
email: koufoude@univ-grenoble-alpes.fr

Spatial Variation of Seismic Ground Motions from European Rock-and Soil-site Dense Arrays

A significant part of spatial variation of seismic ground motions (SVSGM) over very short distances (e.g., dimensions of large engineering structures) is introduced by the effects of complex wave propagation and scattering as the seismic waves encounter heterogeneities at the shallow geology of the site. The stochastic models of SVSGM, expressed as coherency functions, account for the random phase variability of such non-uniform seismic excitations. However, the lack of association of these models with physical site-parameters renders them barely applicable in site-specific risk assessment of infrastructures. The present study, ongoing under the framework of the project EXAMIN (EXperimental Assessment and Modelling of ground motion spatial variability for performance based seismic risk assessment of industrial plants and INfrastructures; http://anr-examin.brgm.fr/) 2017-2019, funded by the French National Research Agency, attempts to shed light on this area of research. We present SVSGM estimated from low to moderate magnitude, local and regional earthquakes, recorded by different dense arrays of seismic stations in Europe. We studied one soft-rock and two soil arrays from the Cephalonia island of Greece, one soil array from the Fucino basin in Italy, along with one hard-rock array from the Saint-Guerin dam site and one soil array from the Grenoble basin in France. We quantified the coherency of shear-wave phases between two seismograms recorded by a pair of stations at 1-20 Hz frequencies. For each array and dataset, we averaged estimations from all pairs of stations having interstation distances falling into 10 m interval ranges between 10 and 200 m. Then, we examined the contribution of different physical parameters (e.g., wave characteristics, shear wave velocity) to the spatial variability of ground motion. Our observations indicate a largely site-dependent nature of spatial coherency and its strong correlation with characteristics of locally generated wave-field, likely to be associated with two-or three-dimensional geometrical effects and near-surface heterogeneities of the sites. Higher coherency is observed at rock sites, the highest being at the hard rock. However, at soil sites, the loss of coherency seems to be dominated by diffracted surface waves (Love and Rayleigh) propagating along the valley-edge directions. We eventually drew a comparison with the existing parametric models and found weaker correlations with our observed decay patterns. An analysis of three orthogonal directions of ground motions showed that the statistics of coherency is not sensitive to the orientation of the horizontal components. However, the vertical coherency, particularly at soil sites, is systematically lower than that of horizontal, and demonstrates an apparent increasing tendency at higher frequencies (>5 Hz) and shorter distances (<50 m), which has not been reported so far in the current literature. The results and understanding acquired through the analysis of these experimental data will be utilized in physics-based simulation of spatially variable ground motion in order to achieve the main goal of the project EXAMIN that is to develop a comprehensive approach for risk assessment of infrastructure networks.