Faiez Zalila

Xavier Crégut

Marc Pantel

Leveraging formal verification tools for DSML users: a process modeling case study

Keywords: Domain specific modeling languages, Model formal verification, Behavioral properties, Translational semantics, Verification feedback

In the last decade, Model Driven Engineering (MDE) has been used to improve the development of safety critical systems by providing early Validation and Verification (V&V) tools for Domain Specific Modeling Languages (DSML). Verification of behavioral models is mainly addressed by translating domain specific models to formal verification dedicated languages in order to use the sophisticated associated tools such as model-checkers. This approach has been successfully applied in many different contexts, but it has a major drawback: the user has to interact with the formal tools. In this paper, we present an illustrated approach that allows the designer to formally express the expected behavioral properties using a user oriented language -a temporal extension of OCL -, that is automatically translated into the formal language; and then to get feedback from the assessment of these properties using its domain language without having to deal with the formal verification language nor with the underlying translational semantics. This work is based on the metamodeling pattern for executable DSML that extends the DSML metamodel to integrate concerns related to execution and behavior.

Introduction

TOPCASED 1 is a project2 started in 2005 in the French "Aerospace Valley" cluster that gathers academic and industrial partners [START_REF] Farail | The TOPCASED project: a toolkit in open source for critical aeronautic systems design[END_REF]. TOPCASED is dedicated to the development of an open source Computer Assisted Software Engineering (CASE) tool for the development of safety critical aeronautics, automotive and space embedded systems. Such developments will range from system and architecture specifications to software and hardware implementation through equipment definition.

TOPCASED provides modeling languages, both domain specific (SAM, EAST-ADL, SAE AADL, SDL 3 and XSPEM 4 ) and general purpose (SYSML, UML, etc.) and associated tools like graphical and textual editors, documentation generators, validation through model animation, verification through model checking, version management, traceability, etc. TOPCASED relies on MDE generative technologies to define the languages and build all these tools for all these languages. It is thus an MDE platform both for building system models and for building the platform itself. MDE technologies used in TOPCASED for defining and tooling languages are centered around Ecore 5 and configuration models taken as inputs by generative or interpretative tools.

Because the TOPCASED toolkit addresses safety critical systems, Validation and Verification activities are of primary importance and should be performed as early as possible at design time on the various models, both to reduce the development costs and to provide higher quality systems.

Validation is performed through model animation [START_REF] Crégut | Generative technologies for model animation in the TopCased platform[END_REF]: the designer builds a model using a graphical editor and can execute it according to scenarios. The runtime data produced by these executions is displayed as decorations of the graphical representation of the model or thanks to a dedicated view. Model animation is thus very similar to source level debugging for software. Scenario driven model execution runs through a single path in the set of all possible executions for the model. The use of several scenarios provides a coverage of the various possible executions but this validation is usually not exhaustive.

On the contrary, verification aims to check whether a property holds for all possible executions of the model. Model-checkers are dedicated tools for that purpose. These tools usually rely on two formal verification languages: one to model the behavior of the system and one to express the properties to check. For example, the TINA toolbox [START_REF] Berthomieu | The tool TINA -Construction of Abstract State Spaces for Petri Nets and Time Petri Nets[END_REF], available in TOPCASED, relies on Time Petri nets (TPN) for the behavior and State-Event Linear Temporal Logic (SE-LTL) for the properties. Thus, the use of such model checking tools requires to translate the system business domain model into an equivalent behavior model in the considered formal verification language and to express the system requirements as properties. Furthermore, results are obtained on the formal side as execution traces and have to be translated back into the system domain. This is a well-known technique called translational semantics. Nevertheless, even if the translations are automated, they are often defined in an ad hoc way, specific to the considered business domain. Furthermore, system requirements are most of the time directly written as formal properties, in the verification tool domain and not in the system domain. Thus, the designer must have a good understanding of: a) the various domain languages; b) the behavior and property languages from the various tools; and c) on the translation scheme used to go from one to the other in both directions. Verification is thus a difficult activity requiring many abilities that are generally not available to the casual business domain designer. Our purpose is to provide methods and tools in order to ease the integration of model checking in MDE toolchain. This integration will provide seamless verification facilities to the business domain designer without requiring him to deal with target verification language and associated model-checkers. We will describe a partly automated MDE driven tool chain for expressing the system requirements in the business domain language, translating the requirements to model checking tools property languages, and translating the failure execution traces back to the designer's world. Our contribution consists in reifying elements involved in the semantics of a DSML in Fig. 1: General approach of a translational semantics with feedbacks order to ease, and partially automate, the different translations that are summarized in Figure 1. More precisely, it includes:

1. the use of MDE technologies on both sides, business and formal verification domains based on a metamodel architecture that combines concerns related to model execution, including runtime information and the stimuli that make the model evolve. This metamodel guides the definition of the translational semantics and simplifies the production of business domain feedback to the end-user. 2. a user dedicated language for the expression of business properties.

automatic translation of business properties into formal verification domain prop-

erties based on the translational semantics. 4. automatic translation of the verification results obtained in the verification domain to the business domain. When a property does not hold, the obtained counter example is presented to the user either as a business domain scenario or a snapshot of the model completed with runtime information.

The approach is illustrated on a case study which concerns modeling of process using a process description language derived from the SPEM OMG Standard [START_REF]Software & Systems Process Engineering Metamodel (SPEM)[END_REF]. The paper is organized as follows. Section 2 presents the case study from the enduser viewpoint. It defines some constraints to assess and the expected feedback. Section 3 describes the formalism used for modeling processes, the language of expression of temporal constraints and extensions made on the DSML to be able to capture verification results. Section 4 presents the formal language and tools. Section 5 describes all required transformations for process verification and verification feedbacks. Section 6 considers related works and the last section concludes.

End-User Concerns

This section presents the business domain -process modelling -considered in the case study and the concerns of end-users. We first present the kind of process models the end-user wants to build. Then we explain the kind of properties he wants to check on his models. Finally we describe the feedback the end-user expects from verification tools in order to get insight on the errors the models may contain. -------------

Business Models

- count = 3 1 startToStart
Fig. 2: A business development process Figure 2 shows an example of a process model. It corresponds to a simplified development process composed of four activities, each represented in an ellipse: Programming, Designing, Test case writing and Documenting. Arrows between activities indicate dependencies: the target activity depends on the source activity. The label specifies the kind of dependency. The word before the "To" is the state that should have been reached by the source activity in order to perform the action on the target activity, action which appears after the "To". For example, the "fin-ishToStart" dependency between Designing and Programming means that Programming can only be started when Designing has been finished. Documenting and Test-CaseWriting can start once Designing is started (startToStart) but Documenting cannot finish if Designing is not finished (finishToFinish). The dependencies put between Programming and TestCaseWriting enforces a test driven development: programming can only starts when test cases are already started and, obviously, test case writing can only be finished when programming is finished in order to take into account test coverage.

Rounded rectangles represent the number of available resources (2 Designers, 3 Developers and 3 Computers). Dashed arrows indicate how many resources an activity requires. Programming needs two developers and two computers. Resources are allocated when an activity starts and freed when it finishes.

These processes are deliberately simplified to avoid overloading this presentation but time constraints or hierarchical decomposition of activities could be added.

User verifications

To validate or to verify a model, the user may check that properties derived from the system requirements hold on that model. We focus on behavioral properties, properties that concern the evolution of the model over time. Static properties are also important for the end user but they can easily be included in the editing tool using for example an OCL checker.

The user may be interested in general properties not specific to a given process model. For example, he may want to check whether a process model may finish or not (P 1 ). A process finishes if all its activities finish while respecting constraints imposed by dependencies and resource allocation. If these properties hold, the user may want to get a terminating scenario and use it to pilot the process execution.

The user may also want to verify properties that are specific to a particular process model. As an example, he might want to know if in all cases Documenting is finished before Designing is finished (P 2 ).

Verification feedback

Once the end user has defined his model and expressed his requirements through properties, he wants to have feedback on the assessment of those properties. If a property evaluates to true, then the requirement is fulfilled. But if a property evaluates to false, the user expects to have feedback in order to understand why the property does not hold. For example, a counter example may be exhibited. Obviously, this counter example should be expressed at the business domain level.

For instance, using the example shown in Figure 2, property P 1 does not hold and there is indeed a deadlock during process execution. The user can be provided with a counter example that explains the deadlock as a scenario like the one of Figure 3 which lists the actions (start or finish) applied on activities. The deadlock is due to the fact that Programming cannot be started because a Computer is missing. If a computer was added, then the P 1 requirement would hold. The property P 2 does not hold. Indeed, it is possible to finish Designing before Documenting is finished. A possible scenario is shown on Figure 4 (counter-example). The user may want to play those scenarios using a model animator like the one developed in the TOPCASED project [START_REF] Combemale | Introducing Simulation and Model Animation in the MDE TOPCASED Toolkit[END_REF].

Business Metamodeling

Metamodels generally focus on business domain concerns and do not take into account other elements required to execute a model. As model execution may be of interest for most of the modeling languages, especially in the context of safety critical systems, we have defined a general solution to describe all the data required to define an execution semantics for any executable modeling language [START_REF] Combemale | Simulation et vérification de modèle par métamodélisation exécutable[END_REF]. It is a kind of metamodeling pattern that may be used from design time to run time. This pattern has been applied in the TOPCASED project to build animators [START_REF] Crégut | Generative technologies for model animation in the TopCased platform[END_REF] that allow validation of SysML/UML State Machine and Activity diagrams or SAM (an automate-based language used by Airbus) models. It is also helpful to ease the definition of forward and backward transformations toward verification languages in order to get back failure scenarios from model checkers.

The pattern advocates to structure an executable DSML metamodel in such a way that the different concerns are stressed: the business domain, the queries a user may ask on a model to assess it satisfies its requirements (i.e. the model business properties), the stimuli that make the model evolve. The corresponding XSPEM metamodel is shown on Figure 5 and detailed in the next paragraphs.

XSPEM Domain Definition Metamodel (DDMM)

A metamodel defines the concepts (metaclasses) of the business domain addressed by the DSML and the relationships between them (references). In the executable meta-model pattern defined in TOPCASED, this reference metamodel is called the Domain Definition MetaModel, DDMM . The DDMM of XSPEM is shown on Figure 5 (package named DDMM at the bottom). It defines the concepts of process (Process) composed of a set of (1) workdefinitions (WorkDefinition) that model the activities (described in section 2.1) performed during the process, (2) worksequences (WorkSequence) that define dependency relationships between workdefinitions and (3) resources (Resource) allocated to activities (Parameter).

Obviously, this metamodel could be extended with well-formedness rules for example using OCL to express constraints not captured by the metamodel definition (names of workdefinitions have to be unique, worksequences should not be reflexive, resources counts should be positive, etc.). This aspect related to the static semantics of the DSML is not in the scope of this paper.

XSPEM Query Definition Metamodel (QDMM) and Formal Expression of Requirements

End users' behavioral properties usually rely on information that is not directly available in the DDMM because they only exist when the model is executed (runtime information). For example, the previous properties rely on the state of a workdefinition: started or finished. As usual, this information has to be reified. Thus, we have chosen to extend the DDMM with a new metamodel which describes the queries the user can conduct on his models. We call it QDMM (Query Definition MetaModel). For XSPEM, queries such as isStarted and isFinished can be applied on a WorkDefinition (see top right of Figure 5). A query isFinished is defined on Process in order to model the end-user business requirement, (P 1 ), defined in section 2.2. To be checked, requirements of section 2.2 have to be formally expressed. OCL is not well suited for that purpose because it only allows the specification of structural properties and some Floyd-Hoare behavioral properties for methods. Nevertheless, it is now a widely known language and a few temporal extensions of OCL have been proposed in order to specify event-based behavioral properties whereas OCL only targets function-based properties. We have chosen to rely on Temporal OCL and especially on the proposal from [START_REF] Ziemann | An Extension of OCL with Temporal Logic[END_REF] as the syntax of this extension is quite natural for OCL users. It introduces usual future-oriented temporal operators such as always, sometimes, next, existsNext as well as their past-oriented duals.

Here-after are the expression of the P 1 and P 2 requirements identified in section 2.2. They rely on the DDMM but also on the queries defined in the QDMM.

XSPEM Event Definition and Trace Management Metamodels (EDMM & TM3)

The DDMM and QDMM allow to express the requirements but we also need to show to the user the results obtained in the formal verification domain. A Snapshot can be expressed using QDMM because it represents all runtime information of interest to the user. To express the scenario corresponding to a counter example (like the one shown in Figure 4) we define two other metamodels that also extend the DDMM. The first one is the Trace Management Metamodel (TM3). It allows definition of a scenario as a sequence of runtime events -a stimulus that makes the model evolve. The TM3 is independent of any DSML. On the contrary, the Event Definition Metamodel (EDMM) is specific to a DSML and defines its runtime events. For instance, runtime events for XSPEM include "start a workdefinition" and "finish a workdefinition".

Formal Level Metamodeling

In order to represent the semantic data and ease the exchange of verification results with business domain models, the metamodeling pattern is also applied on the formal language, TPN in the case of the TINA toolbox [START_REF] Berthomieu | The tool TINA -Construction of Abstract State Spaces for Petri Nets and Time Petri Nets[END_REF]. Like XSPEM, the TPN metamodel is composed of several parts (figure 6). The DDMM describes a Petri net (PetriNet) composed of nodes (Node) that denote places (Place) or transitions (Transition). Nodes are linked together by arcs (Arc). Arcs can be normal ones or read-arcs (ArcKind). The attribute initialtokenCount specifies the number of tokens consumed in the source node or produced in the target one (in case of a read-arc, it is only used to check whether the source place contains at least the specified number of tokens). Finally, a time interval can be expressed on transitions.

The QDMM defines only one query corresponding to the number of token stored in a place (tokenCount). We can define other queries like for example fireableTransition corresponding to the set of fireable transitions in a petri net. The SDMM (State Definition Metamodel) is an implementation of the QDMM, in this case a trivial implementation that defines an attribute for each query.

Finally, the EDMM defines only one event FireTransitionEvent and, obviously, the TM3 is the same as the one presented for XSPEM, as it is DSML-independent. The last steps concern the transformation level of Figure 1. It consists first in defining the translational semantics, that is translating an XSPEM model into a Petri net one. Then, TOCL properties have to be translated into LTL formulae so that they can be checked by the TINA toolbox. Finally, results obtained on the formal verification domain have to be translated back into the business domain. One can notice that business metamodel (or model) TOCL invariants may be expressed to assert that they are preserved by the translational semantics.

All these transformations are defined at the metamodel level. They have thus only to be defined once by DSML and formal domain experts and the end-user can use them for any of the XSPEM models. Several translational semantics may be defined for XSPEM according to the level of details in the execution that we want to model and the kind of properties we want to assess. Thus, we advocate in [START_REF] Combemale | A property-driven approach to formal verification of process models[END_REF] that defining the translational semantics should be property-driven to favor the definition of a minimal semantics, that will allow to answer to the questions the user may ask about his models.

As the QDMM metamodel records all the aspects of interest to the end user, the expert has to verify that all the queries may be expressed using the formal language translation. Thus, the QDMM can be used as a guide to write the translation. For example, in the formal language, one should be able to determine if a workdefinition is started or finished as these queries are part of the QDMM. A WorkDefinition is thus translated into four places characterizing its state (notStarted, started, running and finished) linked by two transitions. These transitions model the actions that we want to observe on a workdefinition: one can start a workdefinition and then finish it. A workdefinition is considered started if it is either running or finished. This is recorded by the place named started.

AWorkSequence becomes a read-arc6 from one place of the source workdefinition (either started or finished) to a transition of the target workdefinition (either start or finish) according to the kind of WorkSequence (linkKind attribute). A resource becomes a place whose initial marking (initialtokenCount) corresponds to its count. Each Parameter element is translated into two arcs, the first one to take resources when the concerned workdefinition starts and the second one to release them when the workdefinition finishes.

Figure 7 contains the Petri net model resulting from the application of the translational semantics on a part of the XSPEM model from Figure 2 (Designing and Documenting workdefinitions and worksequences between them as well as the Developer resource).

The ATL transformation language [START_REF] Jouault | Transforming Models with ATL[END_REF] has been used to implement this translational semantics. First, an ATL module describes the transformation from an XSPEM model to a Petri net model [START_REF] Bendraou | Definition of an executable spem 2.0[END_REF] (not shown here 7 ). Then, an ATL query 8

Translating TOCL to SE-LTL

Translating TOCL properties into LTL implies dealing with three main aspects. First, temporal operators have to be translated. It is straightforward since LTL provides the same kind of temporal operators. Second, OCL operators have to be evaluated, especially those that query the model. Finally, queries specified in XSPEM QDMM have to be translated into LTL. The meaning of these queries obviously depends on the previously defined translational semantics. Thus, while defining the translational semantics, the expert has to provide an ATL library called XSPEM_queries which defines all queries as helpers returning the appropriate LTL formulae. For example the WorkDefinition query "isFinished" becomes a helper "isFinished()" that computes the string "self.name + "_finished"". The TOCL to LTL transformation is composed of two stages (top of Figure 8). The first one is a Higher Order Transformation (HOT-TOCL2ATL) that takes as input a TOCL model and generates an ATL transformation, which is the one executed in the second stage. This second transformation takes as input the business domain model (conforming to XSPEM in our case study) upon which the TOCL properties handled in the first stage will be run to generate the LTL formulae.

This transformation strategy results from two points. First, it is not possible in the first stage to use the XSPEM_queries library because there is no reflexivity in ATL. Second, several TOCL operators can be applied to each element of an input model. So, ATL iteration rules must be generated to traverse the model.

One strong point is that the TOCL2LTL transformation is generic and automated. It is only parametrized by the ATL module that provides the definitions of the QDMM queries for the target formal property language.

Applied on P 1 , the TOCL2LTL transformation produces the following formulae:

<> (Designing _ finished /\ Programming _ finished /\ TestCaseWriting _ finished /\ Documenting _ finished)
The formulae corresponding to P 2 is:

[] ([] (-Designing _ finished) U Documenting _ finished)

Checking SE-LTL properties

Once the Petri net and the LTL formulae have been generated, the selt model checker from the TINA toolbox is used. If the LTL formulae does not hold, it exhibits a counterexample: a specific execution of the model that leads to a state where the property is not satisfied. Figure 9 shows the counter example corresponding to the P 1 property. It consists of a sequence of states. A state is a snapshot of the model showing the places marking. After each state, there is the transition fired to go to the next state. The example shows a deadlock (last transition).

Designer dedicated feedback

Model verification based on a translational semantics provides a significant advantage: the reuse of existing sophisticated model checkers. But there is one significant drawback: results are obtained at the verification level and have to be translated back to the business domain level. This section explains this translation.

Generating PETRINET scenario and trace Using xText, we analyze the output of the selt model-checker and produce a PETRINET scenario and trace using the PETRINET metamodels and the TM3 presented in sections 3 and 4. The PETRINET scenario corresponding to the counter example of Figure 9 is the following:

FireTransitionEvent Designing _ start FireTransitionEvent Designing _ finish FireTransitionEvent Documenting _ start FireTransitionEvent Documenting _ finish FireTransitionEvent TestCaseWriting _ start
The same tool builds the PETRINET models that corresponds to each states of the counter-example (not shown here). Translating PETRINET scenario to XSPEM scenario The PETRINET Scenario is then transformed to an XSPEM Scenario. This transformation converts transition firing events FireTransitionEvent to XSPEM events, either start (StartWD) or finish (Fin-ishWD) a WorkDefinition. The naming convention defined in the Mapping ATL library are used to decode the fired transition names and produce the corresponding XSPEM events and their target workdefinitions. The XSPEM scenario corresponding to the previous PETRINET one is the one shown on Figure 3. To obtain the scenario of figure 4, we have to check the negation of P 2 because we want an example which satisfies P 2 .

The same approach is used to translate PETRINET snapshots (that is PETRINET models) into XSPEM ones. It could be automated using traceability data between source and target models generated during the transformation step.

Related Work

Translational semantics with feedbacking verification results: The main advantage of translational semantics is the reuse of existing tools from the target technical space like model-checkers. Its major issue is that it provides results in the target space that must be translated back to the business domain space. Hegedus et al. [START_REF] Hegedus | Back-annotation of simulation traces with change-driven model transformations[END_REF] propose a method based on a traceability mechanism of model transformations. It relies on a relation between elements of the source (BPEL) and the target (PETRINET) metamodel, implemented by means of annotations in the transformation's source code. The authors propose a technique for the back-annotation of simulation traces based on change-driven model transformations from traces generated by SPIN model checker to the specific animator named BPEL Animation Controller. However, in our approach, we try to generate a scenario (a set of events) that will be animated by a generic animator. In [START_REF] Shah | From UML to Alloy and back again[END_REF], authors use traceability links of the transformation which generates Alloy models from UML. The back-annotation transformation is automatically generated based on these traceability links using a QVT-based implementation. Here, the back-annotation is supported for static model instances, and not for execution traces of Executable DSML models like in our case. In [START_REF] Barber | Providing early feedback in the development cycle through automated application of model checking to software architectures[END_REF], the authors define an approach named Arcade that uses SPIN model checker for evaluating safety and liveness properties of a Domain Reference Architecture that is translated to Promela language. Arcade interprets SPIN counter-example and generates an Architecture Trace Diagram (ATD) that has two dimensions: a vertical dimension that represents time and a horizontal dimension representing SPIN processes. Nevertheless, they do not define a high-level abstraction between business level and formal level. Contrary to our work, we separate the two domains (DSML and formal verification ones) and we hide all formal aspects by translating formal results to business ones. In [START_REF] Pelliccione | Charmy: A framework for designing and verifying architectural specifications[END_REF], Pelliccione et al. present a software tool platform for the model-based design and validation of software architectures, named CHARMY, that offers an extension called SASIM deriving from Theseus approach [START_REF] Goldsby | A visualization framework for the modeling and formal analysis of high assurance systems[END_REF]. Both translate the violation trace from SPIN model checker on a generated sequence diagram and an animated UML state diagrams. vUML [START_REF] Lilius | vuml: a tool for verifying uml models[END_REF] also use the same approach. CHARMY, Theseus and vUML are based on a very ad hoc approach that uses UML diagrams and SPIN model checker. On the contrary, we rely on a generic approach that can be applied to other DSML.

All the above approaches aim at verifying a specific DSML through formal tools by translating business semantics into formal one and by feedbacking formal verification results to the initial business level. However, our work provides a generic approach for the verification of executable DSML. It is based on the explicit definition of the different concerns involved in model execution (runtime information expressed as queries, events) thanks to the executable DSML metamodelling pattern. Based on this pattern and the translation semantics, generic transformations allow to translate user properties to logical formulae and verification results back to business level.

Behavioral property Patterns: To verify BPEL service composition schemas, [START_REF] Yu | Pattern Based Property Specification and Verification for Service Composition[END_REF] proposes a property specification language based on ontologies and named PROPOLS which allows composition of the patterns defined in [START_REF] Dwyer | Property specification patterns for finite-state verification[END_REF]. These patterns are close to TOCL temporal operators and composition corresponds to OCL operators. Rather than relying on a Query Definition MetaModel, a one to one mapping has to be defined for each property item to the corresponding BPEL operation.

In [START_REF] Autili | Graphical scenarios for specifying temporal properties: an automated approach[END_REF], the authors provide a graphical tool named PSC (Property Sequence Chart), to specify temporal properties as an extended notation of a selected subset of the UML 2.0 Interaction Sequence Diagrams. Theseus approach [START_REF] Goldsby | A visualization framework for the modeling and formal analysis of high assurance systems[END_REF] uses SPIDER to translate natural language properties to the property specification language of the targeted analysis tools. Both, PSC and SPIDER are specific approaches used in a specific domain between UML and SPIN. User-oriented property languages, graphical or not, are an important point to make formal verification accessible to end users. TOCL is certainly not the best-suited language despite it is an extension of OCL, a well accepted lan-guage in MDE. Nevertheless, we consider it can be used as a pivot language for more user-oriented languages.

Conclusion

Using the XSPEM case study, this paper has illustrated a method to ease the integration of verification tools for safety and liveness properties on executable models. It relies on the executable DSML metamodeling pattern using a translation to the Time Petri nets as formal verification language providing the semantics. This could be applied to any other kind of formal language providing automated verification tools. We have recently applied it to the FIACRE intermediate verification language [START_REF] Berthomieu | FI-ACRE: an Intermediate Language for Model Verification in the TOPCASED Environment[END_REF], that abstracts several existing verification toolsets such as TINA and CADP in order to factorize common aspects and avoid redefining transformations for all toolsets. This experiment results will be presented in a forthcoming paper. The integration is provided through QDMM extension to the pattern and automated translations on the property side.

This approach has been designed for domain specific languages and this is a key point to keep it simple. It is currently being experimented for several significantly different DSMLs (and sub-languages from general purpose languages) such as data flow models, SAE AADL, SDL, UML and SYSML class, state machine, activity and composite structure diagrams. But, it is still to be shown if it can scale up to more complex languages or to languages that combine different models of computation.

These preliminary experiments allowed a first validation of our proposal for the systematic construction of verification tools for behavioral properties expressed on a DSML. We have chosen to rely on TOCL to express properties at the business domain level because it is close to OCL. However, some early feedback have shown that it is still not well suited to many end users. Therefore, we might need to investigate new useroriented language for expressing behavioral constraints. It is a problem that has already been identified in [START_REF] Ljungkrantz | A formal specification language for plc-based control logic[END_REF]: the authors have defined a new dialect of linear temporal logic more suitable for control engineers.

Finally, we propose to ease the feedback of verification results. We currently rely on naming conventions. We are investigating the explicit construction of the links between the business domain and verification models elements during the downward translation so that they can be used during the upward feedback.

Fig. 3 :Fig. 4 :

 34 Fig. 3: A scenario from P 1

1 Fig. 5 :

 15 Fig. 5: XSPEM Metamodels

  Fig. 6: PETRINET metamodels

Fig. 7 :

 7 Fig. 7: Generated Petri net for parts of Figure 2 5.1 Translation XSPEM 2petrinet

Fig. 8 :

 8 Fig. 8: Architecture of the ATL transformations

FALSE state 0 :Fig. 9 :

 09 Fig.9: Selt output for P 1 checked on example of Figure2

Fig. 10 :

 10 Fig. 10: General approach illustrated with the XSPEM case study

Verification step back results model checker Transformations level formal property Business level Formal level business constraint business model business results formal model formal results

Toolkit In OPen source for Critical Applications & SystEms Development, www.topcased.org

This work was funded by the French ministries of Industry and Research and the Midi-Pyrénées regional authorities through the FUI TOPCASED, ANR OpenEmbedd, ITEA SPICES and ITEA2 OPEES projects

[START_REF] Berthomieu | The tool TINA -Construction of Abstract State Spaces for Petri Nets and Time Petri Nets[END_REF] Specification and Description Language: is an object-oriented formal language developed and standardized by The International Telecommunication Standardization Sector (ITU-T)

Ecore is the metalanguage of Eclipse Modeling Framework, www.eclipse.org/modeling/emf

A read-arc only checks that there is enough tokens in the input place but those tokens are not withdrawn when the transition is fired.

http://combemale.svn.enseeiht.fr/proto/fr.irit.acadie.xspem2tina/


Nicolas Maillet
email: nicolas.maillet@inria.fr

Background: Nowadays, metagenomic sample analyses are mainly achieved by comparing them with a priori knowledge stored in data banks. Even if powerful, such approaches do not allow to exploit unknown and/or "unculturable" species, for instance estimated at 99% for Bacteria.

Methods: This work introduces Compareads, a de novo comparative metagenomic approach that returns the reads that are similar between two possibly metagenomic datasets generated by High Throughput Sequencers.

One originality of this work consists in its ability to deal with huge datasets. The second main contribution presented in this paper is the design of a probabilistic data structure based on Bloom filters enabling to index millions of reads with a limited memory footprint and a controlled error rate.

Results: We show that Compareads enables to retrieve biological information while being able to scale to huge datasets. Its time and memory features make Compareads usable on read sets each composed of more than 100 million Illumina reads in a few hours and consuming 4Gb of memory, and thus usable on today's personal computers.

Conclusion:

Using a new data structure, Compareads is a practical solution for comparing de novo huge metagenomic samples.

Introduction

The past five years have seen the arrival of High Throughput Sequencing (HTS), also known as Next-Generation Sequencing (NGS). These technologies drastically lowered sequencing costs and increased sequencing throughput. They radically changed molecular biology and computational biology, as data generation is no longer a bottleneck. In fact, nowadays a major challenge is the analysis and interpretation of sequencing data [START_REF] Desai | From genomics to metagenomics[END_REF]. HTS democratized access to sequencing to almost all biological labs over the world. It also opened the doors to new techniques such as ChipSeq [START_REF] Johnson | Genome-wide mapping of in vivo protein-DNA interactions[END_REF], ClipSeq [START_REF] Licatalosi | HITS-CLIP yields genome-wide insights into brain alternative RNA processing[END_REF], RadSeq [START_REF] Davey | RADSeq: next-generation population genetics[END_REF] and the topic of this work, metagenomics [START_REF] Wooley | A Primer on Metagenomics[END_REF].

Metagenomics, also known as "environmental genomics", provides an alternative to traditional singlegenome studies for exploring the microbial world. Most microorganisms (up to 99% of Bacteria [START_REF] Amann | Phylogenetic identification and in situ detection of individual microbial cells without cultivation[END_REF]) are unknown and possibly "unculturable". Even if traditional genomics sequencing methods are well studied, they are not suited for environmental samples, because of the need to cultivate clones. By sequencing uncultured genomes directly from environmental samples, metagenomics offers new ways to study this unexplored diversity.

HTS technologies provide fragments of sequences (called reads) of length a few hundred base pairs without any information about the locus nor the orientation on the molecule they come from. In the metagenomic context, an additional difficulty comes from the fact that each read may belong to any species.

Nowadays, it is difficult to assemble complex metagenomes (such as soil or water metagenomes) into longer consensus sequences, because reads from different species may be merged into one chimeric sequence.

Mende and colleagues [START_REF] Mende | Assessment of Metagenomic Assembly Using Simulated Next Generation Sequencing Data[END_REF] showed that for a 400-genomes metagenome, using simulated Illumina reads, 37% of the assembled sequences were chimeric. Thus currently, reads from metagenomes are used to estimate the biodiversity [START_REF] Wang | MetaCluster 4.0: a novel binning algorithm for NGS reads and huge number of species[END_REF] or may be compared to known databases, providing information w.r.t. the current scientific knowledge [START_REF] Markowitz | IMG/M: the integrated metagenome data management and comparative analysis system[END_REF][START_REF] Huson | MEGAN analysis of metagenomic data[END_REF]. Another way to exploit two or more metagenomic datasets is to compare them together, enabling to understand how genomic differences are related to environmental ones (biotopes localizations and/or time spent after an event).

Comparative metagenomics can deal with many aspects, such as sequence composition, i.e. GC content [START_REF] Foerstner | Environments shape the nucleotide composition of genomes[END_REF], and genome size [START_REF] Raes | Prediction of effective genome size in metagenomic samples[END_REF], taxonomic diversity [START_REF] Jaenicke | Comparative and joint analysis of two metagenomic datasets from a biogas fermenter obtained by 454-pyrosequencing[END_REF], functional content [START_REF] Sommer | Functional characterization of the antibiotic resistance reservoir in the human microflora[END_REF], etc. Several methods are currently developed for comparative metagenomics analyses. Some are based on statistical methods to deal with such large number of descriptive variables, e.g. principal component analysis (PCA).

To the best of our knowledge, there is no software designed to compare two or more metagenomic samples at the read level, that is identifying reads that are shared or similar between samples. This can be simply used to compute a similarity measure between samples such as the number or purcentage of similar reads between pairs of samples. When dealing with more than two samples, this would enable among others to classify metagenomics samples based on their raw reads content. One could use the popular tool Blast to align reads in an all-vs-all way, however it is not designed specifically to this task and more importantly it can not cope in time and memory with the size of nowadays metagenomic samples obtained with current sequencing technologies. For instance, with the aim of exploring the diversity of small eukaryotes in the oceans all over the world, the expedition "Tara Ocean" [START_REF] Karsenti | Towards an 'Oceans Systems Biology[END_REF] is generating more than 400 metagenomic samples containing each around 100 million short reads, that will need to be compared to each other.

Here, we introduce a time and memory-efficient method for extracting similar reads between two metagenomic datasets. The similarity is based on shared k-mers (words of length k). In order to fit with current memory capacities, the data structure we use is a modified version of a Bloom filter [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF]. Bloom filters have recently been used in bioinformatics, notably for assembly graph partitioning [START_REF] Pell | Scaling metagenome sequence assembly with probabilistic de Bruijn graphs[END_REF], which enabled to perform metagenomic de novo assembly using 30x less memory. This manuscript presents two main contributions: (I) a new algorithm, called Compareads, which computes the similarity measure between two metagenomics datasets; (II) a new simple but extremely efficient data structure based on the Bloom filter for storing the presence/absence of k-mers in huge datasets. The manuscript is organized as follows: in Section 1, we depict the Compareads algorithm and the new data structure. In Section 2 we provide results both about the data structure and about Compareads, showing the efficiency of our approach in term of computation time, memory and biological accuracy.

Methods

Preliminaries and definitions A sequence is composed by zero or more symbols from an alphabet Σ. In this work, as we are dealing with DNA, Σ = {A, C, G, T }. A sequence s of length n on Σ is denoted also by

s[0]s[1] . . . s[n -1], where s[i] ∈ Σ for 0 ≤ i < n. We denote by s[i, j] the substring s[i]s[i + 1] . . . s[j] of s.
In this case, we say that the substring s[i, j] occurs at position i in s. We call k-mer a sequence of length k, and s[i, i + k -1] is a k-mer occurring at position i in s.

Overview of Compareads Compareads is designed for finding similar sequences between two read sets. This basic operation may appear extremely simple. However, it has to be highly efficient, in term of computation time and memory footprint, in order to scale with huge metagenomics datasets.

In order to perform efficiently this operation, Compareads indexes k-mers and uses a rough but efficient notion of "similar sequences" defined as follows:

Definition 1 (shared k-mer) Two sequences s 1 and s 2 share a k-mer if and only if ∃(i 1 , i 2 ) such that Limiting the indexation space To control the approximation error (see Section 1.4), the indexing phase is interrupted whenever the volume of k-mers in the first reads of B exceeds a fixed value n. The query phase is then performed on the whole A dataset. This phase returns a partial intersection between A and a first chunk of reads from B. The remaining partial intersections between A and the next chunks of reads from B (each representing a volume of n k-mers or less) are sequentially computed, until all the reads from B have been indexed. Eventually, Compareads returns the union of all partial intersections. Note that, in terms of results, this partitioning approach is strictly equivalent to performing a complete indexing of B then a query of all the reads from A. To avoid redundant computations, reads from A considered as "similar" in one of the partial intersections are tagged using a bitvector and are not queried further.

s 1 [i 1 , i 1 + k -1] = s 2 [i 2 , i 2 + k -1].

Time complexity

Let n A and n B be the number of k-mers respectively in set A and set B.

Computing (A ∩ B) is done in time O(n B ) (indexation) + O n A × n B n (query). The n B
n term is due to the limitation of the indexation space. In practice, for instance for a classical Illumina read set, this term is bellow 10.

Ad hoc data structure

The index data structure we use is based on a Bloom filter, specially designed for the task of storing efficiently a huge set of k-mers, while being fast to build and to query. We shortly recall in Section 1.2.1 what a Bloom filter is before describing, in Section 1.2.2, our data structure called BDS.

Bloom filter

A Bloom filter is a probabilistic data structure designed to test the membership of elements in a set [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF].

It consists of an array of m bits, all initialized to zero, and a set of hash functions. Each hash function maps an element to a single position in the array. Each element is associated, through the values of the hash functions, to several positions in the array. To insert an element in the structure, the bits in the array associated to this element are all set to one. The structure answers membership queries by checking whether all the bits in the array associated to an element are set to one. This data structure is probabilistic in nature, as false positives are possible. Even if an element is not in the set, its bits in the array may still be all set to one. This is because the bits associated to an element may independently be associated to other elements. Hence, the Bloom filter returns a wrong answer with non-zero probability. This probability is the false positive rate. An asymptotic approximation of the false positive rate is 0.6185 m/n , assuming n elements are inserted in the m-bits array, and (ln 2 • (m/n)) hash functions are used [START_REF] Broder | Network applications of bloom filters: A survey[END_REF]. False negatives never occur: if an element belongs to the set, the Bloom filter always answers positively. Bloom filters are space-efficient: only (n log 2 e • log 2 (1/ǫ)) bits are required to support membership queries for n elements with a false positive rate of ǫ [START_REF] Broder | Network applications of bloom filters: A survey[END_REF].

The Bloom Data Structure index

In this article, we consider a slightly different variation of Bloom filters: instead of using a single array of bits, each hash function corresponds to a distinct array, disjoint from all other functions. In terms of performance, with uniform hash functions, this variation is asymptotically equivalent to the original definition [START_REF] Broder | Network applications of bloom filters: A survey[END_REF]. To avoid confusion with classical Bloom filters, we refer to this variation as BDS, standing for Bloom Data Structure.

Particular hash functions The hash functions used in this framework is a specific family of functions, which can be efficiently computed on consecutive k-mers. We consider the set of functions which map a k-mer to a bit sequence of length k, where each nucleotide is associated to a bit set to 0 or 1, depending only on its type (A, C, G or T). An exhaustive enumeration, in equations 1 and 2, shows that there exists only 7 functions in this set. We can distinguished two types, the first three, f 1 , f 2 and f 3 , are said to be balanced (equation 1), whereas the other four are said to be unbalanced (equation 2)

f j : Σ k → {0, 1} k : ∀i ∈ [1, k]    f 1 (s)[i] = 0 if s[i] = A or C f 1 (s)[i] = 1 otherwise f 2 (s)[i] = 0 if s[i] = A or G f 2 (s)[i] = 1 otherwise f 3 (s)[i] = 0 if s[i] = A or T f 3 (s)[i] = 1 otherwise (1)
f j : Σ k → {0, 1} k : ∀i ∈ [1, k]        f 4 (s)[i] = 0 if s[i] = A f 4 (s)[i] = 1 otherwise f 5 (s)[i] = 0 if s[i] = C f 5 (s)[i] = 1 otherwise f 6 (s)[i] = 0 if s[i] = G f 6 (s)[i] = 1 otherwise f 7 (s)[i] = 0 if s[i] = T f 7 (s)[i] = 1 otherwise (2)
One important property of these functions is that there is a simple relationship between the hash values of two consecutive k-mers in a read. One can see that the hash value of the next k-mer can be quickly computed, by left-shifting the binary sequence of the previous hash value and appending an extra bit. These functions are not classical hash functions, yet we show that they exhibit good hashing properties when applied to k-mers. In Section 2.1, the performance of these functions is compared with that of a classical hash function in terms of computation time, and false positive rate in the BDS.

The Compareads pipeline

Computing (A ∩ B) is asymmetrical. Indeed (A ∩ B) does not contain the reads from B which are similar to reads in A. For doing this, one needs to compute also (B ∩ A). In practice, for fully and symmetrically comparing two sets A and B we apply a pipeline slightly more complicated than simply (A ∩ B) followed by (B ∩ A). This whole pipeline, designed for reducing a heuristic effect is described in Section 1.4.1.

Similarity measure While comparing read sets A and B, the result provided by Compareads is composed of two sets: (A ∩ B) and (B ∩ A). Then, a similarity measure between the two datasets is computed as follows:

Sim(A, B) = ˛A ∩ B ˛+˛B ∩ A
|A|+|B| * 100, where |X| denotes the cardinality of the set X.

Dealing with false positives

Our approach may generate false positives for two reasons we describe in the two upcoming sections which also expose solutions for limiting these effects.

False positives due to k-mer shared between a read and a dataset

Using t > 1, Compareads algorithm can call similar sequences that do not respect strictly the definition of similarity given in definition 2. Indeed, steps described in Section 1.1 detect reads from A that share at least t k-mers with reads from B. This is less stringent than finding reads from A that share at least t k-mers with at least one read from set B. In fact, the t k-mers found in read A are possibly spread over two or more distinct reads from set B.

This issue can be mitigated by performing the following steps to compute both (A ∩ B) and (B ∩ A): However, the last read from set A is a case of false positive. It contains k-mers spread over distinct reads from B, the latter belonging to (B ∩ A). Thus, even during step 3, these two k-mers remain shared with reads from set (B ∩ A) and are output in (A ∩ B).

1. Compute (A ∩ B),
Note that in practice, the last set (A ∩ B) is obtained by computing (A ∩ B) * ∩ (B ∩ A) instead of simply A ∩ (B ∩ A) (used here for simplifying the reading). This operation provides the same result but is computed faster as

|(A ∩ B) * | ≤ |A|.
As outlined in the example Figure 1, this pipeline may still conserve some false positives. The latter are characterized by the fact that they contain t shared k-mers with at least two distinct reads from the indexed dataset B themselves considered as similar to reads of set A. Even if this side effect is difficult to assess, we show in Section 2.2 that Compareads provides trustfull results, highly similar to a classical approach, on several real datasets.

Bloom filter false positives

As exposed in Section 1.2.2, the BDS index is a probabilistic data structure, that may consider a k-mer as indexed while this is not the case (i.e. a false positive or FP). Here, we analysed the variations of the false positive rate for each hash function and their combinations with respect to the parameter k and the number n of distinct indexed k-mers. This enabled then to choose the parameters and the appropriate combination of functions that give the best tradeoff between memory and false positive rate.

FP probablity for each function Assuming the nucleotide composition of the indexed k-mers and of the query k-mers are unbiaised, we can easily compute the probability, P F P (f i , k, n), for any query k-mer to be a false positive with one of the seven hash functions, f i (see the appendix Section A for details). The expressions of this probablity are presented in equations 3 for a balanced hash function and 4 for an unbalanced one.

∀i ∈ {1, 2, 3}

P F P (f i , k, n) = 1 -(1 - 1 2 k ) n (3)
∀i ∈ {4, 5, 6, 7}

P F P (f i , k, n) = k x=0 k x a x (1 -(1 -a x ) n ) with a x = ( 1 4 ) x ( 3 4 ) k-x (4)
We have plotted in Figure 2a the theoretical FP rate for both types of hash function, and we can see that balanced functions give much less false positives than unbalanced ones. This is due to the fact that balanced functions distributes the hash codes uniformally over the 2 k bit-array, while this not the case for the unbalanced ones.

FP probablity for a combination of functions One important property of the balanced hash functions is that there do not exist two distinct k-mers that have the same couple of hash codes with any two of these functions. This implies that the probability of having a FP with one function does not depend on the result with another function (or almost not, see details in appendix Section A). The probability of FP can then be easily computed as follows:

P F P (f 1 ∩ f 2 ∩ f 3 , k, n) (1 -(1 - 1 2 k ) n ) 3 (5)
This "independence" property implies also that combining these 3 functions in our BDS is very efficient to reduce the FP rate, as can be seen in Figure 2a especially for large values of n

Concerning the unbalanced functions, this independence property is lost, since it is possible to find couples of distinct k-mers that share the same couple of hash codes for at least 2 of the unbalanced functions, or for one balanced function and at least one unbalanced. Therefore we could not figure out the theoretical FP rate. Nevertheless, we estimated it by simulations. We found that empirical results are very close to the formula obtained by multiplying the individual probabilities, i.e. assuming complete independence between all functions (Figure 2a). For details about how empirical results were obtained, see appendix Section B.

Choice of parameters The comparison of these FP rate curves led us to choose the combination of the three balanced functions plus an unbalanced one. This choice is motivated by the fact that unbalanced functions are not essential, as they have a limited effect on the FP rate (Figure 2a). Moreover, using these functions doubles the memory cost as storing four or less functions needs 2 k-1 bytes, while seven functions need 2 k bytes. Thus we choose to consume 2 k-1 bytes memory, having space for four functions, naturally composed of the thee balanced plus one of the unbalanced functions.

For this combination of functions, we plotted the FP rate as a function of n and for several values of k in Figure 2b. The larger k is, the less FP we get for a given number of indexed k-mers. Consequently, for large values of k, more k-mers can be indexed while maintaining a reasonable FP rate. However, the memory allocated to BDS grows with k and larger values of k increases the stringency of our similarity measure. We can see in Figure 2b, that using k-mers of size at least 30 enables to index at least 300 millions of k-mers with less than 2% of false positives.

We used k = 33 and when indexing up to n =one billion distinct k-mers we obtain a theoretical upper bound of 0.13% of false positives (with 3 balanced functions, equation 5). The FP rate is even lower when adding one of the unbalanced function, we estimated it empirically to 0.114%.

Results

Practical performance of the BDS, comparison with other data structures

Since hash functions described in Section 1.2.2 have a fixed range, the memory used by the BDS depends only on the value of k and the number of hash functions used. Recall that each hash function is associated to a dedicated bit array which occupies 2 k bits. Using 7 hash functions, the BDS has a total memory footprint of 7 * 2 k bits and can be stored in 2 k bytes. When using only 4 hash functions, the BDS occupies 4 * 2 k bits and can be stored in 2 k-1 bytes. As shown in Section 1.4.2, using 4 hash functions and k = 33 is a good set of parameters for indexing n = one billion distinct k-mers. Using such parameters, the memory usage of Compareads is 4 GB.

We propose here a comparative analysis of the BDS with other data structures. In next Section we show that classical non probabilistic data structures result in a worst time and memory achievements, while in Section 2.1.2, we show that the BDS is the best suited for the problem of indexing huge amounts of k-mers.

Comparison with non probabilistic data structures: suffix array and hash table

Indexing n characters using the simplest version of a suffix array (not enhanced [START_REF] Abouelhoda | Replacing suffix trees with enhanced suffix arrays[END_REF] and without LCP information) requires 5n bytes of memory [START_REF] Vyverman | Prospects and limitations of full-text index structures in genome analysis[END_REF]. Compared to our set of parameters where n = 1 billion, the memory footprint would be 5 × 10 9 bytes, i.e. 4.66 GB. While this is comparable to the BDS, the query time of the suffix array, O(k log n), is significantly worse.

An hash table can be used to store an exact set of k-mers. Such structure stores the k-mers explicitly, between two samples using Compareads is more than 30 times faster than using Blast for a close total number of similar reads.

For each experimentation, samples were hierarchically clustered based on their pairwise similarity scores and then drawn as a dendrogram. As shown in Figure 4, the dendrogram obtained with the Blast approach (a) is slightly different but the three main branches are the same than with the Compareads approach (b). Interestingly, these branches discriminate three groups of samples corresponding to the three different based on Compareads approach (b, c, d) show a similar organization. Increasing the number of shared k-mers leads to be more stringent and decreases the number of similar reads but do not affect the global organization of the dendrogram, demonstrating the robustness of our similarity measure.

Applying Compareads to Global Ocean metagenomic samples

We tested Compareads on a larger and famous public dataset from the Global Ocean Sampling (The Sorcerer II expedition) [START_REF] Rusch | The Sorcerer II Global Ocean Sampling expedition: northwest Atlantic through eastern tropical Pacific[END_REF]. It is composed of 44 samples from the microbial world of seawater, collected across several thousand of kilometers from the Northwest Atlantic through the Eastern Tropical Pacific oceans and for which an analysis of similarity between samples has been done [START_REF] Rusch | The Sorcerer II Global Ocean Sampling expedition: northwest Atlantic through eastern tropical Pacific[END_REF]. The whole dataset is composed of 44 samples containing each on average 174759 long reads (1249 nucleotides per read on average, Sanger technology).

Compareads compute all of the 990 intersections in 72 hours and half: on average, one intersection was performed in 4 minutes and 23 seconds on a single core of an Intel➤ Xeon➤ CPU X5550 at 2.67GHz.

Results presented in Figure 5 are highly similar to those presented in the original publication [START_REF] Rusch | The Sorcerer II Global Ocean Sampling expedition: northwest Atlantic through eastern tropical Pacific[END_REF], p.418.

Two main groups are well discriminated. The first one, represented in turquoise-blue, groups together almost all samples coming from temperate seawater of the North American East Coast except the 14 one, like in the original study. This group also contains two samples really different from all others: the first contains freshwater and the second hypersaline water. The dark-green part corresponds to samples coming from the north part while light-green one gathers samples from the south part. Orange samples correspond to estuary.

All of those three groups are identical to the original study. The second main part, colored in yellow, groups together datasets of tropical and Sargasso seawater. The dark-blue part aggregates samples coming only from Galapagos Islands. Red square delimitates Sargasso Sea samples. On the original study, the sample 00a is not in this group. According to metadata, the gray part, like in the original publication, is composed of various samples. Finally, purple samples regroup both Caribbean Sea and some Open Ocean datasets, as Those results show that Compareads can also be used on Sanger reads and deliver reliable biological conclusions. Indeed, despite of false positives and the simple definition of similarity, we were able to retrieve the classification of metagenomes according to their geographical origin.

Conclusion

Motivated by de novo comparative metagenomics, this paper proposes two main contributions. The first one is a data structure based on Bloom filters that can index, for instance up to one billion distinct words of length 33 (33-mers) using 4Gb of memory, with an error rate of 0.11%, and that is faster to build and request, to the best of our knowledge, that any other existing data structure. The second main contribution is a software, called Compareads which uses this data structure to efficiently perform de novo intensive comparisons of huge metagenomic datasets generated by High Throughput Sequencers. We have shown that this approach enables to retrieve and classify differences in species content between metagenomic samples.

For this kind of comparison, our approach is much faster than alternative ones such as Blast and thus enables to scale to huge datasets. For instance, we tested the scalability of Compareads on a large oceanic unpublished dataset, from the Tara Ocean expedition [START_REF] Karsenti | Towards an 'Oceans Systems Biology[END_REF]; it is composed of 31 metagenomes and contains overall 3.5 billions of Illumina short reads (108bp). Each intersection was performed in 10 hours and 55 minutes in average using 4Gb of memory. Such features enabled us to compute the 31 * 32 2 = 496 metagenome datasets intersections in 6 days and 10 hours using 50 cores of Intel➤ Xeon➤ CPU X5550 at 2.67GHz.

This would have been unfeasible with any other known existing tools (based on results 2.2, Blast is about 30 times longer and would take more than 6 months to complete this task with the same resources).

Compareads has been conceived for being parallelizable both at fine and coarse grained levels. Future work will consist in implementing a Compareads parallel version exploiting multi-core and GPU chips.

Compareads is released under the CeCILL license and can be freely downloaded from http://alcovna. genouest.org/compareads/.

Author's contributions

  storing the results in a set denoted by (A ∩ B) * . 2. Compute B ∩ (A ∩ B) * storing the results in a set denoted by (B ∩ A).

3 .

 3 Compute A ∩ (B ∩ A) storing the results in a set denoted by (A ∩ B).In a few words, the two output datasets (B ∩ A) and (A ∩ B) are obtained by applying the fundamental operation ( ∩ ) between a query and a read set being itself already the result of the asymmetrical ( ∩ ) operation. This enables to remove some false positives due to k-mers spread over several reads.

Figure 1 :

 1 Figure 1: Representation of the three steps while comparing symmetrically read sets A and B. In each set, reads are represented by horizontal lines. On each read one or two shared k-mers are represented by rectangles.

  (a) FP rate for different combinations of hash functions (b) FP rate for different values of k

Figure 2 :

 2 Figure 2: FP rate as a function of the number of indexed k-mers (in log scale). Plain lines correspond to theoretical predictions, whereas star points correspond to empirical values obtained with simulations. (a) This figure was obtained for k = 33 for balanced and unbalanced functions and some combinations of them. The combination entitled "4 functions" is composed of the 3 balanced functions plus one unbalanced. (b) For several values of k, this figure was obtained for a combination of 4 hash functions: all three balanced plus one unbalanced.

Figure 3 :

 3 Figure3: Comparison of FP rates between classical hash functions and the functions we used in the BDS. FP rate is plotted as a function of the number of indexed k-mers (in log scale), with k = 33. Plain lines correspond to theoretical predictions for the balanced functions (BDS), whereas star points and triangles correspond to empirical values obtained with simulations using respectively BDS functions and classical hash functions. The combination entitled "BDS 4 functions" is the one chosen for Compareads and is composed of the 3 balanced functions plus one unbalanced.

  biological conditions indicated by 1, 10 and 40 in the samples names: 1 corresponds to addition of Carbon in the water, 10 stands for normal condition and 40 for introduction of Nitrogen. Notably, all dendrograms

Figure 4 :

 4 Figure 4: Representation of hierarchical clustering based on pairwise intersections between all samples using Blast (a) and Compareads (b, c, d).

Figure 5 :

 5 Figure 5: Similarity matrix resulting from the comparison of 44 samples from The Sorcerer II Global Ocean Sampling Expedition using Compareads. Grey levels correspond to similarity levels, intersections with more than 50% of similarity are in black. The two main groups, in turquoise-blue and yellow, correspond respectively to north American east coast and tropical samples.

  DL and PP initiated the work. RC and CL provided expertise about Bloom filters datastructures and their statistical aspects. NM and PP made the implementations. NM, CL, DL and PP performed the experiments. All authors participated to the redaction and approved the final manuscript.

  Definition 2 (Similar sequences) Given integers k and t, two sequences s 1 and s 2 are said similar if and only if they share at least t non overlapping k-mers. Compareads computes (A ∩ B) in two steps. The indexing step consists in storing in memory all k-mers having at least one occurrence in the set B. The query step processes reads from set A one by one. For a read r ∈ A, the index is used to test for each k-mer of r if is present in the set B. If at least t non-overlapping k-mers are returned as present, then the read r is inserted in (A ∩ B). The main practical challenge faced by Compareads is to index the possibly huge volume of k-mers contained in B. The data structure must

	In a few words, given two read sets A and B, the goal of the Compareads algorithm is to find the subset
	of reads from A which are similar to a read in B such set being denoted by (A -→ ∩ B). As it is a heuristic (see
	Section 1.4), our algorithm outputs an over-approximation of set (A -→ ∩ B) denoted by (A ∩ B).
	1.1 Computing (A ∩ B)
	therefore fulfill three criteria: it must be quick to build, have a low memory footprint and be quick to request.
	Section 1.2.2 describes the chosen probabilistic data structure, based on a Bloom filter.

Table 1

 1

	, computing one intersection

Table 1 :

 1 The CPU time per intersection and the global CPU time using a single core of an Intel➤ Xeon➤ CPU X5550 at 2.67GHz. Reads Found corresponds to the total number of similar reads in all the 120 intersections.

Acknowledgements

This work was supported by the french ANR-2010-COSI-004 MAPPI Project. Authors warmly thank O. Jaillon from the Genoscope and P. Vandenkoornhuyse from the Ecobio UMR for providing their biological expertise and metagenomic datasets. Additionally, we thank G. Rizk for its help and comments with the data structure and F. Gauthier for the " ∩ " symbol creation.

hence it requires at least n • ⌈ 2k 8 ⌉ • 8 bits (assuming no overhead), i.e. [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF].5 GB for one billion of 33-mers.

Comparison with other hash functions and with a classical Bloom filter

Time comparison with other hash functions The hash functions defined for BDS were designed with speed in mind. In this paragraph, we compare them with a popular and fast hash function (hashlittle2 from http://burtleburtle.net/bob/c/lookup3.c using Jenkins functions). We simulated 1 millon of 100-bp reads, where each nucleotide is drawn uniformly and independently. To simulate the behavior of computing hashes for the BDS, 4 hash values were computed for each 33-mer. For the hashlittle2 function, we simulated this behavior by computing 4 hashes with 4 different initial values. We recorded the time required to compute the hashes for all the 33-mers present in the reads, averaged over 3 executions. Computing the hash with the hashlittle2 function took 13.1 seconds (5.2 MHashes/s), whereas for the BDS hash functions, the same computation took 1.4 seconds (49.8 MHashes/s). Hence, the BDS hash functions are one order of magnitude faster than a classical set of hash functions.

FP rate comparison with other hash functions We can see in Figure 3 that the FP rate of classical hash functions follows the FP rate of our balanced functions (it follows the equation 5 with the exponent 3 being replaced by the number of functions used). However it diverges with more than 3 functions, as we could not add other balanced functions and we added in place unbalanced ones which have higher FP rates. Even if, for more than three functions, classical hash functions produce less FP, the difference with our BDS structure is small: for 1 billion indexed k-mers, combinations of 4 classical functions give 0.01% FP on average, compared to 0.114%. We chose to have a slithly higher FP rate, but with a significant gain in computing time.

Comparison with a classical Bloom filter A classical Bloom filter requires a fixed amount of memory to index n k-mers. Evaluating the Bloom filter memory using formula from Section 1.2.1 for one billion elements, with a false positive rate of 0.114%, yields 1.8 GB of memory. While this is twice smaller than the BDS, a classical Bloom filter would require classical hash functions. However, as shown above, classical hash functions are an order of magnitude slower to compute.

Comparing with a classical approach using Blast

Our approach is an heuristic based on shared k-mers between reads. Here we compare Compareads with a well-established method, Blast [START_REF] Altschul | Basic local alignment search tool[END_REF] that is based on sequence alignment. The dataset used is composed of 15 bacterial metagenomes obtained from fresh water with three different conditions of Carbon/Nitrogen ratio (unpublished data). On average, each sample is composed of 176409 reads with an average of 400

Appendix

A Theoretical details for the false positive rate

As exposed in Section 1.2.2, the BDS index is a probabilistic data structure, that may consider a k-mer as indexed while this is not the case (i.e. a false positive). Here, we tried to express the false positive rate for each hash function that we defined in Section 1.2.2 and their combinations with respect to the parameter k and the number n of distinct indexed k-mers.

False positive probablity for each function Assuming the base composition of the indexed and query k-mers is unbiaised, we can easily compute the probability, P F P (f, k, n), for any query k-mer to be a false positive with one of the seven hash functions, f . This probability depends on the number of distinct k-mers sharing the same hash code. We can notice that for the balanced functions, f 1 , f 2 and f 3 , each 0 and 1 value can come from exactly 2 distinct nucleotides, thus the number of k-mers sharing the same hash code is the same for all k-mers and equals: 2 k . The probability for 2 k-mers to have distinct hash codes is then 1 -2 k 4 k = 1 -1 2 k , and therefore the probability to have at least one k-mer among the n that are indexed sharing the same hash code is:

Note that this corresponds to the false positive probability of any hash function that distributes the hash codes uniformly in a 2 k bit-array, such as those inspired of Jenkins functions, used as a comparison in Section 2.1.2.

As for the unbalanced functions, since the 0 bit-value encodes only one base, the number of k-mers sharing the same hash code depends on the number of 0 in the hash code of the query. For a given query k-mer with a hash code having x 0 the above probability for functions f 4 , f 5 , f 6 and

) n . To obtain the probability for any kmer, we have to sum over the different values of x the latter probability weigthed by the probablity for a k-mer hash code to have x 0. The composition of a given base in a k-mer of length k, assuming unbiased nucleotide composition, follows a binomial distribution, thus we get: ∀i ∈ {4, 5, 6, 7}

We can see in Figure 2 that balanced functions give much less false positives than unbalanced ones. This can be explained by the fact that for unbalanced functions, for a given k-mer with a "normal" composition and thus 25% of 0 in its hash-code, there are many more k-mers with the same hash-code than for a balanced function:

FP probablity for a combination of functions When combining several functions in our BDS, in order to have a false positive for a query k-mer all functions must return a false positive. As concerns the balanced function, we can easily see that for a given kmer, we can not find another k-mer that is a false positive simultaneously for any two 2 of these functions. In other words, there do not exist two distinct k-mers that have the same couple of hash codes with any two of these functions. This implies that the probability of having a false positive with one function does not depend on the result with another function, apart from the fact that the effective number of indexed k-mers that can be a false positive (n in equation 3) is reduced: indeed if x k-mers have the same hash code for one function, these k-mers have a null probability of having the same hash code for another function. Note that this effect can be neglicted given that n is very large. Therefore the product of individual probabilities for each balanced function gives the following upper bound:

Concerning the unbalanced functions, this independence property is lost, since it is possible to find a single k-mer that is a false positive for at least 2 of the unbalanced functions, or for one balanced function and at least one unbalanced. Therefore we could not figure out the theoretical false positive rate, or even an upper bound.

B Empirical estimation of false positive rate

For estimating the false positive rate of the BDS filled up with n k-mers, we made the following experiment.

We generated n+100000 distinct random k-mers, used first for filling up the BDS with n k-mers and then for querying the BDS with the 100000 remaining k-mers. All k-mers being distinct, if the BDS answers "yes" while querying the presence of a k-mer, it is a false positive.

The generation of a huge set of x distinct random k-mers (with x < 4 k ) is not trivial. This was done by dividing the space of 4 k k-mers into x non overlapping blocks, and them by picking up a random k-mer into each block. This method enables to uniformly cover the whole k-mer space.


Olivier Morin

Virginia D ' Auria

Claude Fabre

Julien Laurat
email: julien.laurat@upmc.fr

A high-fidelity single-photon source based on a type-II optical parametric oscillator

Using a continuous-wave type-II optical parametric oscillator below threshold, we have demonstrated a novel source of heralded single-photons with high-fidelity. The generated state is characterized by homodyne detection and exhibits a 79% fidelity with a single-photon Fock state (91% after correction of detection loss). The low admixture of vacuum and the well-defined spatiotemporal mode are critical requirements for their subsequent use in quantum information processing.

The reliable generation of single-photon states is a central resource for the development of quantum information sciences and technologies, including quantum communication and computing [1,2]. For instance, since the seminal proposal by Knill, Laflamme and Milburn [3], single-photons are indeed at the heart of linear-optical quantum computation (LOQC) [4]. Practical implementations however require to generate such states with a low admixture of vacuum as efficient LOQC protocols are constrained by loss thresholds [5]. For instance, the best known figure to date, which applies to cluster state computation, is a 1/2 overall loss tolerance [6], i.e. the product of the source fidelity and detector efficiency has to be above this value. This constraint puts a challenging demand on single-photon generation. Additionally, linear-optical processing cannot increase the fidelity of the state, even with multiple imperfect sources [7]. Moreover, gate implementation requires to make singlephotons interfere with a high visibility and thus to be emitted into a well-controlled spatiotemporal mode.

In this letter, we demonstrate the generation of highfidelity single-photons based on a frequency-degenerate type-II optical parametric oscillator (OPO). Our heralded source combines the two important previous properties: a low admixture of vacuum and a very-well defined spatiotemporal mode due to the OPO cavity.

The generation technique consists in using photon pairs emitted in two distinct modes: the detection of a single-photon on one mode heralds the preparation of a single-photon into the other one [8]. This conditional preparation has been widely used so far with a variety of physical systems, including atomic cascade, atomic ensemble [9] and pulsed single-pass parametric downconversion (PDC) [10][11][12][13]. More recently, continuouswave OPOs, where PDC occurs in cavity, have also been considered for this purpose [14].

In this case, two methods can be used to obtain the two distinct modes. The first one is based on tapping a very small fraction of a squeezed vacuum [15], leading intrinsically to a low count rate. The other method consists in using non-degenerate modes emitted into the same spatial mode of the OPO cavity. Using a type-I non-linear interaction, signal and idler have the same polarization: the only degree of freedom is thus the frequency. The two correlated modes can be then two modes with frequencies separated by multiples of the cavity free spectral range, as demonstrated in [16]. However, this configuration requires for the homodyne tomography to use local oscillators with shifted frequencies. The other possibility is to use a type-II interaction: the frequency-degenerate signal and idler modes have orthogonal polarizations and can be easily separated. This technique, which was not implemented so far, leads to a simpler apparatus and enables us to reach a higher fidelity.

The experimental setup is sketched on Fig. 1. A continuous-wave frequency-doubled Nd:YAG laser (Diabolo, Innolight) pumps below threshold a triply-resonant OPO made of a semimonolithic linear cavity [17]. The input mirror is directly coated on one face of a 10 mmlong KTP crystal (Raicol), with an intensity reflection of 95% for the 532 nm pump and high reflection at 1064 nm. The output coupler, with a 38 mm radius of curvature, is highly reflective for the pump and has a reflection of 90% for the infrared. The free spectral range is equal to ∆ = 4.3 GHz. The internal losses at 1064 nm are estimated to be 0.4%, giving a cavity bandwidth γ = 60 MHz. The triple resonance condition is achieved by adjusting the frequency of the laser and the temperature of the crystal. The OPO cavity is locked on the pump resonance by the Pound-Drever-Hall technique.

At the output of the OPO, the orthogonally-polarized signal and idler modes are spatially separated thanks to a polarizing beam splitter. The idler is directed to frequency filtering elements (conditioning path on Fig. 1) to remove non-degenerate modes. Indeed, the output of the OPO cavity contains a great number of modes separated by the free spectral range: each mode at frequency ω 0 +p∆ is entangled with the mode at ω 0 -p∆, where p is an integer. In order to herald the generation of a singlephoton at the central frequency ω 0 , it is thus necessary to reject all these non-degenerate modes. A first filtering stage is based on a commercial interferential filter (Barr Associates), with a bandwidth equal to 0.5 nm. The remaining modes are then filtered by a linear Fabry-Perot cavity. With a length equal to 0.45 mm and a finesse of 1000, this cavity has a free spectral range of 330 GHz and a bandwidth close to 320 MHz. The bandwidth is chosen in order to be larger than the one of the OPO and the spectral range to be larger than the one of the interferential filter. This two-stage filtering enables us to reach a -25 dB rejection (0.3%) of the unwanted modes.

The experiment is conducted on a cyclic fashion: the filtering cavity is locked during 40 ms, then the measurement period starts for 60 ms.The filtering cavity is locked on resonance by the Dither-and-Lock technique. For this purpose, a counter-propagating beam is injected via an optical switch (Thorlabs OSW12-980E). The beam is rejected at the entrance of the cavity via an optical circulator and detected. During the measurement period, the optical switch is inverted, the filtered mode is detected by a superconducting single-photon detector (SSPD, Scontel) working at cryogenic temperature. The quantum efficiency at 1064 nm reaches 7%. Notably, dark noise is negligible (below 1 Hz), an important feature as false events degrade the fidelity of the heralded state [18].

Given a detection event on the conditioning path, the heralded state is then characterized by quantum state tomography performed via a homodyne detection. The 6 mW local oscillator is provided by the 1064 nm output of the laser, after spatial mode filtering by a high-finesse cavity ('mode cleaner' in Fig. 1). The homodyne detection (55 MHz bandwidth) is based on a pair of high quantum efficiency InGaAs photodiodes (Fermionics, 500 µm, 97%). An optical isolator just before the homodyne detection is required in our setup to remove light from the local isolator backscattered by the photodiodes. For each trigger event, the photocurrent is recorded with an oscilloscope (Lecroy Wavepro 7300A) at a sampling rate of 5 GS/s during 100 ns. The local oscillator phase is swept during the measurement.

As the local oscillator is continuous, post-processing is required to extract the conditional state in the temporal mode f (t). This mode depends on the bandwidth of the OPO and of the filtering path. In our case, the filtering path is much wider and does not affect the temporal mode. The theoretical expression of the optimal mode is thus f (t) = √ πγe -πγ|t| [14]. Overall, each trigger event leads to a single outcome of the quadrature measurement, which is obtained as x = f (t)x(t)dt, where x(t) is the measured difference photocurrent signal from the homodyne detection. Accumulated measurements are then processed with a maximum likelihood algorithm [19], which provides the density matrix of the heralded state and the corresponding Wigner function.

The value of γ is optimized to obtain the highest contribution of single-photon, i.e. the largest negativity of the Wigner function. The resulting value (65 MHz) is slightly different from the OPO bandwidth as the gain of the electronics is not flat over the whole spectrum.

Results are displayed on Fig. 2. Figure 2(a) gives the histogram of the measured quadrature values, while fig. 2(b) and 2(c) show respectively the diagonal elements of the reconstructed density matrix and the corresponding Wigner function. The state is in good approximation a mixture of vacuum and single-photon state. The twophoton component is around 3%, resulting from multiphoton pairs creating by the down-conversion process. The OPO is operated far below threshold (1 mW for a threshold of 80 mW) to limit this contribution. Without any corrections, the single-photon component reaches 78.6 ± 0.5% (Error bar is estimated following the method in [20]). To the best of our knowledge, this is the highest value reported to date. By taking into account the detection losses, we infer a value as high as 91%.

Table 1 gives the various efficiencies affecting the final result. Two contributions have to be distinguished. A first one comes from the overall losses coming from propagation and detection. This includes the propagation efficiency η prop , the mode overlap limited by the visibility V of the local oscillator-signal interference η vis = V 2 , the efficiency of the photodiodes η phot and the electronic noise of the detection η noise (the electronic noise is 20 dB below vacuum noise at the central frequency) [21]. They sum up to an overall detection loss of 1 -η tot =15%, which is taken into account to give the values with correction. The second contribution is more fundamental and cannot be corrected as it is associated with the generation process. It depends on the OPO escape efficiency, given by η OPO = T /(T + L), where T is the transmission of the output coupler and L the intracavity loss. This value is estimated to be η OPO = 96%. Given these estimated parameters, the expected vacuum component is 18%, in very good agreement with the measured value. The heralding rate for single-photon generation is 30 kHz. Given the bandwidth of the OPO, it corresponds to a brightness of 400 photons/s per MHz. This rate is mainly limited by the losses in the conditioning path, which reaches 97%: the quantum efficiency of the SSPD is 7% and the overall transmission (including optical switch and filtering elements) is 40%. Corrected from the losses in this path, the rate can be close to 750 kHz.

In conclusion, we have reported the generation of high-fidelity heralded single-photons, using a configuration based on a type-II optical parametric oscillator below threshold. A 79% fidelity has been demonstrated, mainly limited by the losses in the detection path. Thanks to the OPO cavity, the spatial mode enables to reach high interference visibilities without the need of additional filterings. Moreover, the frequency-degenerate interaction makes the operation much simpler than previous realizations. This practical tool can facilitate the implementation of various new experiments in quantum information processing [20]. The technique developed in [22] could be extended to our efficient scheme if time-gating of the single-photons is required. This work is supported by the ERA-NET CHIST-ERA (QScale) and the BQR from UPMC. V. D'Auria acknowledges the support from the EC (Marie Curie program). C. Fabre and J. Laurat are members of the IUF.

Fig. 1 .

 1 Fig. 1. Experimental setup. A continuous-wave type-II optical parametric oscillator (OPO) is pumped far below threshold. The orthogonally-polarized signal and idler modes are separated by a polarizing beam-splitter (PBS). The idler mode is frequency filtered (Conditioning Path) and detected by a superconducting singlephoton detector (SSPD). Given a detection event, the heralded single-photon is characterized by homodyne detection. FP: Fabry-Perot cavity, IF: interferential filter.

Fig. 2 .

 2 Photon number

Table 1 .

 1 Experimental efficiencies

	η noise η phot	η vis	η prop η tot	η OPO
	96%	97% (98%) 2 95% 85%	96%