Laurent Pierre

Bernard Anselmetti

Comparison of analysis line and polytopes methods to determine the result of a tolerance chain

Keywords: Functional tolerancing, tolerance analysis, Analysis line, Polytop, ISO tolerancing standard, 1. Functional tolerancing context

Functionnal tolerancing must ensure the assembly and the functioning of a mechanism. There are several types of functional requirements: mainly the assembly requirements which are expressed in terms of maximum material conditions and accuracy requirements that impose minimum material conditions. The approach consists of five steps: -Definition of each functional requirement by geometrical characteristics between some functional surfaces of mechanism and determination of the limit values of this characteristic allowing insuring the requirement with an acceptable risk, -Determination of the annotations of each part with the standards ISO of dimensioning and tolerancing, -Choice of the tolerances of each specification of each part, -Verification that the accumulation of the tolerances guarantees the respect for all the requirements, -Optimization of the allocation of tolerances to minimize cost or to maximize the capability of manufacturing processes. In this approach of functional dimensioning, the tolerance analysis is a very delicate stage to take into account the dimensional and angular effects. The problem is to ensure that the studied functional requirement will be well respected considering the tolerances on influential parts. The aim of this paper is to compare two methods of tolerance analysis of a mechanical system: the method of "analysis lines" and the method of "polytopes". The first method needs a discretization of the ending functional surface according to various analysis lines placed on the outer-bound of the face and oriented along the normal of the surface. For each analysis line, the result is the sum of the influences of each junction defects on the search characteristic using pre-established relations for classical types of junctions. The calculation is very fast and it can express the result as a formula according to the tolerances, if necessary with a statistical approach. The second method uses polytopes (finite sets of point in n ). The polytopes are defined from the acceptable limits of the geometric deviations of parts and possible displacements between two parts. Minkowski sums and intersections polytopes are then carried out to take into account all geometric variations of a mechanism. Tools like OpenCASCADE have been developed to quickly realize intersections and Minkowski sums. This comparison is performed on a simple example specified with the method CLIC. The requirement is the assembly condition of a shaft in two bores of two different parts which involves straightness in common zone of the bores. Comparison of these results shows that solutions are identical. Assumptions and model behavior bonds are also identical.

of parts at lower cost, it is very important to allocate the widest possible tolerances, while ensuring that the functional requirements of the mechanism are respected. With this approach, the designer chooses the geometric specifications and tolerances. He must ensure that the combination of the variations permitted by the tolerances will be compatible with each requirement. This step is called tolerance analysis.

Traditionally, tolerance analysis is provided by a tolerance chart model in only one direction which is simply to sum of tolerances of influential parts. Soon as there is an angular effect in the chain of dimensions, designers are still quite poor. Indeed, there are very few commercial software tools for the analysis of tolerances in 3D. The principle is based on the simulation of assemblies:

-Mécamaster [START_REF] Clozel | MECAmaster: a Tool for Assembly Simulation from Early Design, Industrial Approach, dans: Geometric tolerancing of products[END_REF] simulates the assembly of nominal parts but introduces deviations between parts corresponding to tolerances of the support surfaces or clearance in joint.

-Cetol [2], 3DCS [START_REF]DCS -Variation Analysis and Tolerance Analysis[END_REF] generate, by Monte Carlo methods, parts with position, orientation or form defects corresponding to the ISO specifications. Such software measures the characteristic of the requirement on the set of assemblies obtained with simulated parts.

-Anatole [START_REF] Falgarone | An innovative Design Method and Tool for Structural and Functional Analysis[END_REF] characterizes defects of surfaces with small displacement torsor. The equations expressing the displacement of the parts relative to their nominal positions are based on contact between parts to give a formal relationship of part tolerances.

Since 1997, B. Anselmetti developed the CLIC method (a French acronym for "Localization Tolerancing with Contact Influence") [5]. This system proposes the automatic generation of functional requirements and functional specifications in accordance with ISO tolerancing standards. The result of each tolerance chart is described by an EXCEL formula according to tolerances of influential parts. This structure allows optimizing the nominal dimensions of CAD model to maximize tolerances [START_REF] Anselmetti | Part optimization and tolerances synthesis[END_REF].

In this software, tolerance analysis is ensured by the analysis line method. The ending functional surface is discretized by various analysis lines placed on the outer-bound of the face and oriented along the normal to the surface. The principle consists in the sum of the influence of defects in each junction on analysis line using relationship pre-established for all classical types of junctions. The calculation is very fast and expresses the result with a formula function of tolerances with worst case or with a statistical approach.

The method of operations on polytopes in tolerance analysis is a variational approach based on the work initiated by Flemming in 1988. It is based on the operations of sets of geometric constraints [START_REF] Fleming | Geometric relationships between toleranced features[END_REF]. The set of possible positions of a surface within a tolerance zone is characterized by a set of geometric constraints. This allows characterizing the geometric variations of a part according to orientation or position ISO specifications [START_REF]Geometrical Product Specifications (GPS), Geometrical tolerancing, Tolerances of form, orientation[END_REF], [START_REF] Giordano | Mathematical representation of tolerance zones[END_REF], [START_REF] Roy | Representation and interpretation of geometric tolerances for polyhedral objects. II.: Size, orientation and position tolerances[END_REF], [START_REF] Davidson | A new mathematical model for geometric tolerances as applied to round faces[END_REF]. Similarly, a set of geometric constraints characterises all relative positions between two distinct surfaces potentially in contact [START_REF] Fleming | Analysis of Uncertainties and Geometric Tolerances in Assemblies of Parts[END_REF]. Fleming established the correlation between the accumulated defects limits on parts and the Minkowski sum of finite sets of geometric constraints [START_REF] Fleming | Geometric relationships between toleranced features[END_REF]. A summary of these issues is detailed in [START_REF] Srinivasan | Role of Sweeps in Tolerancing Semantics[END_REF]. Giordano showed that the relative position between two parts in contact resulting from a number of potential contacts can be formalized by an intersection operation of a set of geometric constraints [START_REF] Giordano | Clearance Space and Deviation Space[END_REF]. More generally, the variational tolerancing approach is to characterise the relative position between two surfaces of two parts of any system by intersections and Minkowski sums of sets of geometric constraints. This sets of constraints result from ISO specifications expressed in parts and contact characteristics between two parts [START_REF] Teissandier | Operations on polytopes: application to tolerance analysis[END_REF]. Minkowski sum algorithms applied to the problem of tolerance analysis have been developed [START_REF] Wu | Improvements to algorithms for computing the Minkowski sum of 3-polytopes[END_REF], [START_REF] Teissandier | Algorithm to calculate the Minkowski sums of 3-polytopes based on normal fans[END_REF]. A variational method of tolerance analysis using the method based on operations on polytopes is proposed by [START_REF] Pierre | Integration of thermomechanical strains into tolerancing analysis[END_REF] and [START_REF] Pierre | Variational tolerancing analysis taking thermomechanical strains into account: Application to a high pressure turbine[END_REF]. This method was developed to take into account thermomechanical strains.

The aim of this paper is to compare these two methods on an example to show the similarities and limitations.

Definition of requirement

Mechanism description

The mechanism is composed of two rigid parts with two bores (Fig. 1). The shaft is represented by a single cylinder which is mounted with clearance in a bore of the cover and in a bore of the housing. The cover is assembled and screwed onto the housing, before introducing the shaft. It must therefore be ensured that the shaft can be fitted in the mechanism in any case.

Functional tolerancing

Functional tolerancing was performed with the Quick GPS software developed by B. Anselmetti [START_REF] Anselmetti | Quick GPS:A new CAT system for single-part tolerancing[END_REF] in the CATIA environment with Functional Tolerancing Annotation Workshop.

The junction is described by a positioning table of the cover on the housing (Fig. 2).

POSITIONING The housing is the base of this mechanism. The cover is setting-up on the housing by a primary plane and two bores which receive two pins blocked in the housing. There is clearance between the pins and the cover.

A simple matrix shows all the links in the part studied (Fig. 3). The hole of the cover has to be positioned relative to the junction with the housing. References Fig. 3. Link inside cover These data are sufficient for the QUICK GPS software that automatically generates annotations in the 3D CAD model (Fig. 4). The designer can adjust the tolerances. Each part has a main datum reference frame AB and a location of the bore.

Functional requirement

This paper examines more specifically the assembling of the shaft in the sub system composed by housing and cover. The assembling is difficult when the shaft and the bores are at maximum material condition, when the cover and housing are shifted due to the clearance between the pins and the cover (Fig. 5).

The shift is maximum when the distances between the pins are the same on both parts and the pins and the holes in the cover are at minimum material condition with a maximum clearance of 0.04.

The problem is really three-dimensional, because the holes are not in the plane containing the two pins. But to simplify the Fig. 5, the pins are shown in the plane of bores.

Methodology

At maximum material condition, the shaft is considered a perfect cylinder of diameter 29.96. The envelop diameters of the bores are 29.98 mm. If the straightness in common zone of the two holes is r, then the free space for the shaft is 29.98 -r. The straightness has to be less than 0.02 in all cases.

Tolerancing and significant dimensions

0 z t4 M t3 M t2 M A B t1 M A P A B A L y Ø30 0.

Analysis line method

Principle

The analysis lines method was developed in 2004, to calculate the resulting 3D with a statistical approach. The challenge is to globalize the influence dues to translation and angular deviations, whatever the tolerance values. Calculation equations show that the accumulation occurs naturally in certain points and in particular directions called analysis points and analysis directions.

The analysis line method is based on transfer relations that have been established for ten classical junctions (Primary plane / secondary plane / tertiary plane, Primary plane / secondary cylinder / tertiary, Primary plane / 2 secondary holes, Primary plane / n secondary holes; Primary cylinder / secondary plane / tertiary, straightness or flatness in common zone ...).

The linear relationship gives the influence of a given junction in the worst case, directly on the functional requirement. This requires a discretization and successive studies in a limited number of directions. The advantage of this method is that for each study, the analysis direction is known. It is very easy to define the position of the part in the worst case permitted by clearance in the junction.

Other junctions can often be modeled by 6 points of contact with a linear relationship. Some complex connections cannot be treated directly. The solution proposed by Robin Chavanne [START_REF] Chavanne | Phd thesis : Contribution au tolérancement fonctionnel 3D des mécanismes complexes : Synthèse des spécifications et Analyse de tolérances[END_REF] is to use a solver to find the worst case situation and to find the 6 contact points to establish the linear relationship.

Transfer of straightness

The requirement studied thus reduces to straightness between two coaxial bores. The deviation from straightness r is the diameter of a cylinder containing both real axis of the two bores.

The problem is three-dimensional. The analysis line method proposes to study in 8 radial directions 1 f to 8 f . Fig. 8 can be applied to all cases of straightness in common zone of two coaxial cylinders (Fig. 8a).

A location and an orientation specifications are placed on each bore with respect to a common datum reference frame fixed on the housing (Fig. 8b).

Fig. 8c shows that in each planes in fi direction, straightness is proportional to the sum of the maximum displacement   The relationship of transfer is therefore:

        , , ; max ,, C HC C BC i H HC H BC E d C d C EL r E d B d B EL                 fi fi fi fi (1) These displacements   , H dB fi ,   , H dC fi ,   , C dB fi ;   , C
dC fi should be calculated relative to the same reference. The reference is chosen on the median line of both axis B of the datum reference frame AB of the housing. Firstly, the maximal displacement of the points C B and C C of the cover axis are calculated with respect to the nominal cover axis according to the specifications of the bore.

The nominal cover axis is also the median line of both reference axis B of the cover In the worst case, point C B moves in C B . The displacement of point C B is limited by the location tolerance 3 t :

  3 t ,=2 C dB fi (2) Point C C moves in C C .
The inclination of the axis of the bore is limited by the specification of orientation 4 t . The displacement of point C C is:

    3 4 OB OC C t 2 t +L E , L C dC   fi (3)
The calculation is similar to the points H B and H C of the housing axis with the tolerance 1 t in location and 2 t in orientation.

    1 1 2 BC H t2 , , t 2 t L E H H dC dB    fi fi (4)
By symmetry, these relationships are independent of analysis directions fi . Displacements of points B and C the nominal cover axis must now be determined according to the housing.

Influence of clearance in junction housing/cover

The problem is to calculate the displacements of the nominal axis of cover in C B and in C C The cover is setting-up on the housing with a primary plane and two pins. The datum reference frame of the housing is AB. Datum plane A of the housing and datum plane A of the cover are supposedly confused.

The nominal axis of the holes being perpendicular to the primary plane, there is no angular deviations between both nominal axis of the housing and cover.

The shift of the nominal axis is due only to the mobility permitted by the clearance J between the housing and the cover. Fig. 10 shows both axis E and F of the housing holes, and both axis E' and F' of the cover holes in the cover.

The connection between these bores is provided by the pins blocked in the housing. Displacement is maximum when the distance between the holes are identical in housing and the cover and when the pins and the cover holes are at least material condition with a maximum clearance J.

If the analysis line  

, C dB fi cut segment EF , the maximum displacement in the direction fi is equal to 2 J .

If the analysis line  

, C dB fi does not cut the segment EF , the displacement of C B is maximal when E move in E' and when F move in F' in  direction. The maximal clearance between pins and hole of cover is:

The displacement of the point

max_ hole min _pin 5 6 - J D D t t    (6)
Both distances EF L =100 mm and OB L = 60 mm, 1 cos 1.56

 

. Calculated displacements are: For  f1 x and  f5 x :

  ,2 C d B J  fi (7)
For f3 and f7 , cos 1

  :     cos ,2 .56 2 1 C d B J J      fi (8) For f2 , f4 , f6 et f8 , 4   :     co ,2 s 1.10 22 4 C d B J J       fi (9)

Full condition

The full condition giving the straightness in the worst case is the maximum value of 8 relationships in 8 directions fi .



  arctan C OE OB  
These relationships are linear and simple to integrate into a system of equations for the synthesis of tolerances.

Method by operations on polytopes

Transfer of geometric deviations on parts

A real surface resulting from the manufacturing is modelled by a surface of substitution that has the same type as the nominal surface. The deviations between the surface of substitution and nominal surface are called geometric deviations. The surface of substitution i of the part j is noted , ij . Geometric deviations between a surface of substitution 1,1 and a surface of substitution 1, 2 are formalized by a small displacement torsor 1,1/1,2 d   , where the vector 1,1/1,2 ρ is is the rotation vector of 1,1 with respect to 1, 2 while the vector B-1,1/1,2 ε is the translational vector of 1,1 with respect to 1, 2 at point B:

B 1,1/1,2 1,1/1,2 B-1,1/1,2 d       ρ ε (14)
According to the specification of the location of axis 1,1 of the bore of the housing 1 (Fig. 6), the axis of the cylindrical surface of substitution 1,1 is within a tolerance zone ZT . ZT is a cylinder of diameter 1 t  orthogonal to the plane A and centered on nominal axis relating to AB datum reference frame. To ensure that the axis 1,1 of the surface is located in the zone of tolerance, it is necessary to write, to both points C and D, the following relationship where i  n is a unitary vector orthogonal to the axis z and i  is the discretization angle around the axis z :

 

C AB i D AB i i i i i tt tt i i n i n n                                   εn εn n
x y [START_REF] Teissandier | Operations on polytopes: application to tolerance analysis[END_REF] The relation ( 15) is expressed according to the rotation vector and the translation vector of the deviation of the point P middle of the bounded by line segment B and C (Fig. 5):

      11 P-1,1/ 1,
AB AB i AB AB i i i i i tt tt i i n i n n                                       ε CP ρ n ε CD ρ n n x y (16)
The relations ( 16) define a bounded intersection of finite number of closed half-spaces whose boundaries are hyperplanes of 4 [15] intersection. In general, it can formalize the h-representation of a polytope called geometric polytope, noted g AB . This is 4-polytope in which two graphical representations are shown in Fig. 11. The Fig. 11a and the Fig. 11b represent two projections of 4-polytope respectively in  

, x P y  ,
and in  

, y P x  , . , P y  x  1 2 t  1 h t E  1 2 h BC h E L t E    1 2 t  1 2 h BC h E L t E    1 h t E  , Px  y  1 h t E  1 2 t  1 h t E  Fig. 11. Geometric polytope g 1,1/ AB D
The orthogonality of the housing axis 1,1 with respect to the datum A can be characterized by a geometric polytope g 1,1/ A D . The axis of the cylindrical surface of substitution 1,1 is within a tolerance zone TZ . TZ is a cylinder of diameter 2 t  orthogonal to the plane A . To ensure that the axis of the surface 1,1 , is within the orientation tolerance zone TZ , the distance to the axis y between the poins C and D must be less than the dimension of the tolerance zone:

    2 -1,1/ -1,1/ 2 cos sin with
, 0 and ,

C A D A i i i i i tt i i n i n n                        ε ε n n x y (17)
The relation ( 17) is expressed according to the rotation vector and the translation vector of the deviation of the point P middle of the bounded by line segment B and C:

    2 -1,1/ -1,1/ 2 cos sin with
, 0 and ,

C A D A i i i i i tt i i n i n n                        ε ε n n x y (18)
Relations ( 18) formalize the geometric polytope, noted g 1,1/ A D , corresponding to the location of 1,1 with respect to the reference A. This is a 2-polytope whose graphical representation is given in Fig. 12 in  

xy  , . y  x  2 h t E  2 h t E  Fig. 12. Geometric polytope g 1,1/ A D
Both specifications expressed on the axis 1,1 of the housing 1 must be respected simultaneously. The polytope characterizing the position of the axis 1,1 with respect to the datum reference frame AB is the polytope resulting from the intersection between the polytopes g ,

Px  y  2 h t E  1 2 2 2 BC h L t t E    2 h t E  , P y  x  2 h t E  1 2 2 2 BC h L t t E    2 h t E 

Transfer of contact deviations between the parts

Previously geometric deviations of parts were formalised by polytopes. In the following, deviations due to mobility and clearance in the joints will be formalised by polytopes. The joint cover/housing is composed of three elementary joints: a plane pair contact and two ball and cylinder pair contacts. To model the full joint, each of three joints is modelled by a contact polytope.

The planar pair contact between the plane 1,2 cover 2 and the plane 2,2 the housing 2 is without clearance (Fig. 6). The relation [START_REF] Pierre | Variational tolerancing analysis taking thermomechanical strains into account: Application to a high pressure turbine[END_REF] characterised the contact without clearance on the whole contour of the contact surface:

P-1,2/ 2,2 1,2/ 2,2 1,2/ 2,2 0 0 0         εz ρx ρy (19)
Relations ( 20) characterize the limits of displacements between surfaces 1,3 and 2,3 and between surfaces 1,4 and 2,4 of the both ball and cylinder pair contacts (Fig. 6):

  E-1,3/2,3 F-1,
i i i i i i JJ JJ i i n i n n                                   εn εn n
x y [START_REF] Anselmetti | Quick GPS:A new CAT system for single-part tolerancing[END_REF] The maximum clearance corresponds to cases where the diameters of the shafts are smaller and diameters of holes are larger:

max_ hole min _pin 5 6 - J D D t t    (21)
Relation ( 20) is expressed according to the rotation vector and the translation vector of the deviation of the point P:

      1,3/ 2,
i i i i i i JJ JJ i i n i n n                                     EP ρn FP ρn n x y (22)
The polytope of contact characterising the relative position between the housing 1 and the cover 2 results of the intersection between the polytopes defined in [START_REF] Pierre | Variational tolerancing analysis taking thermomechanical strains into account: Application to a high pressure turbine[END_REF] Geometrical deviations in parts and in contacts were respectively determined in sections 4.2 and 4.3. The relative position between the bore of the housing and the cover will be characterized by the calculated polytope defined by the following Minkoski sum:

g c c 1,1/ 2,1 1,1/ / / 2,1 AB AB CD CD  D D + D + D (23)
The functional polytope f 1,1/2,1 D characterizing the straightness between the bore axis of the cover and the housing is defined by the following relation:

f f f 1,1/ 2,1 1,1/ ZT / 2,1 ZT  D D + D (24)
The geometric polytopes f 1,1/ZT D and f 2,1/ZT D characterizing the position of the axis 1,1 and 2,1 with respect to the tolerance zone ZT are defined by the similar manner to the location of the axis 1,1 (see relation (16)). The Fig. 16 illustrates the result of the sum of Minskoski (24) characterizing the polytope calculated

 x  f 1,1/2,1 D 1,1/2,1 D 1 h  2 h  1 V Fig. 16. Inclusion of the calculated polytope 1,1/2,1 D in the functional polytope f 1,1/2,1

D

If the calculated polytope is included in the functional polytope, the straightness of the two bores is inferior to the value defined in the polytope functional. To determinate the value of straightness of the two bores is necessary to minimize the size of the functional polytope. When the size of the functional polytope is minimum, each vertex of the calculated polytope must be include into each half-spaces of the functional polytope. The straightness is defined by these inequalities:

        3
E L a E L L E E L r b E E L L L E r c E L L E L E r E E L L r                                                        () d      (25)
In these four inequalities (25), r is greather into the relation (a) than into the relation (b) and r is greather into the relation (c) than into the relation (d). With this method, the straightness value is the same as the previous method (see [START_REF] Srinivasan | Role of Sweeps in Tolerancing Semantics[END_REF]

Conclusion

Both methods integrate the influence of the geometrical defects of location and orientation, and the mobility allowed by the clearance in the joints. But they consider directly the situation in the worst case defined by the minimum material condition in the joints and the maximum material condition in the bearings. For primary planar joint, the analysis line method considers that the datum reference of one part remains in the tolerance zone of the support part. The polytope method neglects the flatness of planes in contact for superimposing the datum reference and the support surface. These two ways of expressing the assumptions lead to the same behavior. This study shows that geometric models are identical, and that the analysis lines means studying the polytope in particular directions. The polytope is more complex, but gives a global view of all degrees of freedom. The method of analysis line is more direct, but requires considering each requirement in several directions by discretization of the ending surface. Analysis line method directly gives the result as a linear formula of tolerances of influential parts, which subsequently allows to globally optimize the whole system of inequalities to maximize tolerances. Both methods require a lot of rigor and expertise. A specific tool in CAD system is necessary. The setting-up of parts and functional requirements are described by the designer. The tolerance analysis application collects geometric data and the influential annotations. The calculations can be fully automated to display the result in the form of the numerical value of the result to be compared with the limit value of the requirement.

  Fig. 1. Mechanism

Fig. 4 .

 4 Fig. 4. ISO tolerancing of the two parts

Fig. 5 .

 5 Fig. 5. Assembling requirement of the shaft

Fig. 6 .

 6 Fig. 6. Main specifications

Fig. 7 .

 7 Fig. 7. Discretization in 8 analysis directions

  of point H C of the housing axis in fi direction and the maximal displacement   , C dC fi of point C C of the cover axis projected in point C in fi direction. Lengths of bores are respectively C E and H E . The distance between the bores is BC L . The symmetrical condition is presented in Fig. 8d.

Fig. 8 .

 8 Fig. 8. Straightness between two different parts

Fig. 9 .

 9 Fig. 9. Influence of the cover

CBFig. 10 .

 10 Fig. 10. Influence of junction housing/cover

D

  represents the geometric polytope of the location of the surface 1,1 with respect to

D

  . This is a 4-polytope in which two graphical representations are shown in Fig.13. Fig.13aand Fig.13brepresent the projections of the 4-polytope in  

Fig. 13 .

 13 Fig. 13. Global geometric polytope g 1,1/ AB D In a similar manner to the location and orthogonality specifications of the axis 1,1 of the housing 1, the geometric polytope characterizing the position of the axis 2,1 with respect to the datum reference frame CD result from the intersection of the polytopes g 2,1/CD D and g 2,1/C D . This is 4-polytope in which two graphical representations are shown in Fig. 14. Fig. 14a and Fig. 14b represent the projections of 4-polytope in   , x P y

Fig. 14

 14 Fig. 14. Geometric polytope g 2,1/CD D

TABLE Parts

 Parts

	name: housing
		State
	Main reference frame
	PRIMARY FEATURE	SECONDARY FEATURE	TERTIARY FEATURE
	Faces Selection	Faces Selection	Faces Selection
	Simple plane	Set of parallel
		cylinders
	Show feature	Show feature	Show feature
	CONTACT	FITTING PINS
	No projection	Lenght of pin	No projection
	Special datum	Special datum	Special datum
	Breaking feature with:
		Variable distance
	Delete junction	Secondary	Tertiary
	Specif datum in frame
	Delete
			EXIT
	Fig. 2. Setting-

up table of the cover on the housing

  and (22). This polytope called C

								D	/ AB CD	is a 3-polytope which a
	projection in 	P,  , xy P,		is given in Fig. 15.
						  	5 2 2 6 t t 	  		P,x  6 t t  5     P,y	OB EF L L
					     5 6 t t	OB EF L L		  	5 2 2 6 t t 	  
	Fig. 15. Polytope of contact	D	C AB CD /
	4.3. Full condition

  ):

	r		 EL  3 BC 1 4 5 6 C tL t t t t 2 2 E C OB BC EF E L C L         	(26)
	r		 EL  3 BC 1 2 5 6 H tL t t t t 2 2 E OB H H BC EF L E L          


L Magnis

N Petit

GYROLESS BALL: ESTIMATION OF ANGULAR VELOCITY WITHOUT GYROSCOPE

In [START_REF] Magnis | Angular velocity nonlinear observer from single vector measurements[END_REF][START_REF] Magnis | Angular velocity nonlinear observer from vector measurements[END_REF], an algorithm is proposed to estimate the angular velocity of a rigid body using only vector measurements (i.e. measurements of unknown but constant directions). The cases of a single and two vectors measurements have been addressed using similarly structured state-observers.

In this note, we report experiments conducted on one such rigid body, here a wood ball, which is equipped with a smartphone. The direction sensors of the smartphone (magnetometer and accelerometer) are used in the algorithm and the angular rate estimate is compared against the readings from the smartphone gyro. Although neither the magnetic field nor the field sensed by the accelerometer are constant, the estimation algorithm works well.

The experiments were conducted on August 26th, 2015 and have been presented at the Congrès Français de Mécanique 2015. A video can be found on YouTube: https://youtu.be/L8fNeH347C8

Sensors

The smartphone under consideration is an Android Nexus 5 which embeds a triaxis gyros and accelerometers MPU-6500, and a tri-axis magnetometers AK8963. These sensors can be readily seen on the printed circuit board (PCB) of the smartphone, see Figure 1.

The "Sensor Kinetics Pro" Android app allows us to visualize the data produced by the sensors in real-time, to save them and to export them in .csv format. Although we do not have access to a data-sheet of the smartphone PCB, pictures of the components and their location on the PCB allows one to roughly identify the sensor axes orientations. As far as we could detect it, the three sensors are aligned and produce data in the body frame represented in Figure 1. The data are time-stamped. The sampling intervals are approximately 0.01 [s]. They vary with the processor usage.

Experimental results

We test the angular velocity observer on data produced by the aforementioned sensors. The reference vectors are the gravity field and the magnetic field. The smartphone is inserted inside a wood ball comprising a hole. Figure 2 pictures the smartphone partially inserted. The hole is then plugged with foam to ensure Date: 07/12/2015. stability of the smartphone. The rigid body under consideration is the system {ball + smartphone}. This device gives the user a lot of flexibility on the choice of rotation motion, as it can freely roll on the ground, according to its initial conditions.

In details, we implement observer (6) of [START_REF] Magnis | Angular velocity nonlinear observer from vector measurements[END_REF] with a the data produced by the accelerometer, b the data produced by the magnetometer. The tuning parameters are set to α = 1, k = 5. In the observer implementation, the torque governing the dynamics is arbitrarily set to τ = 0, which assumes that the rigid body is in free-rotation. The inertia matrix is (roughly) estimated to

J =   0.0878 0 0 0 0.0878 0 0 0 0.0878   [kg.m 2 ]
In this expression one has neglected the contribution of the smartphone to the matrix of inertia, this corresponds to the inertia of a homogeneous ball of radius 15

[cm] and mass density 0.69.

magnetos gyros + acceleros

Figure 1 .

 1 Figure 1. Left: printed circuit board of the Android Nexus 5. Identification of tri-axis sensors. Right: the body-frame is aligned with the principal axis of inertia of the smartphone.

Figure 2 .

 2 Figure 2. The smartphone is to be inserted in a wooden-ball via a hole.

The ball is then given some initial velocity and left rolling on the ground until it is stopped. Once the data are recovered from the smartphone, they are processed in pseudo real-time (i.e. from initial time to end-time) by the observer. The estimation results can be compared against the measurements of the angular velocity given by the rate gyros of the smartphone. In accordance with the theory, the gyros data are not used in the observer, they only serve as reference for final comparisons.

The estimation results are reported in Figure 3. It appears that the estimates are unbiased (though a bit noisy). Note: The experimental dataset can be requested to the authors by email.


Daniel Gautheret
email: daniel.gautheret@u-psud.fr

SCREENING GENOME SEQUENCES FOR KNOWN RNA GENES OR MOTIFS

Introduction

In the mid 1980s, as genomic sequences in the Genbank database reached their first million nucleotides, a few intrepid biologists turned their attention towards screening DNA sequences for the presence of known RNA genes (mostly tRNA at this time) and RNA motifs such as ribozymes or specific stem-loop structures. They soon realized this task was harder than finding homologues of protein-coding genes, since non-coding RNAs were defined not only by their primary sequence, but by a combination of sequence and secondary structure. Thirty years and hundreds of new RNA families later, this view still holds, albeit with wide variations between individual RNA molecules. For instance, some RNAs have sequence signatures that are stronger than those of most protein-coding genes. Fig. 1A shows an alignment of the Escherichia coli tRNAAla and its homologue in the Bacillus subtilis genome. Although the two species are separated by billions of years of evolution, their tRNAs have retained more than 90% identity. Such homologues may be efficiently identified using a sequence alignment program such as Blast [START_REF] Altschul | Basic local alignment search tool[END_REF]. One may argue, however, that tRNAs are very special molecules: they sustain selective pressure at many sequence positions due to their multiple and universal functions. But experience has shown that a variety of RNA genes, including fast evolving ones such as microRNAs, are also well characterized at the sequence level and can be efficiently detected using sequence comparison programs [START_REF] Menzel | The tedious task of finding homologous noncoding RNA genes[END_REF]. Yet, for many other RNAs, primary sequence information is not sufficient for accurate detection. Certain RNA genes such as the RNase P RNA present little if any conserved sequence, and most RNA structural motifs such as GNRA and E-loops, SECIS elements, U-turns or Rho-independent transcription terminators are composed of a few key nucleotides associated with a larger secondary structure scaffold. Fig. 1B shows consensus structures of the loop E [START_REF] Leontis | A common motif organizes the structure of multi-helix loops in 16 S and 23 S ribosomal RNAs[END_REF] and U-turn [START_REF] Gutell | Predicting U-turns in ribosomal RNA with comparative sequence analysis[END_REF] motifs, whose functions are conferred by a small set of conserved nucleotides embedded into a larger stem-loop structure. Such motifs cannot be identified using sequence comparison and require specialized programs that take into account both secondary structure and sequence constraints. Mining genomic sequences for RNA elements may thus entail using very different approaches depending on the type of element under scrutiny. Freyult et al. [START_REF] Freyhult | Exploring genomic dark matter: a critical assessment of the performance of homology search methods on noncoding RNA[END_REF] published a detailed comparison of the performance of various RNA detection programs using simulated sequence sets, and Menzel et al. [START_REF] Menzel | The tedious task of finding homologous noncoding RNA genes[END_REF] proposed a more qualitative evaluation of three programs based on "real-life" sequence sets. Here our purpose is not to evaluate programs but instead to help readers choose the most appropriate program for different situations and present practical execution examples for each of the programs. This article updates and expands our previous one in this series [START_REF] Legendre | Seeking RNA motifs in Genomic sequences[END_REF]. The proposed tasks are accessible to any junior bioinformatician or biologists with a basic knowledge of Unix commands. We provide links to official download sites for the various bioinformatics programs at the end of the Chapter. We assume that all the required programs are properly installed by a system administrator on a computer running the Linux or another Unix operating system.

Choosing the Right Search Program

RNA search programs can be roughly divided in two classes: descriptor-based and homology-based. Descriptor-based programs allow biologists to describe, using a defined syntax, the sequence and structure of an RNA motif and search this motif in genome databases. The most widely used program in this family is RNAMOTIF [START_REF] Macke | RNAMotif-an RNA secondary structure definition and search algorithm[END_REF]. Homologybased programs constitute a very diverse family. The most popular of them is Blast [START_REF] Altschul | Basic local alignment search tool[END_REF], which uses a single sequence as input and identifies segments of similarity in a sequence database. The HMMER program [START_REF] Bateman | The Pfam protein families database[END_REF] uses a sequence alignment as input and, consequently, is able to exploit variation and conservation in the gene or motif. The last category of homology search programs, exemplified by ERPIN [START_REF] Gautheret | Direct RNA definition and identification from multiple sequence alignments using secondary structure profiles[END_REF] and INFERNAL [START_REF] Nawrocki | Infernal 1.0: Inference of RNA alignments[END_REF], use multiple sequences and a description of secondary structure to perform the search. Knowledge of secondary structure allows these programs to identify parts of an RNA motif that are completely variable in sequence while retaining a constant structure. In Table 1, we present an overview of the qualities and domains of application of the different programs. Any serious RNA search should address the question of secondary structure. Knowledge of an RNA secondary structure opens the way to the most sophisticated and sensitive search programs. It may also provide important clues about the function of an RNA. If a conserved structure is present in the RNA under study, then it must be identified. There are well-known programs for predicting an RNA secondary structure from sequence. However, it is believed that on average about 30% of the base pairs predicted by these programs are incorrect. This false prediction rate is definitely too high if one wants to discover further RNAs based on this information. Much more reliable structure prediction can be achieved using a set of homologous RNAs and seeking a secondary structure that is compatible with all sequences in the set. A typical RNA search begins with a step of sequence collection to create this first training set. This can be done using a purely sequence-based program such as Blast. Once this first set is completed, a common RNA structure is derived using a program that simultaneously folds and aligns the sequences in the set. Until recently, such algorithms were considered as too computationally expansive to be applicable to real-life situations and biologists were reduced to align RNAs using a sequence-based multiple-aligner such as CLUSTALW [START_REF] Chenna | Multiple sequence alignment with the Clustal series of programs[END_REF], which often yield poorly aligned structures. However, recent progress allows reasonably accurate simultaneous alignment and folding. We chose the LOCARNA [START_REF] Smith | Freiburg RNA Tools: a web server integrating INTARNA, EXPARNA and LOCARNA[END_REF] program that presents the advantage of running on a web server. Once a structured alignment is available, it needs to be carefully inspected. Are we sure all RNAs in the set are homologues? Are they complete? Are there obvious misalignments? Some manual curation may be helpful at this stage. Then, one should identify the most conserved sequence and structure elements. This consideration dictates which program can be used to seek further homologues. If no obvious conserved secondary structure stands out, we recommend pursuing with further rounds of BLASTN and the HMMER program. If short sequence/structure motifs appear, then one may use RNAMOTIF to seek further instances of these motifs. If a conserved structure including several stems and loops is visible, then one might resort to more sophisticated programs that will take these structures into account. ERPIN is useful when a core of regular stems is present. INFERNAL is preferred when structures are more variable. The latter is the most general and powerful of all programs used in this Chapter, but it is also the slowest, although recent developments have brought its search speed closer to that of other programs. As further RNAs of the family are discovered using any of the above tools, they can be used to enrich and improve the structural alignment so that further search iterations can be performed. The steps of a typical RNA search can be summarized as in Fig. 2. All search programs accept sequence databases in the Multifasta format. Recent versions have practically no limitation in the size or number of sequences in the input file.

Assessing Search Specificity

Most sequence and motif search programs compute a score and an expectation value (Evalue) for each hit. An E-value tells us how many hits could be obtained by chance with the same score or higher in a random database of same size. An E-value above 1 or even 0.1 is likely to occur by chance. Any E-value below 10 -2 should in principle convince peer scientists that a biologically related motif has been found. However, keep in mind that a low E-value does not demonstrate biological significance, as it is subject to artifacts caused for instance by low complexity regions in genomes. Conversely, hits with a poor Evalue may turn out to be true homologues that have diverged a lot. Also note that Evalues depend on the database size. When a search program is run repeatedly on separate sequence files (e.g. with one file per chromosome), the E-value is valid only for each distinct file. RNAMOTIF is the only program in our protocol that does not automatically compute Evalues. In this case, or when E-values are not trusted for some reason, it is common to use an expected density of False Positives (e.g. FP per Mb) as a measure of specificity. To this aim, we may run the search against a randomized sequence database. Random sequences with a uniform nucleotide composition (25% each) are not advised, since background compositions (single-nucleotide and di-nucleotide frequencies) strongly affect the number of chance hits: for instance, an AU-rich motif is more likely to occur in an AU-rich sequence background. It is therefore important to reproduce in the negative control set the overall compositional biases of the target database. Several programs can do this. We will use shuffle, by S. Eddy (see Section Program Versions and Download Sites), that randomizes a sequence while preserving the single-and dinucleotide frequencies. The randomization procedure can be performed on the target database, so that each genome in the database is shuffled using its own compositional characteristics. When the search database is not decided beforehand, an alternative option is to shuffle the training set itself. This is a very conservative approach, as a randomized training set, owing to its limited sequence space, is more likely to give rise to false positives than any other sequence. If specificity has to be quantified, the negative set should, ideally, be large enough to find at least a few False Positives. By default, the shuffle program produces a sequence of same size as the input sequence. When using a short sequence such as the training set, the shuffling procedure must be iterated to produce a random sequence of the desired size. For instance:

shuffle -d -n 100 trset.fasta > trset.rnd

Here the sequence to be randomized is in file trset.fasta (multi-fasta format) and the output is redirected to trset.rnd. The -d option is used to preserve dinucleotide composition in the random sequence, and -n 100 specifies 100 shuffling iterations. The randomized sequences from each iteration round are concatenated and will thus be 100 times larger than the original sequence.

A Test Case: Looking for Homologues of a Bacterial sRNA

Our test challenge here is to discover homologues of a bacterial small RNA (sRNA) gene named RsaE, a 113-bp RNA first identified in Staphylococcus aureus [START_REF] Geissmann | A search for small noncoding RNAs in Staphylococcus aureus reveals a conserved sequence motif for regulation[END_REF][START_REF] Bohn | Experimental discovery of small RNAs in Staphylococcus aureus reveals a riboregulator of central metabolism[END_REF]. We will try to stay as close as possible to a "real life" situation, which means we assume only one instance of the RNA molecule is known first (here the S. aureus RsaE), and no structure information is available. This is the situation shown in the first box of Fig. 2. The initial, experimentally confirmed, sequence is entered in the Fasta format in a file named rsae.fa (see Section 7, "Supplemental data"), as follows:

>Rsae AUGAAAUUAAUCACAUAACAAACAUACCCCUUUGUUUGAAGUGAAAAAUUUCUCCCAUC CCCUUUGUUUAGCGUCGUGUAUUCAGACACGACGUUUUUUUAU
Our objective is to find all other instances of this sequence in bacterial genomes. Our search database will contain all complete genomes retrieved from the NCBI web site. An initial set of genomes from 1080 species was downloaded in 2010 and reduced to 732 species by retaining a single instance of each species. The sequences are saved in a single multi-Fasta file of about 2.5 Gb named 732bacs.fa (see Section 7). To simplify command lines, we assume all data files are present in the working directory.

BUILDING A FIRST TRAINING SET WITH BLASTN

Blast is the program of choice for a fast and easy identification of a first set of homologues. Biologists are familiar with the NCBI BLAST web interface. Here we use the command line version of the program. An advantage of using command line is that users can scan "private" sequences not yet included in the NCBI sequence databases. Furthermore, playing with parameters is faster with command lines. We assume a recent version of BLAST (such as 2.2.25) is installed on the reader's workstation.

A local BLAST run requires a properly formatted database. Formatting our database file 732bacs.fa is done using:

makeblastdb -in 732bacs.fa -dbtype 'nucl'

This command generates several index files required for BLAST runs. You should not worry about the names of these files, as you will always refer to the Fasta file name (here 732bacs.fa) for further operations. There are several programs in BLAST. Here we use BLASTN, which is used to scan nucleic acid databases for nucleic acid sequences. With our data files, a simple BLASTN command would look like:

blastn -db 732bacs.fa -query rsae-sau.fa

The command takes about 0.7 seconds to execute and produces a lot of output. It is pretty difficult to understand how many hits were exactly recorded. A simpler table output can be obtained using the -outfmt 7 parameter, such as in:

blastn -db 732bacs.fa -query rsae.fa -outfmt 7

Now we can count 12 hits from 12 different species. Are they all statistically significant though? The lowest scoring hit is from the Geobacillus strain WCH70. Its E-value of 3 x 10 - 11 is highly significant. We can control the E-value cutoff using the -evalue parameter, for instance by appending -evalue 0.01 to the BLAST command.

Can we find more distant homologues using BLASTN? This question is commonly addressed by reducing word size. In BLAST, both the database and the query sequences are broken down into k-length words (k-words). BLAST is able to quickly scan these word lists for exact or near-exact matches. Longer alignments with possible gaps and mismatches are generated only after one k-word match is found. This enables a very fast search at the risk of missing some hits that would contain no exact or near-exact k-word match. The default word size for BLASTN is 25. It can be changed using the word_size parameter. Shorter words will increase sensitivity, at the cost of longer runtimes. Let's reduce word size to 10:

blastn -db 732bacs.fa -query rsae.fa -outfmt 7 -evalue 0.01 -word_size 10

The command now takes 2.1 seconds to run, but finds 21 hits from different species, that is nine more than with the default word size. Note how significant these new alignments are. The lowest scoring hit is now found in Geobacillus kaustophilus, and has an E-value of 9 x 10 -10 , which is still very significant. Reducing word size to its minimum value (k=5) extends runtime to nearly 1 minute but does not reveal further solutions. This list of 21 BLAST hits is a good starting point for further analysis. We should now extract the corresponding sequences. The latest version of BLAST provides useful output controls. It is now possible to retrieve automatically the matched sequences and other items using parameter -outfmt 10 followed by the desired item names. For instance replacing -outfmt 7 in the above command with -outfmt "10 evalue sseqid sseq" will output the E-value, database sequence Id and matching sequence region in the database for each hit.

A visual inspection of these sequences reveals that several species share the exact same RsaE sequence. After discarding the redundant hits, our raw training set is composed of 14 distinct sequences. Before going any further, we should pay attention to the species distribution of these sequences. We have found RsaE in Staphylococcus, Macrococcus, Bacillus, Geobacillus, Anoxybacillus and Lysinibacillus. These genera all belong to the "Bacilli" class of the Firmicute phylum. This is an important observation that supports the relevance of the hits. Hits distributed in multiple distant phyla would be a warning sign of possible false positives.

ALIGNMENT AND STRUCTURE PREDICTION

As explained above, establishing the RNA secondary structure is an important step for further homology search. To this aim, we will use the LOCARNA program [START_REF] Smith | Freiburg RNA Tools: a web server integrating INTARNA, EXPARNA and LOCARNA[END_REF] through its Freiburg web server. LOCARNA belongs to a family of programs that simultaneously align and fold a set of RNA sequences. These programs are known to provide better RNA alignments than those relying only on sequence. The LOCARNA server is intuitive and works fine with default parameters. We can just paste the 14 Fasta-formatted sequences into the input window. To produce a Fasta file from the BLAST output, any text editor will do. As LOCARNA does not accept gap characters, we will use the editor to remove gaps too. The Fasta-formatted raw training set is available as Suppl. file trset.fa (see Section 7). The LOCARNA server displays the structural alignment after a few minutes. The output is shown in Fig. 3A with parentheses representing the secondary structure. This predicted structure is remarkably similar to the one established for S. aureus using chemical and enzymatic probing [START_REF] Geissmann | A search for small noncoding RNAs in Staphylococcus aureus reveals a conserved sequence motif for regulation[END_REF] (Fig. 3B). LOCARNA, like any other multiple alignment program, does not produce perfect alignments. For instance, the AAAC in the 5' region of G. kaustophilus RsaE would probably be aligned better with the AAAC segments found in all other sequences. To keep things simple, we will not try to improve this alignment manually, and we will proceed with the one in Fig. 3A. It appears from our alignment that RsaE homologues differ in their 3' regions. Either some RNAs are actually shorter, or this may be an effect of the BLAST program failing to extend the alignment to its actual ending. To settle this, it is necessary to retrieve the corresponding regions from the bacterial genomes and observe their sequence. This can be done using the coordinates from the BLAST output and any genome browser such as the Ensembl Bacterial Genome Browser (http://bacteria.ensembl.org/). We retrieved the missing 3' regions and realigned the sequences using LOCARNA. These 3' regions differed strongly from other RsaE sequences and aligned poorly both at the sequence and structure level (not shown). Therefore we can assume that these 3' regions are either missing or too different from the original RsaE sequence to be kept in our alignment. We also noted that the G. kaustophilus sequence lacked about 15 nt in its 5' region. We retrieved the corresponding segment from the genome sequence and used it for the LOCARNA alignment. The extended sequence aligned nicely with other RsaE sequences although parts of the long 5' stem appeared not to be conserved. The extra segment (italics) was incorporated into Fig. 3A.

The RsaE initial alignment presents a long stem on the 5' side. However, the sequence of this stem is highly conserved and there is little evidence of covariation that would indicate the base pairs are under evolutionary pressure. Therefore we cannot be absolutely sure about this secondary structure. The 3' part, which is missing in half the sequences is a typical Rho-independent transcription terminator, with a stem loop followed by a U-rich region. We may hypothesize that the RsaE sequences lacking this terminator use another termination mechanism.

SEARCHING WITH HMMER

A sequence alignment is a powerful means of conducting homology searches, even in the absence of secondary structure information. The main reason is that an alignment permits to distinguish variable regions that are not crucial for RNA function from conserved, crucial regions, and thus enables weighting matches accordingly. A good way to exploit an alignment is to use the HMMER package [START_REF] Bateman | The Pfam protein families database[END_REF]. HMMER converts the alignment into a profile Markov model that records both the statistics of base frequencies in each alignment column as well as that of consecutive bases in neighbouring columns. Once a profile HMM model is built from the alignments, matches to this model in sequence databases can be sought using the dynamic programming algorithm in HMMER. First, the training alignment should be formatted in the Stockholm format, as shown in suppl. file trset.stk. Here is the command to build the profile HMM:

hmmbuild trset.hmm trset.stk
This creates a file named trset.hmm, which can be used for scanning genome sequences. The search command is then entered as follows:

hmmsearch trset.hmm 732bacs.fa After a runtime of 42 seconds, we obtain 18 solutions. These include new Bacillus and Geobacillus species, as well as an entire new genus, Oceanobacillus, that was absent in the initial training set. This HMMER hit is shown in Fig. 4. The figure legend explains the different output items. The default mode of hmmsearch uses several heuristics to speed up the search, at a certain cost for search sensitivity. It is possible to switch off the heuristics using the --max option. This slows down the search by a factor of 10, but reveals no novel significant hit in our case. Values range between 0 and *, with 0 corresponding to a 0-5% probability, up to * for a 95-100% probability that this residue was produced from the hidden markov model created from the input alignment.

SEARCHING WITH RNAMOTIF

The RNAMOTIF program permits to describe an RNA motif, including sequence and secondary structure constraints, and to find all instances of this motif in a sequence database. Therefore, a major difference from the previous approaches is that we can now use secondary structure information. A typical RNAMOTIF descriptor is shown in Fig. 5A. Let's take a quick look at this language. The first section (parms), says that G:U pairs are allowed in addition to Watson-Crick pairs. The next section (descr) is the actual descriptor, listing all secondary structure elements. There are three types of elements: h5 stands for the 5' strand of a helix, h3 for the 3' strand of a helix, and ss for a single strand. When no other indication is provided, RNAMOTIF assumes pairing between h3 and h5 elements in a purely nested fashion, just as in a parenthesis notation. Non-nested helices (i.e. pseudoknots) require a special tag that modifies this default behavior. Associated to each h5, h3 or ss element is a list of optional parameters that specify minimum and maximum length, conserved sequences, conserved base pairs, etc. By default, RNAMOTIF reports all sequences matching the descriptor. Sometimes, however, additional filters are required that cannot be applied at the "descr" level. An optional score section is used for this purpose. In this section, tests can be performed using a C-like computer language, and matches are rejected if the tests fail. For an in-depth explanation of RNAMOTIF descriptors and scoring functions, refer to the detailed manual included in the program. Let's see how our knowledge of the RsaE RNA can be turned into an efficient descriptor. We will use only the most basic features of the RNAMOTIF language. To build our first descriptor (Fig. 5A), we make the typical assumption that a regulatory RNA has variable sequences in the secondary structure helices and conserved sequences in single-stranded regions, as the latter are supposed to undergo tertiary interactions. This is of course a simplistic view, but we would like to see how it performs. We set helix lengths as in the Fig. 3 consensus (size = 6, 4 and 6 nt, respectively from 5' to 3'). To account for shorter helices, we allow 1 mismatch per helix (mispair=1) and, since the default behavior of this parameter restricts mispairs to internal positions, we use ends='mm' to allow for mispairs at helix ends. We specify conserved sequences in several single-stranded segments. The most conserved single-stranded region is the apical loop where the conserved seven 3'proximal bases are specified by seq="ACCCCUU$". In the fifth "ss" region seq="^A" means that this strand should start with an A. We save this descriptor file as descriptor.txt and we are now ready to look for instances of this motif in our database. The basic RNAMOTIF command then looks like:

rnamotif -descr descriptor.txt 732bacs.fa

The run takes about 1 minute and produces a pretty large output. It is better to use a redirection to store hits in an output file. Also, due to variations in helix sizes and positions, RNAMOTIF tends to generate several overlapping solutions for each actual motif. The utility program rmprune is provided to filter out these overlapping solutions and retain only the highest scoring one at each site. This program is executed in combination with RNAMOTIF using the "pipe" sign (|) as follows:

rnamotif -descr descriptor.txt 732bacs.fa |rmprune > outfile

From now on, we will assume that every RNAMOTIF run is performed using |rmprune > outfile. We can now count 16 solutions produced by the RNAMOTIF run. Sample outputs are shown in Fig. 6. Most hits are from the genera Bacillus, Macrococcus, Staphylococcus, and Anoxybacillus. A few training set sequences such as the Geobacilli are missing, which indicates our descriptor is too strict in some respects. Indeed, we can see in the alignment that shorter helices #2 and #3 in G. kaustophilus do not match the descriptor constraints. Interestingly, only four sequences are from species outside our training set. Considering the large size of the database, this is a small number. These include one hit from Anaerococcus prevotii (Fig. 6A). Albeit distant in the phylogenetic tree, this species belongs to the Firmicutes, as do the other known RsaE sequences, and is therefore an interesting candidate to consider. However, we observed that the sequence just following the motif is not as CU-rich as in other RsaE sequences. The other new hits are from phyla that are very distant from Firmicutes and therefore unlikely to contain an RsaE homologue. We can reasonably assume they are false positives.

To assess the significance of new RNAMOTIF hits, we advise using a shuffled database as a control. Here we will randomize our bacterial genome database while preserving the dinucleotide frequencies, using the shuffle command as this:

shuffle -d 732bacs.fa > 732bacs.rnd

Running the same RNAMOTIF descriptor on the shuffled database 732bacs.rnd produces 6 hits. This is a clear warning that the few hits obtained with real data can be false positives.

Examination of the RsaE alignment suggests that sequence alone may describe RsaE better than secondary structure. The stems are indeed highly conserved in sequence and show no sequence covariation to support their base-pairing potential. Let's see how efficient our RsaE search could be by using solely primary sequence constraints. All RsaE members harbor a constant AAAC motif followed by a highly conserved CU-rich sequence that is repeated in tandem. The descriptor in Fig. 5B describes one copy of the tandem motif in a single line. The syntax for sequence constraints in RNAMOTIF is that of Unix regular expressions. The "^" sign indicates the beginning of the element and the string .\{2,10\} describes a sequence of 2 to 10 unspecified bases, described by the dot character. This simple motif performs surprisingly well. After 75 seconds, RNAMOTIF identifies 27 hits, covering all the species in the training set plus several other species from the same genera and one hit in Pectobacterium wasabiae, an enterobacterium. Enterobacteria are very distant from Firmicutes, which makes this hit an unlikely RsaE candidate.

Running the descriptor in Fig. 5B against the shuffled genome database produces a single hit, which shows that our simple sequence motif is more specific than the complex secondary structure motif in Fig. 5A. Can we now relax this descriptor and find additional candidates? RNAMOTIF allows for mismatches in sequence motifs, using the mismatch=n parameter (n being the number of allowed mismatches). However, the mismatch parameter is compatible only with regular expressions of constant length. Therefore we cannot allow for mismatches in a variable region of size 2 to 10 as above. We thus have to rewrite our descriptor using distinct ss elements, as shown in Fig. 5C. The descriptor in Fig. 5C that tolerates one mismatch in the CU-rich region takes longer to run (2 minutes) than the exact descriptor and produces 71 solutions. In addition to new species from the same genera as above, these include a number of apparent false positives from Cyanobacteria or Enterobacteria, plus new Firmicutes candidates of the genera Streptococcus and Clostridium (Fig. 6B). These results should be considered with caution though, as the false positive ratio has increased significantly. We now find 35 hits in the randomized database (not shown).

A possible control to evaluate the new candidates is to analyze the sequence downstream of the motif. We can see in the training set (Fig. 3) that the downstream sequence is a near exact repeat of the central conserved region. Flanking sequences around a given motif are easy to extract using RNAMOTIF by adding new "ss" elements of desired length with no sequence constraint. We used this to obtain the 30 nt of downstream sequences. This sequence is not UC-rich in Streptococci, but it is in one of the Clostridia, C. phytofermentans (Fig. 6B). Therefore this candidate may be considered further. We can see that RNAMOTIF is a valuable tool to test ideas about what determines an RNA motif. We have developed highly specific descriptors in just minutes, using either primary or secondary structure information. One important limitation of RNAMOTIF though is its lack of an objective scoring function. Although the program accepts user-defined scoring functions, their quality will depend on the user's appreciation of the RNA motif. The subsequent protocols address this issue.

SEARCHING WITH ERPIN

ERPIN circumvents the descriptor design problem by extracting information directly from a structure-annotated alignment and building log-odd score profiles for each stem and single-stranded segment in the alignment. First, we need to convert the alignment into the ERPIN format. Fig. 7A shows how the first lines of our RsaE alignment should look like. The first three lines provide secondary structure information. In this example, the secondary structure is the same as the parenthesized structure in Fig. 3, but each structural element (single-stranded or helical) is assigned a different number that is read vertically. When distant segments have the same number (e.g. 02 or 04 in this example), the corresponding element is understood to be helical, otherwise it is single-stranded. The secondary structure lines can be created manually using a text editor, or automatically from a parenthesized secondary structure file using the parent2epn.pl script provided in the ERPIN distribution. The lines following the structure information contain the aligned sequences in the classical FASTA format.

The initial parenthesized secondary structure had a single long single-stranded segment composed of elements 17 and 18. Element 18 (the terminator) is not present in all sequences while element 17 is present everywhere and is highly conserved. To enable searching for the conserved part while ignoring the other, it is preferable to split this single-stranded segment into two different elements. We thus created element 18 using a text editor. Our ERPIN formatted alignment was saved as trset.epn (suppl. data).

The tview command is useful to visualize constraints in the sequence alignment:

tview trset.epn
The result is shown in Fig. 7B. ERPIN users must decide which part of the alignment is used for the search. A direct search for a multiple-element region could be prohibitive in terms of CPU and memory usage, especially when elements vary a lot in size. Here, the alignment suggests that the region from the 5' side of helix 02 to single strand 17 is the most conserved (asterisks on top). We will thus launch the search using this region (we assume all ERPIN commands below are complemented by -silent > outfile in order to send output to file outfile, devoid of trace text).

erpin trset.epn 732bacs.fa -2,17 -nomask

The minus sign before 2 is used to specify the 5' part of stem 02. The -nomask argument indicates that all elements in the region should be used for the search. We will see later how we can mask/unmask elements. The ERPIN command takes as much as 4 days to be executed, producing 21 solutions. A sample is shown in Fig. 8A. The solutions cover all species from the training set, plus a few closely related species. The search is highly specific, with all E-values around 10 -20 , however it takes a long time and does not find distant homologues. Using the whole structure for search is thus not very practical here (4 days run!), and it can become worse for RNA motifs containing more elements, more gaps or fewer conserved regions than the RsaE RNA. More gaps in single-stranded segments (helices cannot contain gaps) and more variation in the alignment means longer search times and reduced search specificity. One way to reduce runtime is to restrict the search to a smaller conserved region. For instance we can use the region from single strand 9 to single strand 17:

erpin trset.epn 732bacs.fa 9,17 -nomask

The search time is now down to a more reasonable 10 hours and ERPIN finds the same 21 hits as before. However, because of the different search region, the hit sequences and Evalues are different (Fig. 8B). It is relatively easy to further decrease CPU time by implementing a search strategy. The main idea is to perform a stepwise search using masks. In the first step, one should mask most of the selected region and retain only a few key elements. ERPIN will disregard all masked elements and restrict search to the unmasked part. The best parts to unmask during the first step are those conserved elements that occur rarely in the database and are in close proximity to each other so that they are quickly identified. Then, at each successive step, one unmasks additional elements for ERPIN to look for. ERPIN will consider these latter elements only after elements have been identified in the first round. Any number of search steps can be specified on the same command line. This prints out only those hits with an E-value below 1. There are 21 such hits, and the list is the same as before, hence no new significant hit was found. However, one of the two high E-value hits (E = 2.82) is from Oceanobacillus iheyensis, a member of Bacilli that was not in the training set and was also found by HMMER (Fig. 8C). This confirms the strength of this candidate. The output covariance model is written in file trset.cm. In principle, a database scan is possible using this file. However, in order to benefit from E-value calculation and performance-enhancing filters, the model must be processed using the cmcalibrate command, as follows:

SEARCHING WITH INFERNAL

cmcalibrate trset.cm
This command is more CPU-intensive than the cmbuild command. It requires about 40 minutes in our example, but it can take much longer with complex multi-helix motifs. When it terminates, the cm file is changed to incorporate calibration information. We can now run the search over our 732 genome database using the cmsearch command, as follows:

cmsearch trset.cm 732bacs.fa

As cmsearch is relatively slow, we highly recommend using output redirection and the nohup command that enables a program to keep running even if the terminal window is closed:

nohup cmsearch trset.cm 732bacs.fa > outfile &
Here the & sign is used to execute the task in the background. This run takes about 10 hours and finds 26 solutions. In addition to the Firmicute species found by other programs, four identical hits are found in different Listeria species (Fig. 9). Listeria are Firmicutes and E-values around 0.01 indicate that these are serious candidates. Some Listeria hits were found with HMMER but with E-values above 1.

Conclusion

I have presented five different protocols for RNA motif search, based on BLAST, HMMER, RNAMOTIF, ERPIN and INFERNAL. In absolute terms, the usefulness of an RNA search protocol lies primarily in its precision, or balance between specificity and sensitivity. However practicality is of course important too. In the above examples, the INFERNALbased protocol was probably the most effective as it found significant candidates that none of the other approaches identified, while enabling reasonable runtimes and command line complexity. However, the diversity of RNA motif search situations is such that alternative protocols should always be considered. Motifs with little or no conservation information are best tackled using RNAMOTIF, motifs based on a pseudoknot are best dealt with ERPIN or RNAMOTIF, and motifs with highly variable or uncertain secondary structures should be addressed preferentially using HMMER or BLAST. The biologist's eye thus remains instrumental in the process of RNA motif search.

Supplemental data

An archive containing all supplemental files is available at http://rna.

Fig. 1 .

 1 Fig. 1. RNA motifs vary widely in their evolutionary constraints. A. A strong primary sequence constraint is visible is this Blast alignment of tRNAAla from E. coli (top) and B. subtilis (bottom). B. Predominant secondary structure constraints are visible in the consensus loop E and U-turn motifs. Thin vertical lines indicate Watson-Crick base pairs, open circles indicate non-canonical base pairs.

Fig. 2 .

 2 Fig. 2. Outline of a typical RNA motif search protocol. The protocol starts with a BLASTN sequence collection followed by structural alignment (align & fold) and the use of a specialized search software.

Fig. 3 .

 3 Fig. 3. A. An aligned and folded training set of RsaE sequences (Stockholm format) obtained from the initial BLAST screen and aligned using LOCARNA. Parentheses in the bottom line indicate basepaired positions, highlighted by shaded boxes. Bases in italics in "Geobacillus_kaus" were absent in the initial BLAST hit and manually inserted before alignment. B. The S. aureus RsaE secondary structure derived from enzymatic and chemical probing [13].

Fig. 4 .

 4 Fig. 4. Part of an HMMER output. The local hit is shown with 5 lines. The "CS" line is just a reminder of the secondary structure in the .stk file, which is not used by HMMER, the next three lines show the consensus profile, matching residues and target sequence, and the last line shows the posterior probabilities of each aligned residue which is an indication of alignment reliability.Values range between 0 and *, with 0 corresponding to a 0-5% probability, up to * for a 95-100% probability that this residue was produced from the hidden markov model created from the input alignment.

Fig. 5 .

 5 Fig. 5. RNAMOTIF descriptors for RsaE. A. This descriptor uses three stems with variable sequences and six single strands containing conserved sequences. B. This descriptor uses a single 60-nt single strand with a long conserved 3' motif. C. This descriptor splits the previous single strand in three parts. In the last part, a mismatch is allowed in the consensus motif.

Fig. 6 .

 6 Fig. 6. Extracts of RNAMOTIF outputs. Hits are formatted as follows. Line 1: title line of Fasta file entry. Line 2: Sequence name, score (here our descriptors do not include a scoring section hence all scores are 0), strand ("0" indicates plus strand and "1" minus strand), hit position, hit length in nt. Line 3: hit sequence, with stems and single-stranded elements separated by spaces. A. A hit from the descriptor in Fig. 5A. B. Two hits from the descriptor in Fig. 5C.

Fig. 7 .

 7 Fig. 7. ERPIN data formats. A. Typical ERPIN input file (.epn file) showing secondary structure (first three lines) and aligned sequences (next lines). The format is explained in text. B. A .epn file displayed using the tview command. Each secondary structure element, S (single-stranded) or H (helical), is highlighted with element number (read vertically) and sequence conservation (asterisk).

Fig. 8 .

 8 Fig. 8. Extracts of ERPIN outputs. Hits are formatted as follows. Line 1: title line of Fasta file entry. Line 2: Forward (FW) or Reverse (RC) strand, a fixed tag, hit positions, PWM score and E-value. Line 3: hit sequence. Structural elements are separated by dots. Consecutive dots indicate that corresponding elements were not included in the ERPIN command. Gaps inside elements are used to align the hit sequence onto the training alignment. A, B, C: results of different ERPIN runs (see text).

Fig. 9 .

 9 Fig. 9. Extract of an INFERNAL output. Lines 1-4 provide general information about the hit shown below: coordinates of hit in the alignment (query), coordinates of hit in the database sequence (target), score, E-value, P-value and GC-content of hit. The alignment shows the hit sequence (bottom) and the consensus sequence from the training alignment (top). The middle line shows conserved bases. Special characters on top of the alignment show the secondary structure in a bracket notation.

Table 1 .

 1 Domains of application of different RNA search programs

	Program	Good for		Not good for	Pros/cons
	Descriptor-based	Short	motifs	Complete RNA genes	Quick / Requires
	(RNAMOTIF)	(terminators,	hairpin		learning descriptor
		loops, internal loops etc.)		language
	BLAST	Long or short RNA	Short motifs, distant	Quick and easy /
		genes, domains of genes	homologues	boundaries	not
					reliable
	HMMER	Long or short RNA	Short motifs, distant	Quick / a bit
		genes, domains of genes	homologues	more complex than
					Blast
	ERPIN	RNA motifs or genes <	Long RNA genes,	Quick	/
		200 nt, domains of RNA	highly variable structure,	somewhat tricky to
		genes, RNA with important	irregular stems	use
		pseudoknot
	INFERNAL	RNA motifs or genes <	Long RNA genes,	Most sensitive
		200 nt, distantly related	highly variable structure,	method / Slow
		RNAs		pseudoknots
	3. Overview of the RNA Search Procedure

  In this command, we first unmask element 17, then unmask the whole region. The search time is reduced to 26 minutes and this does not change the number of solutions. To increase search sensitivity and find more distant homologues, we now need to play with scores and pseudocounts. Default cutoff scores are set for each element in order to capture the lowest scoring member of the training set. Pseudocounts are used to allow for additional variation in regions that have as yet displayed very little variation, such as strand 17 in Fig.7. Playing with scores and pseudocounts can become relatively complex and is beyond the scope of this Chapter. Interested readers should refer to the Erpin documentation. However, the Perl script erpincommand.pl provided with the ERPIN distribution can suggest a reasonable ERPIN command using cutoff and pseudocounts to achieve a balanced specificity / sensitivity ratio. The script is run as follows:Running this command (replacing DATABASE by our actual database file name) produces 23 hits after a 176 minutes run. Some hits have poor E-values (1 or more). To filter hits, we can use the script readerpin.pl, as follows, assuming ERPIN output was redirected to

	erpincommand.pl trset.epn
	and it suggests the following command:
	erpin trset.epn DATABASE -09,18 -add 09 17 -pcw 51 -cutoff 41.96
	outfile:
	readerpin.pl < outfile -e 1
	For our RsaE
	motif, let's try the following 2-step strategy:
	erpin trset.epn 732bacs.fa 9,17 -umask 17 -nomask

  INFERNAL uses covariance models to describe a structural RNA sequence alignment and subsequently locates instances of this model in genome database. A covariance model is a probabilistic description of the sequence and structure constraints in the training set. It is more general than the model used by ERPIN with the caveat that it cannot take pseudoknot constraints into account. Early versions of INFERNAL required amounts of computing power that kept the program out of reach of most biology labs. However, version 1.0 implemented heuristic filters that increased search speed considerably and rendered INFERNAL searches accessible to anyone with basic command line skills. First the following command should be executed in order to produce the covariance model from the Stockholm structural alignment trset.stk we created above in the HMMER section.

	cmbuild trset.cm trset.stk

  igmors.upsud.fr/suppl_data/ rsae.fasta: initial RsaE sequence from S. aureus. trset.fa: training set of RsaE sequences in Fasta format trset.stk: training set of RsaE sequences and conserved secondary structure in Stockholm format trset.epn: trainings set of RsaE sequences and conserved secondary structure in ERPIN format -732bacs.fa: Fasta-formatted database of 732 bacterial genomes (2.5 Gb uncompressed). 8. Program Versions and Download Sites BLAST 2.2.25 from ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ LOCARNA web server: http://rna.informatik.uni-freiburg.de:8080/LocARNA.jsp HMMER 3.0: http://hmmer.org/ ERPIN 5.4 : http://rna.igmors.u-psud.fr/erpin/ RNAMOTIF: v3.0.7 http://casegroup.rutgers.edu/ INFERNAL : v 1.0.2. http://infernal.janelia.org/ SHUFFLE. Part of SQUID 1.9 library: ftp://selab.janelia.org/pub/software/squid/ 9. Acknowledgements I thank Magali Naville and Antonin Marchais for their help in running the different software and Alain Denise for his assistance with Infernal runs and useful comments on the manuscript.