Elsa Bernard

Laurent Jacob

Julien Mairal

Jean-Philippe Vert

Efficient RNA Isoform Identification and Quantification from RNA-Seq Data with Network Flows

Several state-of-the-art methods for isoform identification and quantification are based on ℓ 1regularized regression, such as the Lasso. However, explicitly listing the-possibly exponentiallylarge set of candidate transcripts is intractable for genes with many exons. For this reason, existing approaches using the ℓ 1 -penalty are either restricted to genes with few exons, or only run the regression algorithm on a small set of pre-selected isoforms.

We introduce a new technique called FlipFlop which can efficiently tackle the sparse estimation problem on the full set of candidate isoforms by using network flow optimization. Our technique removes the need of a preselection step, leading to better isoform identification while keeping a low computational cost. Experiments with synthetic and real RNA-Seq data confirm that our approach is more accurate than alternative methods and one of the fastest available.

Introduction

Over the past decade, quantitation of mRNA molecules in a cell population has become a popular approach to study the effect of several factors on cellular activity. Typical applications include the detection of genes whose expression varies between two or more populations of samples (differential analysis), classification of samples based on gene expression [START_REF] Van't Veer | Gene expression profiling predicts clinical outcome of breast cancers[END_REF], and clustering, which consists of identifying a grouping structure in a sample set [START_REF] Perou | Molecular portraits of human breast tumours[END_REF]. While probebased DNA microarray technologies only allow to quantitate mRNA molecules whose sequence is known in advance, the recent development of deep sequencing has removed this restriction. More precisely, RNA-Seq technologies [START_REF] Mortazavi | Mapping and quantifying mammalian transcriptomes by RNA-Seq[END_REF] allow the sequencing of cDNA molecules obtained by reverse transcription of RNA molecules present in the cell. Consequently, any transcript can be sequenced and therefore quantitated, even though its sequence might not be available a priori for designing a specific probe. In addition to facilitating the study of non-coding parts of known genomes and organisms whose genome has not been sequenced [START_REF] Mortazavi | Scaffolding a caenorhabditis nematode genome with RNA-Seq[END_REF], RNA-Seq technologies facilitate the quantitation of alternatively spliced genes. Genes in eukaryote cells indeed contain a succession of exon and intron sequences. Transcription results in a pre-mRNA molecule from which most introns are removed and some exons are retained during a processing step called RNA splicing. It is estimated that more than 95% of multiexonic genes are subject to alternative splicing [START_REF] Pan | Deep surveying of alternative splicing complexity in the human transcriptome by high-throughput sequencing[END_REF]: the set of exons retained during splicing can vary, resulting for the same gene in different versions of the mRNA, referred to as transcripts or isoforms. Identification and quantification of isoforms present in a sample is of outmost interest because different isoforms can later be translated as different proteins. Detection of isoforms whose presence or quantity varies between samples may lead to new biomarkers and highlight novel biological processes invisible at the gene level.

Sequencing technologies are well suited to transcript quantitation as the read density observed along the different exons of a gene provide information on which alternatively spliced mRNAs were expressed in the sample, and in which proportions. Since the read length is typically smaller than the mRNA molecule of a transcript, identifying and quantifying the transcripts is however difficult: an observed read mapping to a particular exon may come from an mRNA molecule of any transcript containing this exon. Some methods consider that the set of expressed isoforms [START_REF] Jiang | Statistical inferences for isoform expression in RNA-Seq[END_REF] or a candidate superset [START_REF] Huang | A robust method for transcript quantification with RNA-Seq data[END_REF][START_REF] Xing | An expectation-maximization algorithm for probabilistic reconstructions of full-length isoforms from splice graphs[END_REF] is known in advance, in which case the only problem is to estimate their expression. However little is known in practice about the possible isoforms of genes, and restricting oneself to isoforms that have been described in the literature may lead to missing new ones.

Two main paradigms have been used so far to estimate expression at the transcript level while allowing de novo transcript discovery. On the one hand, the Cufflinks software package [START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF] proceeds in two separate steps to identify expressed isoforms and estimate their abundances. It first estimates the list of alternatively spliced transcripts by building a small set of isoforms containing all observed exons and exon junctions. In a second step, the expression of each transcript is quantified by likelihood maximization given the list of transcripts. Identification and quantification are therefore done independently. On the other hand, a second family of methods [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF]Li et al., 2011b;[START_REF] Bohnert | rQuant.web: a tool for RNA-Seq-based transcript quantitation[END_REF]Li et al., 2011a;[START_REF] Mezlini | iReckon: Simultaneous isoform discovery and abundance estimation from RNA-seq data[END_REF][START_REF] Behr | Mitie: Simultaneous rna-seq based transcript identification and quantification in multiple samples[END_REF] jointly estimates the set of transcripts and their expression using a penalized likelihood approach. These methods model the likelihood of the expression of all possible transcripts, possibly after some preselection, and the penalty encourages sparse solutions that have a few expressed transcripts.

The two-step approach of Cufflinks [START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF] is reasonably fast, but does not exploit the observed read density along the gene, which can be a valuable information to identify the set of transcripts. This is indeed a conclusion drawn experimentally using methods from the second paradigm (see [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF]Li et al., 2011b;[START_REF] Bohnert | rQuant.web: a tool for RNA-Seq-based transcript quantitation[END_REF]Li et al., 2011a;[START_REF] Mezlini | iReckon: Simultaneous isoform discovery and abundance estimation from RNA-seq data[END_REF]. To summarize, the first paradigm is fast but can be less statistically powerful than the second one in some cases, while the second paradigm suffers from the exponential number of candidate isoforms and becomes intractable for genes with many exons. The contribution of this paper is to allow ℓ 1 -penalized regression methods from the second family to run efficiently without pre-filtering the set of isoform candidates, although they solve a non-smooth optimization problem over an exponential number of variables. To do so, we show that the penalized likelihood maximization can be reformulated as a convex cost network flow problem, which can be solved efficiently [START_REF] Ahuja | Network Flows[END_REF][START_REF] Bertsekas | Network Optimization: Continuous and Discrete Models[END_REF][START_REF] Mairal | Supervised feature selection in graphs with path coding penalties and network flows[END_REF]. Note that a significantly different approach was adopted by [START_REF] Behr | Mitie: Simultaneous rna-seq based transcript identification and quantification in multiple samples[END_REF] for non-convex ℓ 0 -penalty. Even though the problem they address is NP-hard, they perform simultaneous isoform identification and quantification, without explicitly enumerating all possible transcripts, by using mixed integer programming techniques.

The paper is organized as follows. Section 2 introduces the statistical model (Section 2.1) and the penalized likelihood approach (Section 2.2) we follow. Our model is similar to the one used by [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF], but properly models reads that cover more than two exons, effectively taking advantage of longer reads. We then reformulate the model as a path selection problem over a particular graph (Section 2.3), and present in Sections 2.4-2.6 our method called FlipFlop (Fast Lasso-based Isoform Prediction as a FLOw Problem) for solving it efficiently. Section 3 empirically compares our approach with the state of the art on simulated and real sequencing data. Our experiments show that our approach has higher accuracy in isoform discovery than methods which treat discovery and abundance estimation as two separate steps, and that it runs significantly faster than methods explicitly listing the candidate isoforms. We discuss the implications of our results in Section 4.

Method

Our approach to isoform deconvolution from RNA-Seq data consists of fitting a sparse probabilistic model, like several existing methods including rQuant [START_REF] Bohnert | rQuant.web: a tool for RNA-Seq-based transcript quantitation[END_REF], NSMAP [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF], IsoLasso (Li et al., 2011b), SLIDE (Li et al., 2011a) or iReckon [START_REF] Mezlini | iReckon: Simultaneous isoform discovery and abundance estimation from RNA-seq data[END_REF]. The reads from RNA-Seq data are modeled as a linear combination of isoforms expressions that are estimated by using the maximum likelihood principle. Because the number of candidate isoforms grows exponentially with the number of exons, the above methods are either computationally expensive for genes with many exons (such as NSMAP or SLIDE), or include a preselection step to reduce the number of candidates, which may alter the method accuracy.

The main novelty of our paper is to tackle the sparse estimation problem efficiently without pre-filtering. In the methodological section, we show that the corresponding penalized maximum likelihood estimator can be computed in polynomial time with the number of exons despite the exponential number of candidate transcripts. The key is the use of a non-trivial optimization technique based on the concept of flow in a graph [START_REF] Ahuja | Network Flows[END_REF][START_REF] Mairal | Supervised feature selection in graphs with path coding penalties and network flows[END_REF].

Statistical Model

We consider an extension of the model originally introduced by [START_REF] Jiang | Statistical inferences for isoform expression in RNA-Seq[END_REF] and used in NSMAP for estimating isoform expression for a known set of expressed transcripts. Given a gene of interest, we assume that the list of its n exons is known, and that the reads of the RNA-Seq experiments have been mapped to a reference genome. For the purpose of our work, an exon can either be defined by read alignment softwares as a cluster of reads, or from a pre-defined annotation such as the one provided by the UCSC genome browser1 . In the latter case, exons with alternative 5' donor and 3' acceptor sites are considered as two separate exons. For alternative 5' donor sites, the exon is broken down as one exon ending at the first 5' donor site, and another one starting at this same point and ending at the second 5' donor site (similarly for exons with 3' acceptor sites).

We define a bin to be an ordered set of exons. Each read is assigned to a unique bin, corresponding to the exact set of exons that it overlaps. Our model can involve bins with more than two exons. It is thus more general than the one of NSMAP, where bins are simply exons and exon-exon junctions. This extension of NSMAP is particularly useful for long reads, which often cover more than two exons. We summarize the read information by the counts y 1 , . . . , y q of reads falling in q different bins.

We consider in our model all m possible candidate isoforms consisting of an ordered sequence of exons. Each candidate isoform corresponds to a unique sequence of bins. This sequence can be generated by virtually moving a read along the candidate isoform, and recording the sets of exons that it successively overlaps.

The effective length l i of a bin i is defined as the number of positions in the candidate isoform where reads can start and be assigned to the bin. A simple computation shows that for a bin involving a single exon of length l e , we have l i = l e -L + 1, where L is the read length, while for bins involving several exons, l i = min(l left , Ll int -1) + min(l right , Ll int -1) -L + l int + 1, where l left and l right are the lengths of the leftmost and rightmost exons of the bin, and l int is the total length of the internal exons of the bin. Interestingly, we note that the effective length of a bin does not depend on the candidate isoform it is associated with. A figure illustrating the computation of the effective length is given in section 1 of the supplementary material.

We model read counts as independent Poisson random variables whose means are proportional to the bins' effective lengths and to the total abundances of isoforms associated to each bin. More formally, let us denote by U the m × q binary matrix defined as U ji = 1 if bin i is associated to isoform j and 0 otherwise, and by θ j ∈ R + the expression of isoform j (the expected number of reads per base in isoform j). Thus, m j=1 U ji θ j represents the sum of expressions of all isoforms involving bin i. We expect the observed count for bin i to be distributed around this value times the effective length of the bin l i , and therefore model the read count y i as a Poisson random variable with parameter

δ i = l i m j=1 U ji θ j . For a vector θ = [θ j ] m j=1 in R m + this yields the log-likelihood L(θ) = q i=1 [-δ i + y i log δ i -log(y i !)] , (1)
where the scalars δ i depend linearly on θ.

Maximizing the likelihood (1) allows to quantify the relative abundance of each transcript when the model only includes the list of "true" isoforms present in the sample [START_REF] Jiang | Statistical inferences for isoform expression in RNA-Seq[END_REF]. Since this list is unknown a priori, we present in the next section the sparse estimation approach that can jointly quantify and identify the transcripts using all candidate isoforms, following [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF].

Isoform Detection by Sparse Estimation

Since we do not assume that the list of expressed isoforms-i.e., such that θ j = 0-is known in advance, we endow θ with an exponential prior θ j iid ∼ E(λ) and maximize over all candidate isoforms the resulting posterior likelihood, leading to the estimator θλ = arg min

θ∈R m + [-L(θ) + λ θ 1 ] , (2)
where λ is a regularization parameter, and the ℓ 1 -norm is defined as θ 1 = m j=1 |θ j |. It is wellknown that the ℓ 1 -norm penalty and the non-negativity constraint have a sparsity-inducing effect -that is, lead to estimators θλ that contain many zeroes [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]. The parameter λ controls the number of non-zero elements in the solution θλ , i.e., of selected isoforms, with larger λ corresponding to fewer isoforms. [START_REF] Mezlini | iReckon: Simultaneous isoform discovery and abundance estimation from RNA-seq data[END_REF] claim that the ℓ 1 -penalty is inappropriate for isoform selection, a claim we disagree with. As they note, the sum of true abundances in RPKM weighted by isoform lengths is by definition the true proportion of reads coming from the gene times 10 9 . They conclude that penalizing by j θ j has little effect on the estimate. However, the sum of the estimator θλ weighted by isoform lengths has no reason to be equal to the observed number of reads mapping to the gene. There are several causes for that: model inadequacy, various noise sources, finite sample size, and bias of the estimator. Penalizing this sum therefore modifies the sparsity level of θλ as observed in our and other's experiments [START_REF] Bohnert | rQuant.web: a tool for RNA-Seq-based transcript quantitation[END_REF][START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF]Li et al., 2011a,b).

Note also that ( 2) is better adapted to long reads than the original formulation of NSMAP [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF] thanks to the use of general bins. rQuant [START_REF] Bohnert | rQuant.web: a tool for RNA-Seq-based transcript quantitation[END_REF], IsoLasso (Li et al., 2011b), and SLIDE (Li et al., 2011a) solve a similar problem where the likelihood is a simpler quadratic function, corresponding to a Gaussian model for the read counts. A difficulty with these approaches is that the dimension m grows exponentially in n making (2) intractable when n is large. For example, Li et al. (2011a) restrict themselves to experiments involving genes with less than 10 exons, due to the high computational cost for larger genes. [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF] restrict themselves to genes with less than 80 exons, but only consider candidates with transcription start/polyadenylation sites (TSS/PAS) pairs already observed in annotations, and which involve more than half of the exons of the gene. Other approaches such as IsoLasso include a filtering step to reduce the number of isoforms, similarly as Cufflinks does -in the case of single-end reads, their set of candidates is the set of isoforms returned by Cufflinks. As pointed out in Section 1, this filtering may lead to a loss of power in isoform detection, because it disregards the read density information when constructing the set of candidates. In the next section, we show that, surprisingly, problem (2) can be solved efficiently without pre-filtering the isoforms by using network flow algorithms.

Isoform Detection as a Path Selection Problem

In this section, we reformulate the isoform detection problem as a path selection problem over a particular graph. Remember that a graph G = (V, E) is composed of a finite set of vertices V and edges

E ⊆ V × V . A path is a sequence of vertices v 1 , . . . , v k ∈ V such that (v i , v i+1 ) is an arc in E for all indices 1 ≤ i < k. A graph is a DAG if it contains no path (v 1 , . . . , v k ) with v 1 = v k .
In other words, the graph does not contain any cycle.

We construct an oriented graph G = (V, E) whose vertices are the bins with positive effective length defined in Section 2.1 -each corresponding to an ordered set of exons. An edge connects bin v i to v j if they can be associated to two reads starting at successive positions in a candidate isoform. In general, v j is obtained from v i by removing the first exon of its ordered set or by adding one extra exon at the end of the ordered set, depending on the lengths of the exons composing the bin (see Figure 1). We call starting bins (respectively stopping bins) the bins that can contain a read at the left-most (respectively right-most) position of an isoform. The resulting graph is a DAG generalizing the splicing graph [START_REF] Heber | Splicing graphs and EST assembly problem[END_REF], whose vertices are single exons and edges are exon-exon junctions.

We also consider two new vertices s and t respectively dubbed source and sink, and connected to all starting (resp. stopping) bins. We do not impose any restriction on the set of TSS/PAS sites and each exon can potentially start or end an isoform. Consequently, the source is connected to all bins modeling an exon start, and the sink to all bins modeling an exon end. This leads to the definition of an extended graph G ′ = (V ′ , E ′ ) with V ′ = V ∪ {s, t} and E ′ is obtained by adding to E all edges of the form (s, v) where v ∈ V is a starting bin, and (v, t) where v ∈ V is a stopping bin. This graph construction is illustrated in Figure 1. [START_REF] Montgomery | Transcriptome genetics using second generation sequencing in a Caucasian population[END_REF] use a similar graph structure in the context of estimating the expression of a set of known annotated transcripts.

Let us denote by P the set of paths in G ′ starting from s and ending at t, which are called (s, t)paths. By construction, one easily checks that P is in bijection with the set of candidate isoforms, the path in G ′ corresponding to a candidate isoform being the set of bins that reads generated from the isoform can produce. Based on this one-to-one mapping, we can reformulate the penalized maximum likelihood problem (1)-(2) as follows: we want to find nonnegative weights θ p for each path p ∈ P which minimize:

v∈V [δ v -y v log δ v ] + λ p∈P θ p with δ v = l v p∈P:p∋v θ p , (3)
where the sum p∈P θ p is equal to the ℓ 1 -norm θ 1 since the entries of θ are non-negative. Note In that case, each bin either corresponds to a unique exon, or to a junction between two exons. G ′ in (c) is the resulting graph when the length of exon 3 is smaller than the read length. Some bins involve then more than two exons, here bins (2-3-4) and (2-3-5). The source links all possible starting bins and conversely all possible stopping bins are linked to the sink. There is a one-to-one correspondence between (s, t)-paths in G ′ (paths starting at s and ending at t) and isoform candidates. For example, the path (s, 1, 1-4, 4, 4-5, 5, t) corresponds to isoform 1-4-5. For the sake of clarity, some edges connecting s and t to internal nodes are not represented, and the length of the different bins are assumed to be equal. In (a), one unit of flow is carried along the path in red, corresponding to an isoform with abundance 1. In (b), another isoform with abundance 3 is added, yielding additional read counts at every node.

that we have removed the constant term log(y v !) from the log likelihood since it does not play a role in the optimization. This reformulation is therefore a path selection (finding which θ p are nonzero) and quantification problem over G ′ . The next section shows how (3) can further be written as a flow problem, i.e., technically a constrained optimization problem over the edges of the graph rather than the set of paths in P. A computationally feasible approach can then be devised to solve (3) efficiently, following [START_REF] Mairal | Supervised feature selection in graphs with path coding penalties and network flows[END_REF].

Optimization with Network Flows

A flow f on G ′ is defined as a non-negative function on arcs [f uv ] (u,v)∈E ′ that satisfies conservation constraints: the sum of incoming flow at a vertex is equal to the sum of outgoing flow except for the source s and the sink t. Such conservation property leads to a physical interpretation about flows as quantities circulating in the network, for instance, water in a pipe network or electrons in a circuit board. The source node s injects into the network some units of flow, which move along the arcs before reaching the sink t.

For example, given a path p ∈ P and a non-negative number θ p , we can make a flow by setting f uv = θ p when u and v are two consecutive vertices along the path p, and f uv = 0 otherwise. This construction corresponds to sending θ p units of flows from s to t along the path p. Such simple flows are called (s, t)-path flows. More interestingly, if we have a set of non-negative weights θ ∈ R |P| + associated to all paths in P, then we can form a more complex flow by superimposing all (s, t)-path flows according to

f uv = p∈P:p∋(u,v) θ p , (4)
where (u, v) ∈ p means that u and v are consecutive nodes on p. While (4) shows how to make a complex flow from simple ones, a converse exists, known as the flow decomposition theorem (see, e.g., [START_REF] Ahuja | Network Flows[END_REF]. It says that for any DAG, every flow vector can always be decomposed into a sum of (s, t)-path flows. In other words, given a flow [f uv ] (u,v)∈E ′ , there exists a vector θ in R |P| + such that (4) holds. Moreover, there exists linear-time algorithms to perform this decomposition [START_REF] Ahuja | Network Flows[END_REF]. As illustrated in Figure 2, this leads to a flow interpretation for isoforms.

We now have all the tools in hand to turn (3) into a flow problem by following [START_REF] Mairal | Supervised feature selection in graphs with path coding penalties and network flows[END_REF]. Given a flow f = [f uv ] (u,v)∈E ′ , let us define the amount of flow incoming to a node v in V ′ as f v u∈V ′ :(u,v)∈E ′ f uv . Given a vector θ ∈ R |P| + associated to f by the flow decomposition theorem, i.e., such that (4) holds, we remark that f v = p∈P:p∋v θ p and that f t = p∈P θ p . Therefore, problem (3) can be equivalently rewritten as:

min f ∈F v∈V [δ v -y v log δ v ] + λf t with δ v = l v f v . (5)
where F denotes the set of possible flows. Once a solution f ⋆ of ( 5) is found, a solution θ ⋆ of (3) can be recovered by decomposing f ⋆ into (s, t)-path flows, as discussed in the next section. The use of network flows has two consequences. First, (5) involves a polynomial number of variables, as many as arcs in the graph, whereas this number was exponential in (3). Second, problem (5) falls into the class of convex cost flow problems [START_REF] Ahuja | Network Flows[END_REF], for which efficient algorithms exist. 2 In our experiments, we implemented a variant of the scaling push-relabel algorithm [START_REF] Goldberg | An efficient implementation of a scaling minimum-cost flow algorithm[END_REF], which also appears under the name of ε-relaxation method [START_REF] Bertsekas | Network Optimization: Continuous and Discrete Models[END_REF]. Note that the approach can be generalized to any concave likelihood function, including the Gaussian model used by IsoLasso and SLIDE.

Network flows have been used in several occasions in bioinformatics. [START_REF] Medvedev | Maximum likelihood genome assembly[END_REF] solve a convex cost flow problem on a bidirected de Bruijn graph for maximum likelihood whole genome shotgun assembly. [START_REF] Montgomery | Transcriptome genetics using second generation sequencing in a Caucasian population[END_REF] introduced the formalism of flows for RNA-Seq data; however they did not perform isoform discovery but quantification from a set of known transcripts. Their formulation is a linear program, the dimension of which is the number of candidate transcripts considered, which is not a network flow problem. [START_REF] Singh | FDM: a graph-based statistical method to detect differential transcription using RNA-seq data[END_REF] uses the terminology of flows for RNA-Seq data in the context of testing differential transcription without reconstructing transcripts. Finally, while this paper was under revision, [START_REF] Tomescu | A novel min-cost flow method for estimating transcript expression with rna-seq[END_REF] published a similar method which also uses minimum cost flow techniques for isoform recovery. However, their method only involves bins corresponding to exons and exon-exon junction, and, more importantly, does not solve the penalized likelihood approach. They have therefore no principled way to balance the sparsity of the solution with its likelihood, and even mention that this leads to a NP-hard problem. To our knowledge, our work is the first to show that the sparsity-inducing ℓ 1 penalty can be integrated with the likelihood term in the language of network flow, in order to estimate a flow with large likelihood that can be easily decomposed in a number of paths as small as we wish.

Flow Decomposition

We have seen that after solving (5) we need to decompose f ⋆ into (s, t)-path flows to obtain a solution θ ⋆ of (2). As illustrated in Figure 2, this corresponds to finding the two isoforms from 2(b). Whereas the decomposition might not be ambiguous when f ⋆ is a sum of few (s, t)-path flows, it is not unique in general. Our approach to flow decomposition consists of finding an (s, t)-path carrying the maximum amount of flow (equivalently finding an isoform with maximum expression), removing its contribution from the flow, and repeating until convergence. We remark that finding (s, t)-path flows according to this criterion can be done efficiently using dynamic programming, similarly as for finding a shortest path in a directed acyclic graph [START_REF] Ahuja | Network Flows[END_REF]. We insist on the fact that the flow decomposition returns one solution of the ℓ 1 -penalized estimator given by problem 3. This problem can have several solutions yielding the same objective value, and which are typically sparse in the number of transcripts (see supplementary material, Figure S2). The non-uniqueness of the solution is not an artifact of our network flow approach, but a property of the ℓ 1 -penalized estimator. Algorithms such as SLIDE, NSMAP, or others that explicitly enumerate the candidates and minimize the parameter in the candidate space also return one of several solutions. In parallel, the stability of the estimator is investigated in section 4 of the supplementary.

Model Selection

The last problem we need to solve is model selection: even if we know how to solve (2) efficiently, we need to choose a regularization parameter λ. For large values of λ, (2) yields solutions involving few expressed isoforms. As we decrease λ, more isoforms have a non-zero estimated expression θ j , leading to a better data fit but also leading to a more complex model. A classical way of balancing fit and model complexity is to use likelihood ratio tests. [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF] chose this approach, but we found the log likelihood ratio statistics to be empirically poorly calibrated due to the typically small number of samples units -exons -and the non-independence of the observed read counts. We choose a related approach, which we found better behaved, and select the model having the largest BIC criterion [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF]. An alternative approach taken by Li et al. (2011a) would be to use stability selection [START_REF] Meinshausen | Stability selection[END_REF].

Results

We now compare our proposed method FlipFlop to Cufflinks [START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF] version 2.0.0, IsoLasso (Li et al., 2011b) version 2.6.1 and NSMAP [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF] on both simulated and real data. All experiments were run on a desktop computer on a single core of an Intel Xeon CPU X5460 3.16Ghz with 16Gb of RAM. Reads are aligned to a reference genome by using TopHat [START_REF] Trapnell | TopHat: discovering splice junctions with RNA-Seq[END_REF] version 2.0.6, and the constructed alignment files are used as input to the methods we compare. IsoLasso, Cufflinks, and FlipFlop only use these aligned reads as input, and estimate their exon boundaries and TSS/PAS from read density. NSMAP additionally requires exon boundaries and known TSS/PAS as input. For paired-end experiments, we extended our initial model designed for single-end: a pair of reads is considered as a long single-end read. When the two reads of a pair span bins potentially separated by some exons, we use heuristics based on genomic distances to decide wether or not these exons are spliced. All softwares are used with default parameters, except that for paired-end experiments we provide fragment length mean and standard deviation to IsoLasso, Cufflinks and FlipFlop. Note that all results can be easily reproduced by following the tutorials available at http://cbio.ensmp.fr/flipflop/experiments.html.

Simulated Human RNA-Seq Data

Since little is known about the true set of isoforms expressed in real data, we start our experimental validation with a set of simulations. We use the RNASeqReadSimulator software3 to generate singleend and paired-end reads from the annotated human transcripts available in the UCSC genome browser (hg19). We restrict ourselves to the 1137 multi-exons genes on the positive strand for chromosome 1, corresponding to 3709 expressed transcripts.

We follow the protocol of Isolasso (Li et al., 2011b) and consider that a transcript from the annotation has been detected by a method if it predicts a transcript that (i) includes the same set of exons, and such that (ii) all internal boundary coordinates (i.e., all the exon coordinates except the beginning of the first exon and the end of the last exon) are identical. The objective for each method is to recover a large proportion of transcripts that were used for read generation -high recall -without detecting too many transcripts that were not used to generate the reads -high precision.

Figure 3 shows the precision and recall of the compared methods on single-end and paired-end simulations. Since we expect the difficulty of the deconvolution problem to increase with the number of transcripts of the gene, we stratify the result by this number: each dot represents the precision and recall of one method for genes with a particular number of transcripts in the UCSC annotation. As expected, genes with more transcripts lead to more difficult estimation problems and decreased performances for all methods. Figure 3(a) shows single-end results for different read lengths from 100bp to 300bp and a fixed number of 1 million reads per experiment. FlipFlop clearly takes advantage of longer reads: the longer the read the better the accuracy for all transcript levels. For 100bp long reads, FlipFlop and Cufflinks show similar results, while NSMAP gives slightly better precision and recall for 2 transcript level and degraded results compared to FlipFlop for more than 4 expressed transcripts. These differences might arise from the fact that NSMAP restricts its search to the TSS and PAS observed in the annotation whereas FlipFlop estimates them from reads, and the fact that the two methods use different graphs and model selection techniques. For 300bp long reads, FlipFlop outperforms all other methods as soon as there is more than one expressed transcript. For instance for the 3-4 transcripts level, FlipFlop achieves 75% of precision and 67% of recall, while Cufflinks reaches 74% and 52% and IsoLasso reaches 64% and 51%. This demonstrates that an adapted model for long reads is critical for isoform recovery. NSMAP optimizes a similar Poisson objective function as FlipFlop but only models reads at the exon or exon-exon junction levels; it looses statistical power when the read length increases. Figure 3(b) shows paired-end results for 400bp fragment length, 20bp standard deviation, 1 million read pairs, and read lengths from 100 to 175bp. Although our model is designed for single-end reads and is particularly adapted to long reads, it shows competitive or better accuracy for paired-end reads. Once again, when the read length increases, FlipFlop performance improves proportionally more than other methods. In Figure 3(c), the read length is set to 150bp and the number of simulated reads varies from 1 million to 10 million. Increasing the coverage clearly helps FlipFlop whereas it does not change much for Cufflinks and IsoLasso. Indeed Cufflinks constructs its set of transcripts and estimates their abundances in two separate steps, and the construction of the set of returned transcripts does not take read density into account: it intends to find the smallest set of isoforms covering all the observed reads. IsoLasso is based on penalized likelihood maximization like FlipFlop and NSMAP, but starts from a restricted set of isoforms -the same set returned by Cufflinks for single-end data. Consequently, this family of methods discards some information that can help identifying the set of expressed isoforms. While results in Figure 3 are obtained with default parameters for all methods, Figure S6 of the supplementary shows performances when parameters are first tuned on an independent training set; in that case results are not significantly different. Note that the number of exons of a gene significantly affects the difficulty of isoform reconstruction. We show in the supplementary material (Figure S3) similar results to the ones presented in Figure 3 with a stratification by number of exons instead of number of transcripts. We also detail in section 6 of the supplementary material more realistic simulations that include typical library preparation and sequencing biases using another simulator-the Flux-Simulator [START_REF] Griebel | Modelling and simulating generic rna-seq experiments with the flux simulator[END_REF] which aims at modeling RNA-Seq experiments in silico-and we show results that are consistent with the ones of Figure 3.

Figure 4 shows the mean CPU time taken by each method to perform the deconvolution of genes with different sizes. Genes with more exons tend to have more candidate isoforms and experiments involving such genes are expected to take more time. Therefore, we stratify the observed times by exon number of the genes: each barplot represents the mean time taken by each method for finding the transcripts of genes having a particular number of exons. As expected, FlipFlop is always faster than NSMAP, more than a hundred times faster for genes with more than 20 exons. FlipFlop speed is comparable with Cufflinks, and about 4 times slower than IsoLasso. This is because IsoLasso maximizes its objective over a very restricted set of candidates -in these simulations never more than 9 and around 2-3 on average. Overall, FlipFlop estimates the set of expressed isoforms for 1137 genes in less than 9 minutes, i.e., about 2 genes per second. Note also that the time for data pre-preprocessing (finding exon boundaries and read counts for exons and junctions) is taken into account for all methods except NSMAP.

These simulations confirm several facts. First, methods that identify and quantify transcripts as a single penalized maximum likelihood problem show good performances and take clear advantage of an increase in coverage. Second, correctly modeling long reads allows to greatly improve the accuracy of isoform reconstruction. Third, the proposed network flow strategy allows to solve the penalized likelihood approach quickly even when the set of candidate isoforms is extremely large.

Real RNA-Seq Data

Our second round of experiments involves two independent real human embryonic stem cell data sets. They both contain about 50 million 75bp reads, either paired-end or single-end, with respectively NCBI SRA accession number SRR065504 and ERR361241.

In the experiments of Section 3.1, we generated the reads based on a known set of transcripts. In the present case, the reads come from actual human tissues, and we do not have access to the true set of expressed transcripts. Following [START_REF] Xia | NSMAP: a method for spliced isoforms identification and quantification from RNA-Seq[END_REF] and Li et al. (2011a), we choose to use the UCSC annotation as ground truth in the evaluation. Admittedly, this is not perfect as some expressed transcripts may be missing from the annotation, and some annotated transcripts may not be expressed in this particular experiment. However, agreement of the prediction with the set of known transcripts could be a good sign.

Figure 5 shows precision and recall of each method for different FPKM levels. When considering all transcripts with predicted abundances higher than 1 FPKM, FlipFlop has a higher precision for both the paired-end and single-end data sets, while Cufflinks has a better recall. For transcripts with more than 5 FPKM abundance, all methods have a similar recall, with a slight advantage to Cufflinks, while FlipFlop shows a much better precision. Section 7 of the supplementary gives additional details on the real RNA-Seq data experiments: Figure S8 shows the running time of the compared methods and Figure S9 corresponds to the precision-recall curves obtained for FlipFop when varying the model selection rule.

Discussion

Simultaneously tackling identification and quantitation using penalized likelihood maximization is known to be a powerful approach to estimate the set of expressed transcripts. However, existing ℓ 1 -regularized regression techniques cannot deal with genes that contain too many exons as the set of candidate isoforms grows exponentially with the number n of exons. By leveraging network flow optimization algorithms, we discover a few expressed transcripts among the exponential number of candidates by solving a problem with a number of variables polynomial in n.

We compared our approach to existing ℓ 1 -penalized likelihood maximization methods as well as methods that define expressed isoforms as the smallest set of transcripts covering all observed reads; the latter methods perform abundance estimation in a separate step. We observed on simulation data-where the true set of expressed transcripts is known-that, unlike the second set of methods, penalized likelihood maximization methods take advantage of an increase in read coverage. Moreover, we show that correctly modeling long reads is of primary importance for isoform recovery. Our approach, which models reads covering any number of exons, outperforms other methods for 300bp long reads. We believe this is an important improvement as RNA-Seq technologies are moving forward longer reads. Our FlipFlop method has also shown to be competitive with state-of-the-art methods on real single-end and paired-end human stem cells data, especially for transcripts whose abundance was significant. In addition, the runtime of our method was comparable to the runtime of the second set of methods, and orders of magnitude faster than existing software for penalized likelihood maximization.

We believe these results have important practical implications. In addition to the obvious gain in time when estimating the expression of transcripts for a single gene and a single sample, our approach makes the task amenable in a reasonable amount of time for all genes in a large number of samples. This is a necessary step for high-throughput differential expression studies at the transcript level, a direction we are planning to explore in future work. Differential expression studies were until now restricted to gene level studies, i.e., ignoring the transcript level information, to cases where the set of expressed transcripts was known in advance or to methods which were not using the read density to estimate the set of expressed transcripts -a less efficient approach as illustrated in our experiments. Furthermore, accurately estimating the transcript level expression for all genes of all samples in a study may lead to improvements in molecular based diagnosis or prognosis tools, as well as in clustering of samples, e.g for cancer subtype discovery. The ability of our approach to deal with splicing graphs with potentially hundreds of nodes also paves the way to efficient de novo transcript identification, where we do not restrict ourselves to annotated exons within a single gene.

  Graph G ′ when all exons are bigger than the read length. Graph G ′ when the length of exon 3 is smaller than the read length.

Figure 1 :

 1 Figure1: Illustration of the graph construction for a gene with 5 exons. The original splicing graph is represented in (a). The 5 exons are represented as vertices and an arrow between two vertices indicates a junction. The nodes of graph G ′ in (b) and (c) are bins with positive effective length denoted by gray square, as well as source s and sink t represented as circles. G ′ in (b) is the resulting graph when all exons are bigger than the read length. In that case, each bin either corresponds to a unique exon, or to a junction between two exons. G ′ in (c) is the resulting graph when the length of exon 3 is smaller than the read length. Some bins involve then more than two exons, here bins (2-3-4) and (2-3-5). The source links all possible starting bins and conversely all possible stopping bins are linked to the sink. There is a one-to-one correspondence between (s, t)-paths in G ′ (paths starting at s and ending at t) and isoform candidates. For example, the path (s, 1, 1-4, 4, 4-5, 5, t) corresponds to isoform 1-4-5.

  Reads at every node after adding another isoform.

Figure 2 :

 2 Figure2: Flow interpretation of isoforms using the same graph as in Figure1(b). For the sake of clarity, some edges connecting s and t to internal nodes are not represented, and the length of the different bins are assumed to be equal. In (a), one unit of flow is carried along the path in red, corresponding to an isoform with abundance 1. In (b), another isoform with abundance 3 is added, yielding additional read counts at every node.

  (a) Single-end reads with different lengths (100, 200, 300bp) (b) Paired-end reads with different lengths (100, 125, 150, 175bp) (c) Single-end reads with a fixed 150bp length and an increasing amout of material (1, 5, 10 million)

Figure 3 :

 3 Figure 3: Precision and recall of compared methods on simulated reads from the UCSC annotated human transcripts.

Figure 4 :

 4 Figure 4: Average CPU times in milliseconds (logarithmic scale) for the compared methods to process a gene from human simulated 100bp single-end reads.

Figure 5 :

 5 Figure 5: Precision and recall of compared methods on human embryonic stem cells data.

Centre for Computational Biology -CBIO, Mines ParisTech, Fontainebleau, France,

Institut Curie, Paris, France,

INSERM U900, Paris, France,

LBBE, Lyon, France,

LEAR Project-Team, INRIA Grenoble -Rhône Alpes, France

http://genome.ucsc.edu/

The function (5) can be decomposed into costs Cv(fv) over vertices v. The general convex cost flow objective function is usually presented as a sum of costs Cuv(fuv) over arcs(u, v). It is however easy to show that costs over vertices can be reduced to costs over arcs by a simple network transformation (seeAhuja et al., 1993, Section 2.4). Note that all arcs have zero lower capacities and infinite upper capacities.

http://alumni.cs.ucr.edu/ ~liw/rnaseqreadsimulator.html

Acknowledgments: We thank Wei Li for sharing the source code of IsoLasso for pre-processing RNA-Seq reads, and Matahi Moarii and Edouard Pauwels for insightful discussions. Fundings: This work was supported by the European Research Council [SMAC-ERC-280032 to J-P.V., E.B.]; the European Commission [HEALTH-F5-2012-305626 to J-P.V., E.B.]; the French National Research Agency [ANR-09-BLAN-0051-04 and ANR-11-BINF-0001 to J-P.V., E.B.]; the Stand Up to Cancer grant [SU2C-AACR-DT0409 to L.J.]; the National Science Foundation [SES-0835531, CCF-0939370 to J.M.]; and the CNRS [project FlipFlop, under program PEPS-BMI].

R package from the Bioconductor web site (http://www. bioconductor.org/) and more information is available at


Armelle Bauer

Florent Paclet

Violaine Cahouet

Ali-Hamadi Dicko

Olivier Palombi

François Faure

Jocelyne Troccaz

Interactive Visualization of Muscle Activity During Limb Movements : Towards Enhanced Anatomy Learning

Keywords: Anatomy Learning, Biomechanical Simulation, Real-Time, Augmented Reality, Embodiment, Motion Capture and Reconstruction

Figure 1: Interactive display of muscle activations during knee flexion. Left: a user moving in front of the Kinect sensor.

Introduction

Learning functional anatomy and understanding anatomical kinematics is important for medical studies, sport education as well as for education more generally. However, this learning process remains tedious and difficult due to the complex-ity of the domain which includes a large amount of structured, static and dynamic notions. Understanding the musculoskeletal system and its mobility usually requires 3D visuospatial abilities [HCD * 14,HCR * 09]. An important difficulty comes from the fact that this is still taught using static ob-jects from books, drawings, photos, or the dissection of cadavers [GCB * 07]. With the emergence of new technologies such as 3D modeling and animation, Augmented Reality (AR) systems, etc, new approaches of learning through interactions have been proposed [FBE * 13, PR13, KBSC14] to ease the learning process and the construction of mental representations of anatomy. Our work is based on two assumptions :

• Space visualization is fundamental to the understanding and the learning of complex articulated 3D shapes.

• Our motor system influences our cognition (Embodiment theory)

Assuming that the use of one's own body could make the learning process of musculoskeletal biomechanics more efficient led us to sketch the "Living Book of Anatomy" (LBA), an AR system to visualize the relevant anatomical knowledge of the trainee body in action. The key idea is to capture a user action (a limb motion for instance) and to visualize its effect on the internal anatomy within an AR paradigm : for instance, the internal structures are superimposed on the image of the trainee. Developing the LBA requires a novel combination of technologies including motion capture, anatomical knowledge modeling, user-to-model mapping, AR rendering. In this paper, we tackle the challenge of real-time motion capture and muscle activation display using a novel combination of state-of-the-art approaches.

Our biomechanical model is created using an ontology of anatomical structures (MyCorporisFabrica) [PUF * 14, DGFP12] connected with a complete set of 3D geometrical model of the human body [BUD * 14]. User motion is captured using a Kinect sensor associated with standard software and complemented with additional filtering. To avoid intensive, ill-posed inverse dynamics computations, muscle activation is queried in tables based on pose and velocity. Muscles are rendered in real-time as deformed 3D shapes using vertex blending methods, or as colored lines to highlight activation. Our key contribution is the real-time visualization of the human body musculoskeletal system kinetics and muscle activity.

The remainder of the paper is organized as follows. In Section 2, we briefly survey the related work covering medical learning and real-time applications. In Section 3 we then describe in detail the proposed pipeline and the input data we are working on. We present and discuss our results in Section 4. We finally conclude and sketch future work in Section 5.

Related Work

Motion capture of the human body has been a subject of great interest in computer graphics, vision, biomechanics, films and medical applications, and force space reasons we can only scratch the surface of this domain. People generally try to use marker-less motion capture systems be-cause of the lower price and the simple setup [CSC14]. This approach has also gained popularity in the biomechanics community [CSC14, CGMA10, CMC * 06]. Our project is based on it in order to provide an application everyone could use easily and at a low cost. Some recent papers analyze motion capture using Kinect sensors and compare it with the real motion captured with more sophisticated hardware [Red13, EPnB * 13].

A lot of studies were carried on about anatomy learning [DRC * 06, PR13, KBSC14, SM08], and 3D models have become increasingly popular [HCD * 14,HCR * 09,GCB * 07]. Augmented reality (AR), which consists in superimposing data on visible objects, has been increasingly used in the medical domain [DHLD * 02, JBC08, Spe12]. Our project is inspired by the work of Nassir Navab on the Magic Mirror project [FBE * 13, BKBN12]. We extend it using biomechanical data related to user motion.

Materials and Methods

Our method combines data from three sources:

• Muscle activation data for specified movements, measured during a preliminary phase (Sec.3.1)

• Real-time user motion capture and muscle activation evaluation (Sec.3.2)

• A 3D avatar displayed in the same pose as the user, and highlighting muscle activation (Sec.3.3)

Figure 2 summarizes the pipeline of the method. In this experiment we focus on muscle activation during knee flexion/extension for anatomy learning purpose. We used My-CorporisFabrica [BUD * 14] to select the entities (muscles and bones) involved in the knee flexion and extension.

Measuring activation

To evaluate musculoskeletal movement, a motion capture system is sufficient. However, muscle and/or tendon tension is difficult to evaluate due to redundancies in the musculoskeletal system. Indeed there are more actuators (muscles) than degrees of freedom in joints. This leads to an under-determined system of equations, resulting in an infinity of possible muscular coordination with different synergies and co-contraction patterns to achieve the same task. In this work, motion capture and identification are applied to the flexion-extension of the lower limb, and our goal is to evaluate and display the corresponding muscle activations.

To solve the redundancy problem, we use experimental measurement data to feed a biomechanical model of the lower limb based on anatomical and functional knowledge, in order to compute knee muscle activations including agonist and antagonist co-contraction. The acquisition of the knee flexion muscle activation was made only on the right segment of the human body, since the muscle activation should be the same for each leg. Database construction An experimental acquisition is performed on a single healthy subject (48yo / 186cm / 77kg) to provide muscle activation patterns corresponding to flexion and extension. The motion capture is performed during squat and no-load lower limb flexion-extensions using 12 Vicon cameras sampled at 250Hz. At the same time, the electromyographic activity of the mono and bi-articular flexor and extensor knee muscles is recorded at 2000 Hz (total 10 muscles) using a Biopac MP150 device.

The subject is standing on a force plate (force plate data not used in this work) and performs right leg flexion/extension cycles either in squat or no-load conditions (see picture 3) at three (low, medium and high) velocities. The subject performs 6 flexion/extension cycles in each condition. Pairs of silver-silver chloride surface electrodes were used to record the raw activity from the 10 main muscles producing the knee torque : Vastus Medialis (vast_med), Vastus Lateralis (vast_lat), Rectus Femoris (rect_fem), Tensor Fasciae Latae (tfl), Sartorius (sartorius), Gracilis (gracilis), Biceps Femoris (bi_fem_lh), Semi Tendinosus (semiten), Gastrocnemius (med_gas and lat_gas). The electrodes are positioned parallel to the muscle fibers

Data analysis

The Y and Z lateral displacements of each marker are filtered using a fourth-order Butterworth filter with zero phase lag and a cut-off frequency of 4 Hz. Intersegmental angles are calculated at ankle, knee and hip joints according to the convention shown in Figure 4. The derivatives of angular position are calculated at each time t using a centred finite difference scheme. Finally, kinematic values are averaged to give a mean cycle for each condition (Figure 5). The data is first demeaned then filtered using a Butterworth filter (4th order, band pass 20-450Hz) with no lag essentially to remove movement artefacts. RMS is calculated using a 250ms sliding window with overlapping to get the envelope of the EMG signal. RMS values is normalized using Maximal RMS value (RMSmax) detected across all trials for each muscle i :

RMS l (t) = RMS i (t) RMSmax i
Finally the normalized RMS signal is averaged to provide a mean cycle for each condition, and is filtered using a Butterworth filter (4th order, low pass 4Hz) with no lag to obtain an activation coefficient a i for each knee muscle. The activation of deep muscles (Biceps femoris: short head (bi_fem_short), semimembranosus (semimem), vastus intermedius (vast_int)), which EMG activities cannot be measured using surface electrodes, is estimated using synergist neighbor muscles. This data can then be used to provide suit- able information for anatomical structures animation in the LBA.

Motion capture

We use a Kinect with NiTE and the OpenNI SDK [Ope] to track the user actions. The process starts with a calibration (Sec.3.2.1). We then enter the main application loop. At each time, the user tracking provides us with joint angles (Sec.3.2.2), and the motion is analyzed (Sec.3.2.3).

Calibration

In the process of Kinect calibration we find the subject parameters and we transform a reference 3D anatomical avatar to match the user in the learning environment. We use the Tpose (3D avatar and Kinect initial position) to calibrate the user. Figure 7 we use a ratio between the width and height scale factors. Using these scale factors we find the transformation matrix to deform the reference avatar to fit the user global size, as illustrated in Figure 8.

Tracking

After the calibration we compute smooth rotation interpolations between the frames, while maintaining the distances between the joints. We only take into account the joints needed to control our system : hip joint, knee joint and foot joint.

NiTE re-computes skeleton joints at each frame based on the depth image. Unfortunately, the Kinect sensor offers an imprecise depth image which leads to limits in precision for the skeleton joints. To avoid inconsistencies during the tracking we save the initial skeleton at initialization time, and we enforce constant distance constraints between the joints during the tracking. Before applying our hierarchical model, we smooth the individual bone trajectories using a 5x5 convolution matrix. The filter is coupled with Dual-quaternion Linear Blending [Ken12, KCZO07] for multiple rotation interpolations. Each joint is composed of a position (3D vector) and an orientation (quaternion). We separate the joint positions and orientations to take only into account the joint rotations in the hierarchical joint model. This results in smooth rotation interpolations between frames and avoids jiggle effects.

Motion analysis

To avoid complex inverse dynamics computations, the muscle activations are retrieved from recorded data, based on joint angles and velocities. To fit our experimental data (Sec.3.1) we have to compute the angle and velocity of the tracked user as well. The current version considers only the hip. To compute the user angle we transform the hip quaternion to Euler angles and we peek the angle related to the x axis. The angular velocity is computed using finite differences.

Figure 10: experiment and Kinect : hip angle (rad) evolution during a mean cycle (flexion/extension)

The angular velocity is negative during flexion and positive during extension. Thus based on the angular velocity we know if the user is performing flexion or extension. Then, we use the hip angle to find activation per muscle using a linear interpolation of the experimental data. In the current version, we do not consider the value of the velocity.

Motion and activation visualization

The mapping from the user to the model/avatar is performed by attaching bones and muscles to the skeleton using Linear c The Eurographics Association 2014.

Blend Skinning. The muscles can be visualized in two ways: 1D muscles (muscle line of action) [The03] and 3D muscles.

• 1D (line of action) Hill muscles : to see all muscles behaviors (even underneath muscles like the "vastus intermedius"). Using anatomy transfer [AHLG * 13], we adapt the muscles of the OpenSim [AWLD10, DAA * 07] model to our system.

• 3D muscle meshes : for more realistic muscles visualization. Our 3D anatomy avatar provides meshes of lower limb muscles.

We display muscle activation using a smooth color map from blue (no activation) to red (maximum activation). In theory, muscle activation is a percentage value between 0.0 and 1.0.

To better highlight muscle activations, we choose to use the "experimental" maximum instead of "theoretical" maximum to 1.0. For a complete visualization we display all the lower limb entities. Figure 12 shows the visualization (1D and 3D muscles activation) for 4 of our users. More results are in the video attached.

Results and discussion

Our framework is implemented in C++ and runs on a commodity laptop PC. The real-time capture and activation display run at 25 Hz, with a lag of less than 0.1 s. The visual feedback is currently performed on the screen of the PC, as shown in the accompanying video.

The motion sequences were acquired for 6 young adult men with an average height of 1.81m, an average of weight of 85.3kg, and 4 young adult women with an average height of 1.60m and an average weight of 55.7kg. All of them without any type of injury or physical incapacity in their extremities. To get uniform results we work with Kinect sequences made in the same room and with the same lighting conditions. Scaling the reference model takes about 4 seconds, then the interactive visualization begins. The user performs flexions and extensions and sees in real-time his/her muscle activations on the screen. We have studied different kinds of flexion-extension movements with various velocities and loading conditions to create a set (library) of averaged lower limb movement references to feed the LBA. This library of muscle activation patterns is used to provide the feedback about the knee muscles for functional anatomy learning. In future work, this experimental process will be greatly improved using additional measurements of muscle activation in maximal voluntary contraction situation. Using an isokinetic ergometer, these activations will be related to knee moment level in various isometric, concentric and eccentric contractions. This knowledge will allow us to more accurately normalize EMG signals acquired during flexion-extension movements. Finally, the development of a biomechanical model of the lower limb coupled with inverse dynamics, optimization procedure and EMG information could produce a better muscle tension estimation in all kind of flexion-extension movement.

Gender Age Height

The Kinect sensor provides an imprecise depth image which leads to limits in precision for the skeleton joints. This leads to an approximate motion feedback for the user. A better filtering of the Kinect skeleton [WZC12, MCA07], or a completion of depth sensor data are good response to this problem.

Our user calibration is simplistic while the final goal of our project is to provide a really realistic user-specific avatar. Computing our avatar using real-time anatomy transfer [AHLG * 13, GRP10] and anthropometric measurements [START_REF] Espitia-Contreras A | Development of a Kinect-based anthropometric measurement application[END_REF] would be a good alternative. An other alternative would be to combine the Kinect skeleton and silhouette [RVHT12, CPLFR13].

Conclusion

Our framework is ready to validate the theory of Embodiment. In this paper we presented a first innovative application for the visualization of human body kinetics in realtime, by displaying muscle activation. We believe that with accurate anatomically-based models and realistic motion, learning anatomy will be eased.

The AR visual feedback will soon be implemented. In the c The Eurographics Association 2014. short term, we want to improve the personalization of the avatar to reinforce embodiment. We also hope to visualize information on other limbs and their motions. In the longer term, we plan to automatically detect the user motion and deliver knowledge accordingly.

Figure 2 :

 2 Figure 2: Pipeline of the application.

Figure 3 :

 3 Figure 3: Photos of the experiment.

Figure 4 :

 4 Figure 4: Three-link model of the lower limb.

Figure 5 :

 5 Figure 5: Hip, knee and ankle angles (rad) evolution during a mean cycle.

Figure 6 :

 6 Figure 6: Activation coefficient of knee flexor muscles during a mean cycle.

Figure 7 :

 7 Figure 7: NiTE tracking skeleton (15 bones).

Figure 8 :

 8 Figure 8: user-specific calibration (user, reference avatar, user-specific avatar)

Figure 9 :

 9 Figure 9: on the left skeleton tracking without filter, on the right skeleton tracking with filter

Figure 11 :

 11 Figure 11: key positions of our sequence from the Kinect point of view

Figure 12 :

 12 Figure 12: users : color and depth map, 1D muscles visualization, 3D muscles visualization

system in action with visualization of muscle activations.

			User Tracking
	USER	User Calibration (Kinect)	(Kinect) Skeleton : 15 joints	Filtering Temporal Consistensy	Selection of the sub-sytem (here : hip, knee, joint)
			(absolute coordinates)
					Hierarchy System
					(relative coordinates)
	AVATAR	User-Specific Avatar (Transform current Avatar)	Fixed Avatar (Save current Avatar)
					Compute transform
					(from initial to current joints)
			Compute Angle
				&
			Angular velocity
		Experimental	Muscle Activation :
		Average data	Per muscles
					Full

Muscle types : Hill Muscle (1d), Mesh Muscle (3d)

Apply Transform to Avatar

Table 1 :

 1 Kinect sequence users exact data

				Weight
	F	27	1.59	48
	F	23	1.60	50
	F	24	1.60	58
	F	31	1.63	67
	H	30	1.72	74
	H	27	1.76	110
	H	29	1.82	80
	H	26	1.83	72
	H	25	1.85	100
	H	24	1.91	76

c The Eurographics Association 2014.

Acknowledgments

Many thanks to Estelle Charleroy and Laura Paiardini our Graphist Artists for their great job on the video. To thanks Jonathan Dumon for allowing us to use the MOCA experimental platform. To all our kinect subjects. This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025).