Marine Bruneau

Thierry Mottet

Serge Moulin
email: serge.moulin@univ-fcomte.fr

Maël Kerbiriou

Franz Chouly

Stéphane Chrétien

Christophe Guyeux

Stéphane Chretien

A clustering package for nucleotide sequences using Laplacian Eigenmaps and Gaussian Mixture Model

Keywords: DNA clustering, Genomics, Laplacian Eigenmap, Gaussian mixture model

come

Introduction

As the amount of available genetic sequences increases drastically every year, accurate methods to deeply study them are strongly demanded [START_REF] Valot | Hocquet, What it takes to be a pseudomonas aeruginosa? the core genome of the opportunistic pathogen updated[END_REF]. Among these methods, clustering is a very powerful tool that helps to understand relations between sequences. It can be used, for instance, to classify 16S RNA sequences into OTUs [START_REF] Hao | Clustering 16s rRNA for OTU prediction: a method of unsupervised bayesian clustering[END_REF], which are standard proxies for microbial species (clustering is here a way to identify an a priori unknown number of "species"). Clustering is used too to define taxa within groups of species represented by their DNA sequences. Other utilizations of sequence clustering in genomics encompass the study of sub-populations within the same species [START_REF] Torroni | Native american mitochondrial dna analysis indicates that the amerind and the nadene populations were founded by two independent migrations[END_REF], the discovery of possible hidden variables that can explain differences between such sub-populations, and so on [START_REF] Suzek | Uniref: comprehensive and non-redundant uniprot reference clusters[END_REF].

However, most of the times, only generic methods for clustering a matrix of similarity scores is applied, and methods that are more specific to DNA sequences are still waited. Furthermore, the few existing methods that focus on sequence clustering mainly target on speed and scalability, and they need a strong similarity between sequences to produce accurate clusters. Furthermore, an a priori knowledge on the number of clusters is required, for instance by providing a similarity threshold cutoff. The objective of this article is to provide a new method that relaxes such constrains, allowing the sequences to be really divergent, and returning an a posteriori of the optimal number of clusters in an efficient manner.

One important issue in any clustering procedure is to find an appropriate embedding of the data under study, that will make the respective salient features in each of the groups more clearly delineateable. However, as often stated in the Machine Learning (ML) literature, high dimensional data are often much too scattered in order for an off-the-shelf method to be able to work properly. This phenomenon, often referred as the "curse of dimensionality" [START_REF] Bellman | Dynamic programming[END_REF][START_REF] Bellman | Adaptive control processes: a guided tour[END_REF], explains why several embeddings have been proposed recently. Some of these embeddings are very over-parametrized and can thus only be implemented in the supervised setting. This is the case for methods based on neural networks (e.g., auto-encoders). Other methods need less parameters and are more suitable to unsupervised learning, which is the case of our proposal. Among many nonlinear embedding methods, the Laplacian Eigenmap [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] approach has been extensively studied from both the theoret-ical and the application viewpoint [START_REF] Spielman | Spectral graph theory[END_REF]. More involved methods relying on Semi-Definite programming have also appeared recently with higher separation power in practice than spectral methods, see e.g., [START_REF] Chrétien | A semi-definite programming approach to low dimensional embedding for unsupervised clustering[END_REF].

In this article, the aim is to establish the practical efficiency, for DNA sequence clustering, of the combination of a plain Laplacian Eigenmap approach coupled with a Gaussian Mixture based clustering. This particular choice is motivated by its remarkable computational efficiency, even in the case where the objects to classify are really divergent. The overall procedure can be divided in three steps.

1. Compute a similarity matrix between each pair of DNA sequences, i.e., provide a matrix W of size n × n, where n is the number of sequences, which is such that W i,j increases with the "similarity" between sequences number i and j. 2. Diagonalise the Laplacian matrix of W . By such an operation, DNA sequences are mapped to elements of a given vector space, whose dimension is much smaller than the sequence lengths. This reduction of the problem dimension is a key element that usually has a great impact on both visualization and clustering. The combination of the two stages above is often referred as the "Laplacian eigenmap" approach. 3. Cluster the transformed data using a Gaussian Mixture Model (GMM [START_REF] Day | Estimating the components of a mixture of normal distributions[END_REF]).

By using a ready to use Python package designed at this occasion, and freely available on-line, we have demonstrated the accuracy and efficiency of this approach for DNA sequence clustering. Indeed, the methodology has firstly been tested on a sample of 100 ND3 genes (DNA sequences) from Platyhelminthes and Nematoda species that have been downloaded from the NCBI website [START_REF]National center for biotechnology information[END_REF]. The classification obtained via this approach has been compared with the phylogenetic tree of these species obtained by a likelihood maximization method using PhyML [START_REF] Guindon | PHYML online-a web server for fast maximum likelihood-based phylogenetic inference[END_REF]. Obtained clusters are consistent with both clades appearing in the phylogenetic tree and the NCBI taxonomy. In particular, this clustering perfectly separates the Nematoda and Platyhelminthes phyla.

To evaluate the method further, extensive simulation experiments have secondly been run on synthesized data: n DNA sequences have been randomly generated, and random mutations and block deletions have been applied on them, leading finally to N > n sequences. The objective was then to group these N sequences within n clusters (n is not a priori known), in such a way that all the elements in each cluster are originated from the same initial DNA sequence. On this set of synthesized data, our proposal has outperformed the two other state-of-the-art software for DNA clustering, namely Cd-hit-est and BLASTClust.

The remainder of this article is structured as follows. The three stages of the proposed method are detailed in the next section. Numerical results are then presented: the application example involving a real dataset is provided in Section 3.1, while the simulation based procedure is detailed in Section 3.2. A general discussion about the proposed approach is provided in Section 4. This research work ends by a conclusion section, in which the contribution is summarized and intended future work is provided.

The Clustering Method

Laplacian Eigenmap

The so-called Laplacian Eigenmap [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] is an original method for embedding data living in a structured set into a k-dimension vector space. The main information needed to compute the eigenmap is a matrix containing the value of the measured similarity between pairs of data. The main motivation for such embeddings is dimension reduction and visualization. Moreover, spectral methods often exhibit the nice property of separating clusters.

The matrix of similarity

The first step in the construction of a good embedding is the creation of a similarity matrix W . This matrix is intended to measure the similarity between each pair of sequences by providing a number ranging between 0 and 1. The main assumption on W is that the greater the similarity is, the closer are the sequences to each other.

In order to create this similarity matrix, a multiple global alignment of the DNA sequences is first run using the MUSCLE (Multiple Sequence Comparison by Log-Expectation [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF]) software. Then, an ad hoc "Needleman Wunsch distance" [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF] is computed for each pair of aligned sequence, and with the "EDNAFULL" scoring matrix. This distance takes into account that DNA sequences usually face (1) mutations and (2) insertion/deletion. Note that, by using MUSCLE as first stage of this matrix computation, we operate only one (multiple) sequence alignment, instead of n(n-1) the edit distance). Introducing Muscle leads to a real acceleration in the construction of the similarity matrix.

Let us call M the distance matrix obtained by this way. M is then divided by the largest distance value, so that all its coefficients are between 0 and 1. W can finally be obtained as follows:

∀ i, j ∈ [[1, n]], W i,j = 1 -M i,j ,
in such a way that W i,j represents the similarity score between sequences i and j.

Operations on W

Once the similarity matrix has been constructed, the next step is to create the normalized Laplacian matrix, as follows [START_REF] Chen | Resistance distance and the normalized laplacian spectrum[END_REF]:

L = D -1/2 (D -W )D -1/2 ,
where W is the similarity matrix defined previously and D is the degree matrix of W . That is to say, D is the diagonal matrix defined by:

∀i ∈ [[1, n]], D i,i = n j=1 W i,j .
L being symmetric and real, it is diagonalisable in a basis of pairwise orthogonal eigenvectors {φ 1 , ..., φ n } associated with eigenvalues 0 = λ 1 λ 2 ... λ n . The Laplacian Eigenmap consists in considering the following embedding function:

c k 1 (i) =      φ 2 (i) φ 3 (i) . . . φ k 1 +1 (i)      ∈ R k 1 ,
where c k 1 (i) is the coordinate vector of the point corresponding to the i th sequence. In other words, the coordinate vector of the point corresponding to the i th sequence is composed of the i th coordinate of each of the k 1 first eigenvectors, ordered according to the size of their eigenvalues. The choice of the k 1 cutoff is a crucial step and one usually proceeds as follows. The ordered eigenvalues are plotted, and we stop at the index where the increase in the eigenvalue is negligible: the number of eigenvalues that are not discarded is k 1 . For instance, in our program, we have chosen to set k 1 as the first time the difference between the k th and (k + 1) th value is lower than 0.01.

Note that W can be seen as a weighted adjacency matrix of a graph, where nodes are the DNA sequences while edges are labeled by the degree of affinity between their adjacent nodes. In the literature, the Laplacian matrix is often described as constructed from the weighted adjacency matrix of such a graph rather than constructed from a similarity matrix. These definitions are equivalent.

Gaussian Mixture based clustering

The final step is performed by applying Gaussian Mixture based clustering (GMM, [START_REF] Day | Estimating the components of a mixture of normal distributions[END_REF]) to the point cloud. Gaussian Mixture Models belong to the class of unsupervised learning schemes [START_REF] Friedman | The elements of statistical learning[END_REF], and allow to distribute the data points into different clusters without a priori assumption about the clusters interpretation. One of the very useful features of model based clustering is that the model allows to use information criteria in order to estimate the number of clusters using Akaike Information Criterion (AIC [START_REF] Akaike | A new look at the statistical model identification[END_REF]), Bayesian Information Criterion (BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF]), or Integrated Completed Likelihood (ICL [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]) well-known criteria. The mathematical assumption of a GMM is that the point cloud follows the distribution:

k 2 i=1 τ i N (µ i , Σ i ),
where k 2 is the number of clusters, τ i is the probability for a point to be in cluster i, and N (µ i , Σ i ) is the normal distribution of mean µ i and covariance matrix Σ i . GMM parameters are computed with the Expectation-Maximization (EM) algorithm [START_REF] Mclachlan | Finite mixture models[END_REF]. Notice that the EM algorithm may converge to singular distributions exponentially fast [START_REF] Biernacki | Degeneracy in the maximum likelihood estimation of univariate gaussian mixtures with EM[END_REF]. However, degenerate situations can be easily discarded and consistent estimators can be easily obtained in practice. Gaussian Mixture models are still a topic of current extensive research, both from the statistical perspective [START_REF] He | Laplacian regularized gaussian mixture model for data clustering[END_REF][START_REF] Wang | High dimensional expectationmaximization algorithm: Statistical optimization and asymptotic normality[END_REF] and the computational one [START_REF] Yi | Regularized EM algorithms: A unified framework and statistical guarantees[END_REF].

The Bayesian Information Criterion has been chosen to determine the most relevant number of clusters k 2 to be considered. The BIC, which is a criterion for model selection, is defined as follows:

BIC = -2 ln(L) + ln(n)p,
where L is the likelihood of the estimated model, n is the number of observations in the sample, and p is the number of model parameters. This criterion allows us to select a model whose validity is based on a balance between the value of the model's likelihood (fidelity term) and the number of parameters to estimate (complexity term). The likelihood of the model increases with k 2 as well as the number of parameters. The selected model will be by default the one that minimizes this criterion. In the proposed package, the user can also set the number of clusters manually.

The clustering software

The Python program corresponding to the algorithm described in this section is freely available online1 . The main function of this package provides a clustering from nucleotide sequences. Its prototype meets the following canvas: clustering = Gclust(liste, nbClusters='BIC', drawgraphs=True, nbEVMethod = 'delta', nbEVCutOff = 'default', AddToNamesOfOutputs = ''): where:

• list is an in-memory image of a fasta file containing the sequences associated with their names. The fasta file must be configured as the "ND3.fasta" file in our github repository.

• nbClusters is the number of clusters desired by the user. By default, the program applies the BIC criterion to determine it. The user may also choose to use the AIC criterion by writing "nbClusters = AIC".

• drawgraphs is an optional Boolean value to produce some graphics. If drawgraphs = TRUE, a two dimensional clustering of data is plotted (cf. Figures 4,5, and 6), as well as the graphical representation of similarities (as in Figure 1).

• nbEVMethod is the method chosen to determine the number of considered eigenvectors k 1 . The user can choose between 3 methods, usually reported in the literature.

-If nbEVMethod = 'delta' then k 1 is the lowest value such that λ k 1 +1 -λ k 1 < δ where δ is a constant to be fixed by the user. 'delta' is the default method, and with δ = 0.01.

-If nbEVMethod = 'energy' then k 1 is the lowest value such that

k 1 i=1 (λ max -λ i ) ≥ C × n i=1 (λ max -λ i )
where C is a constant to be chosen by the user (default C = 0.9)

-If nbEVMethod = 'log' then k 1 is the rounded value of log(n).
• nbEVCutOff is the constant δ if nbEVMethod == 'delta' or the constant C if nbEVMethod == 'energy'. By default, nbEVCutOff = 0.01 if nbEVMethod == 'delta' and nbEVCutOff = 0.9 if nbEVMethod == 'energy'.

The GClust function can be, for instance, launched as follows within a Python script: list = open('MyFileName.fasta').read() groups = gclust(list, drawgGraphs=True) The output is a list of k 2 lists, each list containing the references of the sequences grouped in a particular cluster.

Module and package dependencies

As explained in Section 2.1.1, the DNA sequence alignment software used during the similarity matrix stage is MUSCLE. More precisely, the "muscle v38" function of cogent package has been used [START_REF] Knight | Pycogent: a toolkit for making sense from sequence[END_REF]. The Gaussian Mixture Model, for its part, is performed using the GMM function of sklearn.mixture package [START_REF] Buitinck | Api design for machine learning software: experiences from the scikit-learn project[END_REF].

Numerical evaluations

Evaluation on real genomic data

The proposed method has been applied to a sample of 100 DNA sequences from the mitochondrial gene ND3 taken from various species of both Platyhelminthes and Nematoda. Figure 1 graphically displays the similarity matrix that was obtained from the data. As expected, the largest similarities are obtained on the diagonal, that is, when the sequences are compared to themselves. However, blocks seem to appear as well. This justifying the prior intuition that the embedding method separates the clusters if there happens to be more than one in the dataset. Figure 3 shows the Bayesian Information Criterion of the Gaussian Mixture Models for various number of clusters. One can see that this BIC reaches 4, to the one formed by the first and third eigenvectors in Figure 5, and to the plane formed by the second and the third eigenvectors in Figure 6. In these graphs, clusters 0, 1, 2, and 3 are represented in blue, red, yellow and cyan respectively. A FASTA file was then written, in which each nucleotide sequence has been labeled according to its taxonomy and its cluster number. The tax- onomies have been found thanks to the efetch function (sub-package Entrez2 , Python package Bio). We then have used the online software PhyML (maximum likelihood method for phylogenetic tree reconstruction) with default options [START_REF] Guindon | PHYML online-a web server for fast maximum likelihood-based phylogenetic inference[END_REF], to build a tree based on the same mitochondrial ND3 gene that we have used during clustering. The tree, whose leaves contain taxa names and cluster id, has been displayed using FigTree software [START_REF] Morariu | Automatic online tuning for fast gaussian summation[END_REF]. This tree is depicted in Figures 7 and8.

The efetch function provides 8 taxonomic levels for each species in our sample. In Figures 7 and8, for readability reasons, we only show the taxonomy between levels 4 and 6. Note first that this clustering perfectly separates Platyhelminthes and Nematoda phyla. Indeed Trematoda and Cestoda are two classes of Platyhelminthes. Cluster 1 corresponds exactly to the Platyhelminthes, while Clusters 0, 2, and 3 represent the Nematoda. Cluster 0 is only composed of Spirurida, it contains 10 of the 12 members of this taxon. More precisely, when considering the seventh taxonomic level, we found that Cluster 0 contains nine Filarioidea and one Thelazioidea, while the Spirurida that are not in this cluster are a Dracunculoidea and a Physalopteroidea. Cluster 3, for its part, corresponds exactly with the Trichocephalida taxon. Finally, in addition to recover taxonomy, the clustering agrees very well with the tree obtained via PhyML, as shown in Figures 7 and8. Note that the taxonomy, based on morphology and nuclear genome, fully agrees with the PhyML tree and our clustering, which are both based on the mitochondrial genome. This fact suggests that this mitochondrial gene evolves like the nuclear genome.

Tests on simulated data

In addition to the analysis of real data processed in the previous section, we also performed tests on simulated data. This new series of tests allow us to assess the accuracy of our tool and to compare this accuracy with other existing tools. The data are simulated according to the following steps:

• We generate a "global root". That is to say, a random word of 500 letters {A,T,C,G} representing a sequence of 500 nucleotides. • From this global root, we generate 3 "cluster roots" which are each an evolution of the global root for which each nucleotide has mutated with probability 0.05.

• From each cluster root, we generate a cluster of 10 sequences. Each sequence is an evolution of the cluster root for which each nucleotide has mutated with probability 0.05.

• Each of the obtained sequences undergoes a block deletion of 2.5% to 7.5% of its nucleotides with a probability of 0.5. That is to say, for each sequence:

-A random draw is performed in order to decide on whether a block deletion is carried out or not. The deletion is carried out with probability 0.5.

-A uniform random draw is performed to determine how many nucleotides are deleted between 2.5% and 7.5% of the sequence size.

-A uniform random draw is performed in order to decide on where the block deletion starts.

The procedure defined above makes it possible to generate a simulated base of 30 sequences. We repeatedly applied this procedure 25 times with different random seeds (seed = 0 ... 24) in order to simulate 25 different databases. Figure 9 shows the similarity matrix obtained via these simulations (with a random seed of 0). The clusters are sorted in each of the simulated databases (as can be seen on the similarity matrix). Thus, the ideal outcome we would like to obtain in applying our tool to these databases is, for example:

[['S11', 'S12', 'S13', 'S14', 'S15', 'S16', 'S17', 'S18', 'S19', 'S20'], ['S21', 'S22', 'S23', 'S24', 'S25', 'S26', 'S27', 'S28', 'S29', 'S30'], ['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10']] In order to evaluate the performance of our tool to find the right clustering, we have created a distance that corresponds to the number of sequences that need to move to get the perfect clustering. For instance, the following clustering:

[['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10'],

['S21', 'S23', 'S24', 'S25', 'S26', 'S27', 'S28', 'S29', 'S30'],

['S11', 'S12', 'S13', 'S14', 'S15', 'S16', 'S17', 'S18', 'S19', 'S20'], ['S22']] has a distance from the perfect clustering equal to 1.

Table 1 summarizes the distances obtained with our tool for each of the simulated database. This table also summarizes the distance obtained with the "CD-Hit-est" and "BlastClust" tools.

Tests on simulated data with other tools

"CD-Hit-est" and "BLASTClust" are two of the most used clustering tools for nucleotide sequences. According to [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences[END_REF]: "Cd-hit-est is a greedy incremental clustering algorithm. Briefly, the sequences are first sorted in decreasing length. The longest sequence becomes the representative of the first cluster. Then, each remaining sequence is compared with the representatives of the existing clusters. If the similarity with any representative is above a given threshold, it is grouped into that cluster. Otherwise, a new cluster is defined with that sequence as the representative." And, according to [START_REF]Using blastclust to make non-redundant sequence sets[END_REF]: "The program [BLASTClust] begins with pairwise matches and places a sequence in a cluster if the sequence matches at least one sequence already in the cluster. In the case (...) of nucleotide sequences, the Megablast algorithm is used."

In both cases, the user has to provide a similarity threshold. In the case of CD-hit-est, this threshold indicates the minimum similarity that a sequence must have with the reference of a cluster to integrate this cluster. In the case of BlastClust, this threshold indicates the minimum similarity that a sequence must have with at least one other sequence in the cluster. This necessity to specify the threshold for similarity makes a important difference with our tool. Indeed, our tool determines the number of clusters automatically using BIC. value for this similarity (e.g., 0.9 for CD-hit-est) but it is simple to find situations in which this value is not adapted. Indeed, this value does not fit the data as BIC does.

We have not found a way to choose similarity thresholds based only on the simulated databases. Thus, to apply CD-hit-est and BLASTClust to the simulated data we have deliberately biased our inputs in favor of these programs: we have sought similarity thresholds that minimize the distance to the ideal clustering (as if this ideal clustering was known in advance) in the case of the first simulated database. Then we have applied these thresholds to the 25 simulated bases.

In the case of CD-hit-est, we have tested different similarities between 0.82 and 0.90. Indeed, we have calculated that, given the parameters of the simulations, the ideal threshold should be in this area. The results are shown in Table 2. On can see that the distance to the ideal clustering is minimized for a similarity threshold of 0.84. When studying these CD-hitest clustering in detail, one can see that, even with this ideal similarity of 0.84, some sequences are isolated and do not fit into their ideal cluster. If the chosen similarity threshold increases, this number of isolated sequences increases. On the other hand, if the chosen similarity threshold decreases (below 0.84), some clusters merge together.

The case of BLASTClust is a little more complicated. In this case, we have varied four parameters detailed below:

• -S <threshold> similarity threshold if <3 then the threshold is set as a BLAST score density (0.0 to 3.0; default = 1.75)

if ≥ 3 then the threshold is set as a percent of identical residues (3 to 100)

• -L <threshold> minimum length coverage (0.0 to 1.0; default = 0.9)

• -b <T-F> require coverage as specified by -L and -S on both (T) or only one (F) sequence of a pair (default = TRUE)

There are thus two binary parameters (-b on the one hand and the fact that S < 3 or S ≥ 3 on the other hand) and two continuous parameters (S et L).

The two binary parameters constitute 4 configurations (b = 'T 'and S < 3 ... b =' F' and S ≥ 3). For each of these 4 cases, we have tested BlastClust on a 10,000 points grid (100 possibilities for S × 100 possibilities for L). We have observed the minimum distance to the ideal clustering obtained on these 40,000 points (it is 3). Then we have counted which of the four configurations contain the largest number of points that minimize the distance to the ideal clustering (it is b = 'F' and S < 3). After this, we have selected the point of this configuration closest to the average of the points that minimize the distance to the ideal clustering among the points that minimize the distance to the ideal clustering. The setting obtained is: S = 1.91 (BLAST score density), b = 'F', L = 0.49 and p = 'F' (p = 'F' means that we work on nucleotide sequences, not on proteins).

After we had obtained these parameters for CD-hit-est and BLASTClust, we finally applied these two tools to our [START_REF] Knight | Pycogent: a toolkit for making sense from sequence[END_REF] 1. It uses a statistical criterion (in this case, the reputed BIC) to determine the number of clusters to consider. Thus, it can propose a number of clusters adapted to the database without the user having to provide any a priori. 2. As shown in Section 3.2, it allows a better reconstruction of the ideal clustering, even when the similarity threshold is chosen advantageously for CD-hit-est and BLASTClust. 3. It allows the user to plot useful graphical representations of the clustered data.

Finally, we note that CD-hit-est works only for intra-cluster similarities larger than 75%. One drawback of our tool as compared to CD-hit-est and BLASTClust is the computation speed. According to [START_REF]Cd-hit user's guide[END_REF], CD-hit-est is the fastest of these two programs. We have not been able to rigorously confirm this remark in our study, simply because we used BLASTClust on our computers while we used CD-hit-est online.

To sum up, the main features of our tool are different from that of CD-hitis and BLASTClust. For users wanting to cluster their dataset into meaningful taxa, which makes it possible to apprehend the evolution, our tool might be of greater interest. On the other hand, if the objective is to reduce the size of the dataset by removing duplicates (or retain only one sequence per group of close sequences), it is more appropriate to use one of the other tools with a similarity of 100% (or close to 100%). In particular, CD-hit-est is written so as to provide a representative sequence for each cluster.

Possible alternatives with the same caneva

Various options are possible to perform the analysis we have presented previously, some of them being listed below.

Similarity matrix

As stated previously, the multiple global alignment step is performed first before computing similarities, using MUSCLE software [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF]. Among the most extensively used methods, we can choose MAFFT [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF] too, as well as ClustalW or ClustalX [START_REF] Larkin | Clustal w and clustal x version 2.0[END_REF]. In addition, instead of defining the similarity matrix W as W i,j = 1 -M i,j , it could be possible to consider W i,j = 1 M i,j or W i,j = e -M i,j .

Number of considered eigenvectors

The number of eigenvectors to keep is another point to investigate. As explained in Section 2.1.2, it is usually advised to check graphically when the increase of eigenvalues is reducing. In this article, we have chosen to consider as default proposal: k 1 such that δ = λ k 1 +1 -λ k 1 < 0.01. This criterion has led to k = 4 in the case study, which seems acceptable according to the considered taxonomy.

Some authors in the literature proposed to compute k 1 as the logarithm of n [START_REF] Matias | Notes de cours: Analyse statistique de graphes[END_REF]. This method has been implemented as an option, as specified in Section 2.3. In addition, the "energy" method defined in Section 2.3 has been implemented too, in order to propose something similar to what is done in usual Principal Component Analysis (PCA). Other approaches can be considered to solve this problem, which is still an open one.

Number of clusters

We have chosen to consider the BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] to determine the optimal number of clusters, which is a common choice for this type of problem. An alternative may be to use the Akaike Information Criterion (AIC, [START_REF] Akaike | A new look at the statistical model identification[END_REF]). The principle of calculating the AIC is the same as the BIC, since the goal is to maximize log-likelihood penalized by the number of parameters (or more precisely, to minimize the number of parameters to which the log-likelihood is subtracted). AIC formula is the following:

AIC = -2ln(L) + 2 × p,
where L is the likelihood of the estimated model and p the number of model parameters.

When ln(n) ≥ 2, that is to say n ≥ 8 (which is always the case in practice), BIC penalization is larger than that of AIC. Thus the number of clusters obtained by BIC is lower or equal to the one obtained by AIC. BIC is said to be "more conservative" than AIC. The choice between these two criteria can be dependent on how stringent the clustering is desired. The user of the GClust function may choose to use the AIC rather than the BIC as specified in Section 2.3. The user of the GClust function is also able to choose the number of clusters of his or her choice.

Conclusion

In this work, we have proposed a new method of nucleotide sequence clustering. This clustering is produced by a methodology combining Laplacian Eigenmap with Gaussian Mixture models, while the number of clusters is automatically determined by using the Bayesian Information Criterion. The proposed methodology was applied to 100 sequences of mitochondrially encoded NADH dehydrogenase 3. The resulting clusters appeared to be coherent with the phylogeny (gene tree obtained with PhyML) as well as with the NCBI taxonomy. In addition, further tests have also been carried out on fully simulated data. These tests showed that our methodology allows to recover the expected clusters with greater accuracy.

One possible extension for future work could be to investigate more deeply the impact of parameters in the obtained clusters. The effects of using a different similarity matrix, or choosing a different dimension of the image space in Laplacian Eigenmap, or the number of desired clusters, could be investigated. Moreover, our tests on real data allowed us to watch our tool in action while performing a clustering of species into taxons. It might be interesting to also test the ability of our tool to classify 16S RNA sequences into OTUs on real data. Another avenue could consist in adapting the code to the very similar problem of protein clustering.

Figure 1 :

 1 Figure 1: Similarity matrix

Figure 2

 2 Figure2displays the 14 first ascending eigenvalues, labeled λ i , and satisfying 0 = λ 1 ≤ λ 2 ≤ ... ≤ λ 14 . According to the proposed methodology, we have chosen k 1 as the minimal index k such that λ k+1 -λ k is lower than or equal to 0.01. In this case study, we found k 1 = 4.Figure3shows the Bayesian Information Criterion of the Gaussian Mixture Models for various number of clusters. One can see that this BIC reaches

Figure 2 :

 2 Figure 2: Curve representing the first 14 eigenvalues

Figure 3 :

 3 Figure 3: Bayesian Information Criterion of the Gaussian Mixture Models

Figures 4 , 5 ,

 45 Figures 4, 5, and 6 represent the point cloud divided in 4 clusters. This point cloud is projected into the plane formed by the first and second eigenvectors in Figure4, to the one formed by the first and third eigenvectors in Figure5, and to the plane formed by the second and the third eigenvectors in Figure6. In these graphs, clusters 0, 1, 2, and 3 are represented in blue, red, yellow and cyan respectively.

Figure 4 : 2 Figure 5 :

 425 Figure 4: GMM clustering in the plane formed by the eigenvectors 1 and 2

Figure 6 :

 6 Figure 6: GMM clustering in the plane formed by the eigenvectors 2 and 3

Figure 7 :

 7 Figure 7: First part of the phylogenetic tree (Platyhelminthes)

Figure 8 :

 8 Figure 8: Second part of the phylogenetic tree (Nematoda)

Figure 9 :

 9 Figure 9: Similarity matrix of the simulated clusters (seed = 0)

Table 1 :

 1 Of course, CD-hit-est and BLASTClust provide a default Distance from the perfect clustering

	random Distance Distance	Distance
	seed	our tool CD-hit-est BlastClust
	0	0	7	3
	1	0	11	19
	2	2	7	1
	3	10	13	20
	4	0	8	3
	5	0	7	5
	6	0	8	2
	7	0	7	9
	8	0	7	12
	9	1	10	19
	10	0	6	5
	11	0	7	10
	12	1	6	5
	13	1	10	2
	14	0	7	9
	15	0	7	5
	16	0	6	4
	17	0	11	17
	18	0	12	16
	19	0	6	17
	20	0	8	17
	21	10	8	17
	22	0	7	1
	23	0	7	2
	24	0	8	6
	total	25	201	226

  simulated datasets. The results are shown in Table1with those of our tool.

	Similarity Intra-Cluster	0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.90
	Distance to Ideal Clustering 10	8	7	9	9	9	10	10	13

Table 2 :

 2 Search for the best similarity threshold for CD-hit-est

	4. Discussion
	4.1. Comparison with other tools
	Compared to CD-hit-est and BLASTClust, our tool has several advan-
	tages.

(pairwise) alignments in the classical Needleman Wunsch algorithm (that usually contains two stages: finding the best pairwise alignment, and then compute

https://github.com/SergeMOULIN/clustering-tool-for-nucleotide-sequences-using-Laplacian-Eigenmaps-and-Gaussian-Mixture-Models

A package that provides code to access NCBI over the world wide web.


Paul Langevin

G Pascoli
email: pascoli@u-picardie.fr

The Sagnac effect and its interpretation by Paul Langevin L'effet Sagnac et son interpretation by

The French physicist Georges Sagnac is nowdays frequently cited by the engineers who work on devices such as ring-laser gyroscopes. These systems operate on the principle of the Sagnac effect. It is less known that Sagnac was a strong opponent to the theory of special relativity proposed by Albert Einstein. He set up his experiment to prove the existence of the aether discarded by Einsteinian relativity. An accurate explanation of the phenomenon was provided by Paul Langevin in 1921.

Résumé

Le nom de Georges Sagnac est aujourd'hui très connu des ingénieurs travaillant sur les systèmes de navigation, tant maritime qu'aérienne, qui exploitent sa découverte de 1913. Ce que l'on sait moins est que ce physicien était un farouche opposant à la relativité, théorie révolutionnaire qui avait été élaborée par Albert Einstein quelques années auparavant en 1905. L'expérience de Sagnac a été pensée pour prouver l'existence de l'éther luminifère. Son interprétation correcte (relativiste) a été fournie par Paul Langevin en 1921.

1.

Historical development

It would seem useful to recall the context of the the discovery of the Sagnac effect. Around 1910, a very strong conservatism reigned concerning the paradigm of absolute space and its counterpart, the aether, the hypothetical medium of propagation of light. To this statement we must add the presence of a very large number of sceptics, both physicists and philosophers, who opposed a non-sequitur to the new vision of space and time that was implied by the very recent theory of Einsteinian relativity. Some key points should also be highlighted.

Newton and the absolute space (end of the 17th century)

For Isaac Newton the concept of absolute space (or pure extension) exists, regardless of the presence or not of any material. In Newton's worldview, absolute space was a pre-existing frame, empty and immutable, a fundamental background to any moving body. But its undefined nature made it a kind of metaphysical entity, close to transcendence. In the Scholium, at the beginning of the Principia (Definitiones), Newton expatiated on time, absolute or relative space and motion [START_REF] Newton | Philosophiae Naturalis Principia Mathematica, Definitiones, Scholium[END_REF]. He suggested that an absolute motion can be distinguished from a relative motion. To justify his argument, he introduces a thought experiment that is still famous. A bucket filled with water is subjected to a fast rotation. There are two phases : i. The bucket rotates but the surface of the water remains flat, ii. The movement is finally transmitted to the water. It is then observed that the latter leaks from the center and ascends towards the walls, the free surface taking the form of a paraboloid. During the first phase there is indeed a relative movement water/bucket, but no centrifugal force appears. It is only in the second phase, when the movement of rotation is communicated to the water, as these forces develop. For Newton the conclusion was that these forces are not correlated with the relative movement of water/bucket, but are generated by the movement of water in relation to absolute space. A second thought experiment was also devised by Newton. Two spheres are connected by a rope and the assembly is rotated. The rope tightens and one can measure the tension that develops there. Newton's conclusion is the same.

Ernst Mach and the relativity of motion (end of the 19th century)

The concept of absolute space as introduced by Newton was strongly criticized by Ernst Mach in his book published in 1893 [START_REF] Mach | The Science of Mechanics[END_REF]. Mach assumed a positivist vision from the outset. He noted that Newton's experiment simply demonstrates that the shape of the water surface (a paraboloid) is not induced by the relative movement of water in relation to the immediate environment (the bucket in particular). But Newton did not extend his thought sufficiently far. It did not take account the distant masses (represented by the stars). According to Mach it is legitimate to attribute the deformation of the surface of the water to a rotation, but not relative to an absolute space as Newton had suggested, but to a relative rotation of the water and massive distant bodies. Suppose that the bucket were fixed and that the whole Universe put in rotation. Would the surface of the water take the form of a paraboloid ? For Mach the answer is yes. We know that Mach's ideas had a considerable influence on the development of Albert Einstein's reasoning, especially during the first five years of the 20st century. But most physicists did not share Mach's point of view, and were very opposed to getting rid of the concepts of absolute space and aether.

The luminiferous aether

The idea of a medium for the propagation of light waves goes back to Christian Huygens (1690). Unfortunately owing to Newton's great aura, Newton's corpuscular theory, exposed in his treatise Opticks (1704), completely dominated for more than 100 years and took precedence over Huygens' wave description. But Huyghens hypothesis enjoyed a sudden revival at the very beginning of the 19th century, thanks to Augustin Fresnel and Thomas Young. Since then the aether was considered as the absolute frame of reference, especially in relation to which Maxwell's equations (dating from 1862) had to be written and the celerity of the light was the constant c. This aether was automatically identified with the absolute space of Newton. But the Earth moves in this absolute space. Does it carry the aether with it as it moves ? To answer this question considered as fundamental by the physicists at that time, a number of experiments were carried out in the first half [START_REF] Arago | [END_REF]4] and at the end of the 19th century [5,6]. In spite of the fact that all these experiments, taken together, showed that the aether was an illusive medium, it must be remarked that as late as 1910, Hendrick Lorentz and Henri Poincaré had still not adopted a very clear position concerning the status of the aether [START_REF] Lorentz | Das Relativitätsprinzip[END_REF]. For a more complete analysis of the debate on the light-bearing aether, see [START_REF] Darrigol | [END_REF][START_REF] Martinez-Chavanz | L'Expérience de Sagnac et le disque tournant[END_REF][START_REF] Lalli | [END_REF][START_REF] Lalli | Styles of Thinking in Science and Technology[END_REF].

The experiment of Sagnac

In 1899, Georges Sagnac 1 had developed a theory of the existence of a motionless mechanical aether [START_REF] Sagnac | [END_REF]. His aim was to explain within this theoretical framework all the optics phenomena, and especially the Fresnel-Fizeau experiment for the drag of light in a moving medium [4]. In 1910, he conceived a rotating interferometer for testing his ideas and to see the optical whirlwind effect, as he called it (Figure 1).

The effect discovered by Sagnac, and published in the Comptes Rendus de l'Académie des Sciences in 1913 [13,14] 2 , has led to many applications in positioning technologies (Global Positioning System, Galileo, etc.). The sensors (gyrolasers and gyrometers) that measure angular velocities with respect to an inertial frame of reference, exploit the physics of the Sagnac effect. Generally numbering three, these sensors form a part of the inertial stations, in particular those used to determine the orientation of all types of vehicles : aircraft, boats, submarines, space probes, etc. These devices have also been applied to geodesy and sismology [16,17,18]. Since the 1960s, a large number of studies have been devoted to the narrative of the discovery of this effect, its interpretation and to its practical consequences [19,20,21,22]. On 10 October 2013, a symposium held at the Fondation Simone-et-Cino-del-Duca (Paris) was organized to celebrate the centenary of the Sagnac effect [23]. Even more and ironically enough, the Sagnac device initially conceived to prove that the special relativity is wrong, is today used for testing the general relativity. For instance, in the GINGER (Gyroscopes IN GEneral Relativity) project a three-dimensional array of large-size ring-lasers will be conceived for measuring in a ground-based laboratory the de Sitter and Lense-Thirring effects [24]. Let us consider a rotating disk on whose rim two light beams circulate in opposite directions relative to each other 3 . Both beams initially start from the same point, located on the rim of the disk. In order to simplify the reasoning but without loss of generality, it will be assumed that the source and the receiver are combined and both attached to the disk. Let R be the radius of the disk and ω its angular velocity. This velocity is assumed to be 2. A more extended version has also been published by Sagnac in the Journal of Physics the following year in 1914 [15].

3. The generalization to any closed (not-circular) loop is immediate.

constant and the disk rotates in the counterclockwise sense for the inertial observer, L linked to the laboratory. For this observer, the receiver travels the distance ∆l (circular arc) at the speed v = Rω. The beam in the sense (+) must therefore travel the distance l + ∆l before to reach the receiver. The duration of the path of the light ray, t + , is deduced from the equality (Figure 2)

t + = l + ∆l c = ∆l v (1)
or t + = l c-v . The same reasoning applies to the other ray, moving in the opposite sense (-), and the traveled distance is l -∆l. The journey time is t -= l c+v . The difference in the journey times is therefore

∆t L = t + -t -= 2lv c 2 -v 2 = 4πR 2 ω c 2 (1 - R 2 ω 2 c 2 ) -1 4Aω c 2 (2)
where A = πR 2 denotes the area of the disk. The difference in duration of travel induces a phase shift between the two waves, hence the appearance of interference fringes (Sagnac made a differential measurement by rotating the disk in one direction and then in the other). The phase shift (expressed in wavelength) is given by the Sagnac formula

∆ϕ = 8πAω λc (3) Figure 2
What is the reasoning of a non-inertial observer (D) related to the disk ? According to Georges Sagnac, the light has a medium of propagation which is the aether. In the reference frame linked to the laboratory (observer L) the propagation velocity (celerity) is c. On the other hand, for an observer linked to the rotating disk the light speed propagates differently, according to whether the beam propagates counterclockwise (the direction of rotation of the disk by convention), or in reverse. It is then respectively c -v and c + v (Figure 3). For the observer D, we have for the difference in time of course

∆t D = l c -v - l c + v 4Aω c 2 (4)
Figure 3 in perfect agreement with the result found by the observer L. Numerically the area of the polygonal interferometer was 0.0860 m 2 , the rotation rate of order 2 Hz and the resulting fringe δϕ = 0.07 ± 0.01. In his paper, Sagnac also predicts that a similar effet would be observed from the diurnal Earth rotation (the Earth taken as a rotating turntable) using a very large interferometer 4 . For Sagnac this result proves the existence of the aether. The argument of this author appears a priori quite convincing. The existence of a propagation medium for light waves seems to be confirmed by his experiment. But in science an hypothesis cannot be validated only if it allows to explain all the experiments, not one in particular. Other experiments definitively invalidate the hypothesis of the aether. Perhaps most famous is the crucial experiment of Michelson and Morley [5]. The aim of this experiment was to measure the displacement of the Earth with respect to the aether. Using an interferometric device, they showed that the velocity of the aether wind, if any, was less than 5 km/s (we should remember that the velocity of the Earth in the solar system is 30 km/s, and therefore much greater). The procedure has obviously been repeated several times within a year (in case at a given time, and by an extraordinary coincidence, the 4. Let us note that the so-called Sagnac effect had been approached by Sir Oliver Lodge as early as 1893 [25]. This effect had also been glimpsed by Harres in 1911 [26], but this author does not seem to have paid much attention to it. When Sagnac realized his experiment, the idea was therefore in the air.

Earth's velocity at the time of measurement had been zero with respect to the aether), but no seasonal effects have been found. The Michelson-Morley experiment has been repeated with ever-increasing accuracy. The Michelson-type experiment carried out by Georg Joos in 1930 allowed him to give a bound to the velocity of the aether wind of less than 1.5 km/s [27]. More recently experiments have been carried out using various sophisticated devices, such as lasers and masers, optical resonators and microwaves, etc. [28]. The conclusion is always the same : no measurable aether wind. We can thus conclude that the relative velocity of the Earth with respect to the aether is experimentally found to be very close to zero, with high precision. But if we only focus on Michelson-Morley's experiment, it is still perfectly legitimate to imagine that the Earth completely drags the aether in its translatory motion (Fresnel had already suggested in 1818 the hypothesis of a partial dragging of the aether with moving substances [29]). In this case the negative result of the Michelson experiment has a trivial explanation. Unfortunately however, such an assumption cannot be maintained. There are two immediate reasons for this. First of all, it is inconsistent with the aberration of fixed stars (as we know, during a year the stars describe a small ellipse on the background of the sky. This effect cannot occur if the aether is fully dragged by the Earth). Secondly the experiment of Sagnac was repeated by Michelson and Gale, in 1925 [30], but this time taking the Earth as a rotating disk (as already suggested by Sagnac himself). These authors observed a displacement of the fringes of interferences, as had Sagnac in his own experiment. This positive result undoubtedly confirms that the Earth does not drag the hypothetical aether in its rotation (it is therefore illogical to admit that it drags this medium in its translation). The only acceptable conclusion that can be drawn from these two experiments, Michelson-Morley, on the one hand, and Michelson-Gale on the other hand, is that the hypothesis of the existence of a medium of propagation for light is not tenable (unlike sound, the light can spread in the vacuum). In the classical context, it is clear that the Sagnac effect cannot at all be explained.

The relativistic interpretation

Special relativity has successfully passed tests of thousands of experiments of all types (for example, the experiments of Kennedy -Thorndike [31] and Ives -Stilwell [32]). If it is clearly no longer questionable today, the situation was quite different at the beginning of the twentieth century 5 . An experiment such as Sagnac's was still used as a strong argument by the opponents of the theory of relativity, who were still quite numerous at the time, including Sagnac himself. The scientific community as a whole only accepted relativity after some time. So a lot of pseudoparadoxes issued from relativity, which seemed to arise here and there, were solved rather slowly. The first convincing explanation of the Sagnac effect, in the framework of relativity, was provided by Paul Langevin in 1921 [36] 6 . It appears rather amazing that the correct relativistic interpretation of the Sagnac effect took eight years. A seemingly obvious reason is that Sagnac's experiment was not very much discussed in the scientific literature even in France after the discovery of 1913 7 . Conscious of this situation, in 1919, Sagnac published five papers on his work in the Comptes Rendus [37]. The paradox is that his ideas were nevertheless borne by a French group of strong antirelativists. In 1919, Sagnac was even rewarded with the Pierson-Perrin Prize for his achievements on this topic (first for the experiment, seen as a rebuttal of the relativity principle, the constancy of light, and also for having proven the reality of absolute space and time) [38]. Eventually, in 1921, the mathematician (and relativist skeptic) E. Picard asked P. Langevin to demonstrate the effect in the framework of the relativity, perhaps guessing it was impossible [39]. Again any such expectations and within a couple of months, Langevin derived the effect from the General relativity [36]. Let us specify however that the sagnac effect was already derived from Special relativity by Von Laue one year before Langevin [40]. Many derivations were in fact proposed all along the 20th century issued from both Special relativity and General relativity [START_REF] Rizzi | Relativity in Rotating Frames[END_REF][START_REF] Tartaglia | [END_REF], even though the Sagnac effect is usually deemed to be a Special relativistic effect [START_REF] Malykin | Ring Interferometry[END_REF].

5. Remember that Einstein was not awarded the Nobel Prize for the relativity, but for the photoelectric effect in 1921 ! The Nobel Committee has never given this prestigious award to anybody, in recognition of the discovery of the relativity. The name of the eminent mathematician Henri Poincaré, the other important contributor to this theory, had been suggested (by a committee led by Georges Darboux in 1905 [START_REF] Archives | université de Lorraine)[END_REF]) ; but Poincaré did not receive this prestigious award either. If he had, would Albert Einstein have jointly received the award ? It is a remarkable fact that Poincaré presented the great trails of his work on the special relativity in the Comptes Rendus in 1905 [START_REF] Poincaré | [END_REF]. In the same year, a lesser known physicist, Albert Einstein, independently presented in the Annalen der Physik, similar results, but more complete [35]. It should be noted that Poincaré, too conservative, never totally got rid of the aether, whereas Einstein considered this concept superfluous.

6. Also let us note that Paul Langevin is recognized as one of the first physicists to have defended the relativistic conceptions in France in numerous seminars ; in particular the contributions of Einstein in the field, better known as they had been eclipsed in this country by the work of Poincaré on the same topic.

7. The quotations were very scarce between this date and the 1960s where the effect became well known for its applications in gyrolasers (first built in 1963). The first significant paper in english relating to the Sagnac effect is the study of Post published in 1967, which contained also some historical information [19].

The reasoning, both simple and elegant, of Langevin can be shown in a few lines. This illustrious physicist begins by noting that this is an experiment of the first order in v c (more precisely in Rω c , if R denotes the radius of the disk and ω the angular velocity, assumed to be constant), contrarily to the experiment of Michelson-Morley which involves a second-order effect. We start from the Minkowskian metric associated with the central observer, denoted by O, located at the center of the rotating disk. In his article Langevin used the cartesian coordinates, we modernize the presentation a little by choosing the cylindrical coordinates, which are better suited to the analysis of the problem (by eliminating the third coordinate z which is superfluous). We write the minkowskian ds 2 , expressed in these coordinates

ds 2 = c 2 dt 2 -dr 2 -r 2 dθ 2 (5)
The observer O then performs the global transformation of the coordinates

(θ ∈ [0, 2π]) (t, r, θ) -→ (t, r, θ + ωt) (6)
imposing the constraint r < R < c ω . This transformation means that the observer O (inertial) now uses a coordinate system that accompanies the disk in its rotation. From the inertial observer O point of view, the source, the detector and the non-inertial observer D (assumed to be at the same place) are located on a radius with a fixed orientation, for instance θ = 0. Neglecting a small second-order term in Rω c , the ds 2 (relativistic invariant) takes the form

ds 2 = c 2 dt 2 -dr 2 -r 2 dθ 2 -2r 2 ωdθdt (7)
Following Langevin, this expression is called the metric of the rotating disk.

It will be noted the presence of a cross term in dθdt, source of the impossibility of synchronizing clocks, uniformly distributed around the periphery of the disk and connected thereto (a relativistic phenomenon which is the essence of the Sagnac effect). The path in space-time of the light rays is determined by the cancellation of the ds 2 and by putting dr = 0 and r = R (the spatial trajectory of the light rays is circular and follows the periphery of the disk). We obtain an equation of the second degree in dt. The solutions (Rω c) are written

dt ± = R 2 dθ c 2 ω ± c 2 R 2 + ω 2 R 2 dθ c 2 ω ± c R (8)
It is then integrated on the periphery of the disk (noting that dt ± > 0, that is dθ > 0 for dt + and dθ < 0 for dt -). It yields

∆t = t + -t -= 4πR 2 ω c 2 (9)
We find again the relationship given by eq. 4.

Figure 4 illustrates this situation. When the disk is in rotation, the light beams are represented in thick gray lines, labeled respectively + for the beams traveling in the trigonometric direction, -for the ray moving in the opposite direction (in dashed lines are the beams for the disk at rest). For the fixed observer L, the difference in velocities between the light beam + (resp. -) and the observer D is c -v (resp. c + v). If the speed of the light beams with respect to this observer is c, the same value is definitely measured by the observer D linked to the disk, whatever the direction of travel of the considered beam. The points D + and D -shown in figure 4 are identified 8 .

Figure 4 8. But then we may wonder why, since the waves have the same speed in a direction and in the other for the observer D, is there a difference in the travel times measured by this observer. The subtlety is that the points D+ and D-reported in figure 4 are actually distinct and are identified only by abuse. A close analogy can be made with the screw dislocations seen in material sciences. A more mathematical approach can be envisioned with help of the complex function Lnz. This function is continuous and differentiable throughout C, but a cut must be created along any line joining 0 and infinity. The cut is made necessary because when the z phase jumps to ±2π, the Lnz function does not return to its initial value. From a topological (and not merely geometrical) point of view, we must therefore see the two regions represented in Figure 4, ] -2πR, 0[× t and [0, 2πR[× t as a two-leaf foliation of the space-time of the rotating disk.

The universality of the Sagnac effect

The Sagnac effect appears universal. More specifically speaking, one can cite the experiment carried out by Hafele and Keating in 1971 [45]. It consisted of the following. Two planes each equipped with an atomic clock flew around the Earth, one in the direction of the Earth's rotation and the other in the opposite direction (the two clocks are synchronized at the start). On arrival, it is observed that the two clocks indicate different times (in fact it is simply a revisited Sagnac experiment). A practical conclusion is that one cannot synchronize a collection of clocks distributed all along the Earth equator. The effect in question is taken into account in the GPS calculations. The phase shift of the Sagnac experiment was also measured in the case of coherent beams of relativistic neutrons travelling in opposite directions along the cicumference (closed path) of a rotating disk [START_REF] Hafele | [END_REF]45]. More generally, since the first atom interferometry experiment in 1991, measurements of rotation using the Sagnac effect have been currently performed with atoms [46]. Finally the quantum mechanical phenomenon predicted by Y. Aharonov and D. Bohm [47] can also be perceived as some kind of Sagnac effect (by replacing the light beams by charged particle beams (charge q) in the experiment, formally substituting the rotation ω by a magnetic field B and doing 8π λc -→ q in the relationship (3)) [48].

Conclusion

In Sagnac's view, his experiment had been initially designed to prove the existence of the aether, the assumed medium of propagation of light waves. It is now known that this hypothetical medium, with its clearly very contradictory physical properties, does not exist. In 1921 Langevin gives the appropriate (relativistic) explanation in a very short note to the Comptes Rendus. The Sagnac effect is today at the basis of virtually all of our terrestrial and satellite inertial guidance systems. These technological developments were far from being imaginable when this famous experiment was conducted.

Figure 1 :

 1 Figure 1 : Sagnac interferometer mounted on a turntable (reproduced from his paper of 1913)

Acknowledgments

We are indebted to Professor J. Villain for carefully reading the manuscript and providing us with detailed comments and suggestions.


G Pascoli
email: pascoli@u-picardie.fr

An elementary pedagogical derivation of a lot of metrics, seen in the very large context of metrics theories, is supplied, starting from very basic algebra (especially without Christoffel symbols and tensors) and a few formulas issued from Special relativity.

Introduction

More especially some of these metrics are usually studied in the framework of General relativity. However, in a rigorous context, we know that the resolution of the Einstein equations of General relativity needs a rather evolved background in tensorial calculus (generallly seen at the graduate level). On the contrary, the aim is here to develop in undergraduate students the sense of the heuristics, starting from their knowledge acquired at their proper level, allowing them to very rapidly reach the same results but without cumbersome calculations. In a more extended context this type of method, even though not optimal, is very often sufficient to speed up the process of finding a solution.

The static black hole (Schwarzchild metric)

We begin with the metric of the external region (assumed to be empty) surrounding a static body, with spherical symmetry and of masse M . The ds 2 is given by the well-known expression (Schwarzchild metric) [START_REF] Misner | Gravitation[END_REF]

ds 2 = 1 - r S r c 2 dt 2 - dr 2 1 -r S r -r 2 (dθ 2 + sin 2 θ dϕ 2 ) ( 1
)
where the spherical coordinates has been used (r, θ, ϕ) and the time coordinate is denoted by t. These coordinates are linked to the observer located at the infinity. The Schwarzchild metric is thus just the point of view of this special observer who studies the gravitational configuration surrounding the mass M . We have put as usually r S = 2GM c 2 . It is one of the very few exact solutions of the General relativity, a highly nonlinear theory of physics. A heuristic method makes it possible to find it very quickly 1 . An observer who evolves in free fall does not feel the gravity. For him the metric is minkowskian. Certainly this is strictly true only when the gravitational field is uniform, which is far from the situation encountered here, except if we consider an infinitesimal region. We attach to this observer a meter-stick represented by the infinitesimal space (radial) coordinate dρ and a clock represented by the infinitesimal time interval dτ (leaving aside the angular variables θ and ϕ). We have

ds 2 = c 2 dτ 2 -dρ 2 (2)
The elements (dτ , dρ) and (dt, dr) are linked by the relations defining the proper time and proper length, let

dt = dτ 1 -v 2 c 2 dr = dρ 1 - v 2 c 2 (3)
These expressions are directly issued from Special relativity. In Newtonian mechanics, the free falling velocity of any particle toward a point attractor, of mass M , and located at the origin of coordinates r = 0,

is v = -2GM
r (assuming that the particle in question starts from infinity with zero initial velocity). An intuitive idea is to abrutly take this value for v in the equations (3) 2 . Substituting (3) in (2), one immediately finds the relation (1) (assuming that the transverse lengths to the movement, orthoradial and azimuthal, are not altered). Even if such a procedure could appear unsound, the aim is to obtain a trial metric in a very simple way. This metric must then be checked using a Mathematica-type calculator, which can easily be manipulated even by undergraduate students.

The rotating black hole (Kerr metric)

The static metric that we just described cannot be an accurate representation of the outer region at the ultimate stage of a massive star at the end of evolution. Indeed, a star is always rotating. When it collapses, its kinetic moment is conserved and its rotation therefore correlatively increases. A more realistic model must therefore introduce, in addition to the unique parameter of a static spherical object, that is its mass M , a second parameter representing the angular momentum (it is usually denoted by J). The metric is no longer Schwarzchild's. A pure case concerns the exact metric attached to a rotating (not charged) black hole. As we know this object is defined by only the two parameters M anf J. It has long been sought (if such an exact metric existed were unknwon). It was not until 1963 (about fifty years after the discovery of the Schwarzchild metric) that an exact solution was provided by Roy Kerr [START_REF] Kerr | [END_REF]. A reasoning on symmetry makes it possible to greatly simplify the search. First, we can assume that the rotation is around a fixed axis, for example the Oz axis. Since the metric revolves around this axis, its coefficients must be independent of the azimuth angle ϕ. If the motion is stationary, these coefficients must not depend on time either. The g ij are such that g ij ≡ g ij (r, θ). There are four coefficients which define the Schwarzchild metric g tt and g rr , g θθ , g ϕϕ . How many are needed to define, in all generality, the outer metric of a body assumed to rotate in a stationary manner ? A track is suggested to us by the following heuristic reasoning 1 . It is assumed that near the rotating object the coordinate network (r, θ, ϕ) participates at the rotation (Fig. 1). The effect is called frame -dragging. Obviously the effect in question will be more and more important as r decreases.

Figure 1 -The Lense-Thirring approximation First, we place a local frame-dragged observer O l located at radial distance r in the outer region of the rotating object whose center is in r = 0. This object is not necessarily a black hole at this stage. We can think to the outer region of a star or even of a planet like the Earth. In this case r r S . For this observer, the azimuthal velocity measured for a moving particle is rsinθ dϕ dt -rsinθω 2 . It is therefore advisable to make in the ds 2 the direct substitution

rsinθdϕ -→ rsinθdϕ -rsinθωdt (4)
At this level we must give the expression of ω 3 .

ω = ω S ( r S r ) 3 (5)
where ω S is an arbitrary constant. By putting J = M r 2 S ω S (this quantity is just the angular momentum of a homogeneous thin ring of radius r S and mass M put in rotation around its normal axis), the latter expression can still be written

ω = r S J M r -3 (ω = 2GJ c 2 r 3 ) (6)
Now substituting the expression (4) in equation ( 1), we obtain for the ds 2 , in the weak field approximation and at the first order in J (we put α = J M c )

ds 2 = 1 - rS r c 2 dt 2 -(1 + rS r )dr 2 -r 2 (dθ 2 + sin 2 θ dϕ 2 ) + 2rSα r sin 2 θdϕcdt (7)
where the non-zero coefficients are g tt , g rr , g θθ , g ϕϕ and the non-diagonal coefficient g ϕt = g tϕ , which are five in all. It is legitimate to think that it is still the minimum number of non-zero coefficients needed to represent the metric of a rotating black hole. This expression was obtained by J. Lense and H. Thirring from 1918 (but in a rigorous way, using the formalism of the General relativity [3]). It is correct at first order in J. The phenomenon of frame-dragging has been proven experimentally by the exploitation of data from the LAGEOS mission [4] and, more specifically, by the Gravity Probe B satellite [5].

-The Kerr metric

If the (inertial) coordinate network is frame-dragged, it seems natural to assume that it is also flattened by a centrifugal effect. The coordinatespheres (r = Const) then become flattened confocal spheroids, each of them being still labeled with r = Const, see Fig. 2. The coordinates r, θ, ϕ from the current point P , plotted on the spheroid drawn in figure 3, are actually those of the point P , projection of P on the inscribed sphere. These deported spherical coordinates are the Boyer-Lindquist coordinates [6].

N o rotation Rotation x = r sinθ cosϕ x = √ r 2 + a 2 sinθ cosϕ y = r sinθ sinϕ -→ y = √ r 2 + a 2 sinθ sinϕ z = r cosθ z = r cosθ Figure 2 The minkowskian ds 2 (ds 2 = c 2 dt 2 -dx 2 -dy 2 -dz 2 ) is written in this case 4 ds 2 = ( r 2 Σ )(1 + a 2 r 2 )(cdt -asin 2 θdϕ) 2 (8)
- sin 2 θ Σ ((r 2 + a 2 )dϕ -acdt) 2 -( Σ r 2 ) dr 2 1 + a 2 r 2 -Σdθ 2
where we are put Σ = r 2 + a 2 cos 2 θ. Two elementary manipulations are now needed to modify the latter expression to make it compatible with the limiting cases (7) ( r S r 1 and α 2 negligible) and (1) (α = 0). To do this, we first identify α and a (the flattening of the coordinate-spheres is due to rotation), then we add the gravitational content (not present in ( 8)), that is we replace 1 + a 2 r 2 by 1 + a 2 r 2 -r S r . We eventually obtain

ds 2 = ∆ Σ (cdt -asin 2 θdϕ) 2 - sin 2 θ Σ ((r 2 + a 2 )dϕ -acdt) 2 - Σ ∆ dr 2 -Σdθ 2 (9)
putting ∆ = r 2 (1 -r S r + a 2 r 2 ). The quadratic form (9) is the Kerr metric associated to a rotating black hole (with the restrictions which are imposed : stationarity, axisymmetry and asymptotically flat metric). The intuition has guided us toward the accurate result, but of course we still have to validate this expression by showing that it checks Einstein's equations in vacuum (the ultimate referee !). Kerr's metrics is very often provided without any support, even at an advanced level, creating a dose of perplexity among students. The simple procedure developed above can help to dissipate this unrest.

The rotating disk of Sagnac

The sagnac experiment is well known [7,8]. The Sagnac effect can be derived from Special relativity without explicitly appealing to General relativity [9]. However the nature of the metric attached to a rotating disk has long been a controversial topic [START_REF] Rizzi | Relativity in Rotating Frames[END_REF]. A simple method of deduction was early proposed by Langevin [8]. Let an inertial observer O be located at r = 0 of a polar coordinate network (r, θ). We write the minkowskian ds 2 , expressed in these coordinates

ds 2 = c 2 dt 2 -dr 2 -r 2 dθ 2 (10)
The observer O wishes to accompany the rotating disk (radius R) in its motion. Then he performs the global transformation of the coordinates (θ ∈ [0, 2π])

(t, r, θ) -→ (t, r, θ + ωt) (11)
imposing the constraint r < R < c ω . This transformation means that the observer O (always inertial !) now uses a coordinate system that accompanies the disk in its rotation. The ds 2 (relativistic invariant) takes the form

ds 2 = c 2 1 - r 2 ω 2 c 2 dt 2 -dr 2 -r 2 dθ 2 -2r 2 ωdθdt (12)
By rearranging of the terms, we can still rewrite the latter expression into the form

ds 2 = c 2 1 - r 2 ω 2 c 2 (dt - r 2 ω c 2 dθ 1 -r 2 ω 2 c 2 ) 2 -dr 2 - r 2 dθ 2 1 -r 2 ω 2 c 2 (13) Introducing a new interval of time dT = 1 - r 2 ω 2 c 2 (dt - r 2 ω c 2 dθ 1 -r 2 ω 2 c 2 ) (14)
We finally obtain

ds 2 = c 2 dT 2 -dr 2 - r 2 dθ 2 1 -r 2 ω 2 c 2 (15)
The purely spatial part, which can now be extracted from it, is

dσ 2 = dr 2 + r 2 dθ 2 1 -r 2 ω 2 c 2 (16) whose the curvature scalar is R = -6 ω 2 c 2 (1 -r 2 ω 2 c 2 ) -2
< 0. The spatial part is hyperbolic. However, the preceding reasoning is true if the rotational velocity ω is such as ω 2 c 2 R 2 . When the rotational velocities are much higher, the centrifugal forces cannot be neglected. Let us select the point of view of an observer, located at r and sliding without friction in a pipe radially directed and linked to the rotating disc. This observer has the impression to evolve in free fall, not toward the center as in the situation of paragraph 1 (free fall toward a point attractor), but on the contrary toward the outer region.

Resuming the heuristic reasoning of paragraph 1, but with v = rω (taken as radial velocity of free fall in the pipe), we directly deduce

ds 2 = c 2 (1 - r 2 ω 2 c 2 )dt 2 - dr 2 1 -r 2 ω 2 c 2 -r 2 dθ 2 (17)
with the restriction Rω < c where R = M aj(r) (in the case of the rotating disc the limiting value c ω for the radius plays the role of a Schwarzchild radius). However this result is just a part of the affair. Like what we did in paragraph 3 about the Kerr metric, we still have to make compatible (12) and (17). Following the heuristics, a very simple and intuitive manner consists to superimpose the two expressions ( 12) and (17) to create an emergent metric with the simultaneous properties of the two original ones (but that they did not possess separately). Eventually, we get

ds 2 = c 2 (1 - r 2 ω 2 c 2 )dt 2 - dr 2 1 -r 2 ω 2 c 2 -r 2 dθ 2 -2r 2 ωdθdt (18)
We can call this form the dematerialized rotating disk metric. The purely spatial part (doing dt = 0) is then

dσ 2 = dr 2 1 -r 2 ω 2 c 2 + r 2 dθ 2 (19)
whose the curvature scalar is constant and equal to R = 2ω 2 c 2 > 0. The spatial part is in this case elliptic (the geometry of the sphere). The latter geometry also shares some analogies with the well-known conformal Beltrami model, even though the latter one has a constant but negative curvature (the geometry of the pseudosphere or tractricoid where R -→ -R). Let us mention that the metric (18) is no longer minkowskian (contrarily to (12) that is always minkowskian in spite of its appearance !). Another important point to specify is that (18) cannot be deduced from the general relativity because the mass and gravitation are absent. Thus the procedure exposed above is worth for a demateralized disk without mass and gravitation, that is a rather academic (and unphysical) situation. In fact, a very high velocity ω brings a very large kinetic energy or equivalently a very large mass. In this case the situation drastically changes. The inertial frame-dragging must be taken into account. Another point is that the disk can no longer maintain an uniform rotation. This one becomes automatically differential. The metric (18) must be replaced by the Lense-Thiring one with a variable ω (7). Let us notice that the sign of the crossed term in dθdt is then reversed. This reverse of sign is easily explained owing to the fact that in the Sagnac experiment there is no frame-dragging for a dematerialized disc, whereas this effect is naturally existing in a massive disc. The coordinate systems are also different. For the metric (7) the inertial observer is located at infinity, while he is at the origin of coordinate system for the expression (18). In parallel the expression (7) is asymptotically flat at infinity and divergent at the origin, while (18) is flat at the origin and divergent at r = ω c .

conclusion

We have shown that the metric attached to a static black hole (Schwarzchild metric), that of the rotating black hole (Kerr metric) or that of the dematerialized rotating disk (Sagnac experiment) can be derived with very simple physics accompanied by some heuristic arguments. Such a treatment appears available for gently introducing General relativity to students at the undergraduate level.

1.

A heuristic step obliterates justification and rigor to go straight to the goal ! 2. An explanation is needed here. This idea is supported as follows. The corresponding lagrangian for radial infall is

L = c 2 (1 - r S r ) ṫ2 -(1 - r S r ) -1 ṙ2
where we have put ṫ = dt ds and ṙ = dr ds . We write now the Euler-Lagrange equation for r

d ds ( ∂L ∂ ṙ ) - ∂L ∂r = 0
or more explicitely

d 2 r ds 2 - 1 2 r s r 2 1 -r s r ( dr ds ) 2 + c 2 2 r s r 2 (1 - r s r )( dt ds ) 2 = 0 *
On the other hand, the Schwarzchild line element can still be written 1. This step is assuredly non-rigorous, even unsupported in any way. It is a simple scaffolding which is then removed once the adequate metric is built ! 2. We use here the classical formula for velocity composition : relative velocity = absolute velocity -velocity of the dragged frame (Galileo -Newton). The time remains Newtonian and cdt -→ cdt (we neglect any term expressed as a function of the square of the speed of rotation). Similarly, the product of the speed of the particle by the rotational speed is assumed to be negligible compared to c 2 ).

1 = c 2 1 - r s r ( dt ds ) 2 - 1 1 -r s

3.

The speed of rotation must decrease rapidly when we move away from the rotating object (just like a whirlwind in a river is limited to a small singular region). On the other hand, the exponent (the strong decrease in r -3 ) is more difficult to justify. Let's say that the extra term added in the metric must not be felt at infinity (which would be the case if the decrease was slower). A physical solution requires that the metric be asymptotically flat at infinity. 4. It is appropriate to express dx, dy, dz according to dr, dθ, dϕ. On the other hand, we have symmetrized on the pair of time-space variables (dt, dϕ).

r

  Substituting for dt ds (issued from * * ) in the equation * gives thus dr ds → 0 for r → ∞ and inward free fall ⇒ dr < 0). The Euler-Lagrange equation for t is An astronaut at rest, placed at a position labelled by the space coordinate r and at a time coordinate labelled by t, measures the free fall speed of a particle. At this stage, we assume then that the intuitively deduced ds 2 (1) is accurate. For this observer the proper length is then dρ * = ( 1 -r S r ) -1 dr and the proper time dτ * = ( 1 -r S r )dt ). The measured velocity is given byv (r) = dρ *By identifying, we find v (r) = -c r s r This relationship is formally equivalent to that obtained in newtonian mechanics even though by a luckily coincidence ! Section 3

NotesSection 2