Leander Schietgat

Celine Vens
email: celine.vens@kuleuven.be

Ricardo Cerri

Carlos N Fischer

Eduardo Costa

Jan Ramon

Claudia M A Carareto

Hendrik Blockeel

A machine learning based framework to identify and classify long terminal repeat retrotransposons

Transposable elements (TEs) are repetitive nucleotide sequences that make up a large portion of eukaryotic genomes. They can move and duplicate within a genome, increasing genome size and contributing to genetic diversity within and across species. Accurate identification and classification of TEs present in a genome is an important step towards understanding their effects on genes and their role in genome evolution. We introduce TE-LEARNER, a framework based on machine learning that automatically identifies TEs in a given genome and assigns a classification to them. We present an implementation of our framework towards LTR retrotransposons, a particular type of TEs characterized by having long terminal repeats (LTRs) at their boundaries. We evaluate the predictive performance of our framework on the well-annotated genomes of Drosophila melanogaster and Arabidopsis thaliana and we compare our results for three LTR retrotransposon superfamilies with the results of three widely used methods for TE identification or classification: REPEATMASKER, CENSOR and LTRDIGEST. In contrast to these methods, TE-LEARNER is the first to incorporate machine learning techniques, outperforming these methods in terms of predictive performance, while able to learn models and make predictions efficiently. Moreover, we show that our method was able to identify TEs that none of the above method could find, and we investigated TE-LEARNER's predictions which did not correspond to an official annotation. It turns out that many of these predictions are in fact strongly homologous to a known TE.

Introduction

Transposable elements (TEs) are DNA sequences that can move and duplicate within a genome, autonomously or with the assistance of other elements. The field of TE annotation includes various steps such as the identification and classification of TE families. In this article, we focus on these activities since accurate identification and classification of TEs enable researches into their biology and can shed light on the evolutionary processes that shape genomes [START_REF] Wheeler | Dfam: a database of repetitive DNA based on profile hidden Markov models[END_REF].

TEs in eukaryotes can be classified according to whether reverse transcription is needed for their transposition (Class I or retrotransposons) or not (Class II or DNA transposons). A consensus for a universal TE classification has not been reached yet [START_REF] Pie ´gu | A survey of transposable element classification systems-A call for a fundamental update to meet the challenge of their diversity and complexity[END_REF], but this lack of consensus does not affect the focus of our study. Here, we will follow the hierarchical system proposed by Wicker et al. [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF], which includes the levels of class, subclass, order, superfamily, family and subfamily. Fig 1 presents an illustration of Wicker's hierarchy considered in our study. Class I is composed of five orders: LTR retrotransposons, DIRS-like elements, Penelope-like elements (PLEs), long interspersed nuclear elements (LINEs) and short interspersed nuclear elements (SINEs). Similar in structure to retroviruses, LTR retrotransposons have long terminal repeats (LTRs), two normally homologous non-coding DNA sequences that flank the internal coding region and that range in size from a few hundred base pairs to more than 5 kb. Superfamilies within an order are distinguished by uniform and widespread large-scale features, such as the structure of protein or non-coding domains and the presence and size of the target site duplication (TSD). Families are defined by DNA sequence conservation and subfamilies on the basis of phylogenetic data. Class II is divided into two subclasses, which are distinguished by the number of DNA strands that are cut during transposition. Subclass 1 consists of TEs of the order TIR, which are characterized by terminal inverted repeats (TIRs). Subclass 2 groups the Helitron and Maverick orders.

Methods identifying TEs in a genome are homology-based, employ structural information or do not use prior information at all about the TEs to be identified [START_REF] Andrieu | Detection of transposable elements by their compositional bias[END_REF][START_REF] Caspi | Identification of transposable elements using multiple alignments of related genomes[END_REF][START_REF] Bergman | Discovering and detecting transposable elements in genome sequences[END_REF]. The latter methods, known as de novo repeat discovery methods, search for example for repeats in the genome. A widely used method for TE identification is REPEATMASKER [7]. This tool screens a query sequence searching for repeats, taking into account their similarity with sequences from a reference library, using an optimal pairwise alignment algorithm. CENSOR [START_REF] Jurka | Censor-a program for identification and elimination of repetitive elements from DNA sequences[END_REF] works similarly as REPEATMASKER but uses BLAST for the comparison. Afterwards, both REPEATMASKER and CENSOR remove overlaps and defragment detected repeats. Loureiro et al. [START_REF] Loureiro | Boosting the Detection of Transposable Elements Using Machine Learning[END_REF] show that machine learning can be used to improve the identification of TEs. They assessed a set of (non-machine learning based) identification methods and learn a classifier that combines their predictions to determine whether a sequence is a TE or not. Another classifier predicts the best method to determine the exact boundaries of a TE. In their analysis, both REPEATMASKER and CENSOR were the most accurate tools. While Loureiro et al. demonstrate the benefit of using machine learning models to improve predictions, they do not use machine learning to obtain the predictions, which we address in this article.

A few methods have been proposed to classify TEs. LTRDIGEST [START_REF] Steinbiss | Fine-grained annotation and classification of de novo predicted LTR retrotransposons[END_REF] evaluates a list of LTR retrotransposons generated by another tool called LTRharvest [START_REF] Ellinghaus | LTRharvest, a efficient and flexible software for de novo detection of LTR retrotransposons[END_REF], annotating these sequences w.r.t. the protein domains and other structural characteristics that were found in them. LTRDIGEST can then be used for de novo (unsupervised) classification, i.e., finding groups within the LTR retrotransposons without any predefined classification scheme. To evaluate whether the resulting groups represent known LTR retrotransposon superfamilies, Steinbiss et al. [START_REF] Steinbiss | Fine-grained annotation and classification of de novo predicted LTR retrotransposons[END_REF] have matched representative sequences of the groups to a reference set of known transposon sequences using a fixed set of rules. LTRSIFT [START_REF] Steinbiss | LTRsift: a graphical user interface for semi-automatic classification and postprocessing of de novo detected LTR retrotransposons[END_REF] takes the LTRDIGEST output and clusters the candidate sequences. It then tries to find patterns of shared cluster membership that might indicate multiple TE families, e.g. different Copia-like, Gypsy-like or Bel-Pao families. It is a generic tool that uses sequence clusters to find family-specific patterns, based on the LTRDIGEST detected features. These patterns are then used as evidence for family discrimination. TECLASS [START_REF] Abrusan | TEclass-a tool for automated classification of unknown eukaryotic transposable elements[END_REF] classifies TE sequences into Class I and Class II TEs. The Class I elements can further be classified into LTRs and non-LTRs, and the non-LTRs are classified into the SINE or LINE orders. This classification is obtained by a hierarchy of binary classifiers based on machine learning support vector machines, using oligomer frequencies as features. REPCLASS [START_REF] Feschotte | Exploring repetitive DNA landscapes using REPCLASS, a tool that automates the classification of transposable elements in eukaryotic genomes[END_REF] consists of three independent classification modules: a module based on homology information, a module that searches for structural characteristics such as LTRs or TIRs, and a module that searches for target site duplication. The three modules provide classifications at different levels of granularity, typically at the subclass or order level, sometimes at the superfamily level. Finally, an integration module aims to compare, rank, and combine the results of the three modules providing a single tentative classification. PASTEC [START_REF] Hoede | PASTEC: an automatic transposable element classification tool[END_REF] also uses multiple features of TEs to classify TE sequences: structural features (TE length, presence of LTRs or TIRs, presence of simple sequence repeats, etc.), sequence similarities to known TEs, and conserved functional domains found in HMM profile databases. It provides classifications on the order level, including all orders from the classification hierarchy defined by Wicker et al. [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF], whereas TECLASS and REPCLASS only consider a subset of the orders. Importantly, none of the above classification systems is able to provide classifications for LTR retrotransposons at the superfamily level. One exception is a recently introduced method called LTRCLASSIFIER [START_REF] Monat | LTRclassifier: A website for fast structural LTR retrotransposons classification in plants[END_REF], which performs both annotation (i.e., identifying structural elements) and classification (but not identification) for plant genomes, and returns predictions for the Copia and Gypsy superfamilies.

In this article we introduce TE-LEARNER, a framework for the identification of TEs of a particular order, and for the classification of these TEs on the superfamily level. TE-LEARNER consists of three steps. First, based on the characteristics of the order under consideration, it extracts from the genome a set of candidate sequences, which may include false positives. Second, it automatically annotates these candidates with features. Finally, the features are given as input to a machine learning model, which predicts whether a given candidate sequence is indeed a TE of the considered order, and if so, predicts its superfamily.

In particular, we present TE-LEARNER LTR , an implementation of this framework for LTR retrotransposons, which include the superfamilies Copia, Gypsy and Bel-Pao [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF]. As features we consider the occurrence of conserved protein domains, which help TEs perform the transposition process. The machine learning method we apply is random forests. This last step is essential, since the model of the three superfamilies contains the same protein domains [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF]; for Gypsy and Bel-Pao some domains even occur in the same order.

As LTR retrotransposons have a high abundance in the genomes of Drosophila melanogaster [START_REF] Bergman | Recurrent insertion and duplication generate networks of transposable element sequences in the Drosophila melanogaster genome[END_REF] and Arabidopsis thaliana [START_REF] Devos | Genome Size Reduction through Illegitimate Recombination Counteracts Genome Expansion in Arabidopsis[END_REF][START_REF] Lin | Sequence and analysis of chromosome 2 of the plant Arabidopsis thaliana[END_REF], and as these genomes are well annotated, they present the ideal candidates for evaluating how well our proposed method can identify and classify the LTR retrotransposons without using any prior information about these genomes. We present an extensive quantitative analysis on D. melanogaster and A. thaliana comparing the obtained results to three widely used methods (each dealing with one of the two tasks considered) and we show that TE-LEARNER LTR outperforms the state-of-the-art methods w.r.t. predictive performance and runtime.

The novelty of our proposed method w.r.t. the available methods lies mainly in three aspects. First, in contrast to the other methods, which focus on one task, here we consider the tasks of identifying and classifying TEs together. Second, we propose a general framework for these tasks. Even though the implementation we provide in this article focuses on LTR retrotransposons, our framework can be extended to other TE orders. Third, in contrast to classification methods such as LTRDIGEST, LTRSIFT, REPCLASS, PASTEC and LTRCLASSIFIER, our method is not based on a predefined set of rules. Instead, we exploit the strength of machine learning to automatically derive rules from the available data, with no need of prior knowledge. Our framework is the first step towards completely automatic identification and classification of TEs in superfamilies.

Methods

Framework

We address the following problem: given an unannotated genome, find subsequences in it corresponding to a particular order from the classification scheme [START_REF] Wicker | A unified classification system for eukaryotic transposable elements[END_REF], and predict their superfamily. We propose the following three-step framework, called TE-LEARNER:

1. The genome is split into subsequences, that become the candidate TE sequences.

2. Every candidate TE sequence is annotated with features related to the TE order considered.

3. Every candidate is represented by its features and fed into a machine learning model. The model predicts for every candidate the probability that the sequence belongs to a specific superfamily of the order considered.

We now discuss TE-LEARNER LTR , one particular implementation, for every step in detail, focusing on the LTR retrotransposon order. In Step 1 we use a simple splitting strategy to obtain subsequences of the genome. The features used in Step 2 are conserved protein domains known to occur in LTR retrotransposons, and the machine learning model used in Step 3 is a random forest. Note the modularity of the framework: every step can be implemented independently of the other steps. For instance, an alternative implementation could use an LTR pair detection tool in Step 1, annotate the candidates with oligomer frequencies in Step 2, and apply an artificial neural network in Step 3. Any machine learning classifier can be used, as long as it outputs a probability.

Step 1: Obtaining candidate LTR retrotransposons from the genome. First, the genome needs to be cut into subsequences, which together with their conserved protein domain information are used as input to the machine learning model. In our approach, we use a sliding window to generate all subsequences of a particular length (10,000 nucleotides in our tests) with an overlapping between them that avoids important regions potentially being cut (1,000 nucleotides in our tests). For this first step, one interested, for example, in potential full-length LTR retrotransposons could alternatively use existing tools [START_REF] Ellinghaus | LTRharvest, a efficient and flexible software for de novo detection of LTR retrotransposons[END_REF][START_REF] Xu | LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons[END_REF] to build the sets of candidates based on their LTR pairs.

The output of this first step is a list of potential LTR retrotransposon candidates, which will be searched for conserved protein domains (Step 2) and assigned probabilities of belonging to a particular superfamily (Step 3). Even though this first step generates a lot of false positives, the following steps will correct for these.

Step 2: Predicting conserved protein domains. The second step consists in screening every candidate TE, obtained previously, for the presence of conserved protein domains that can be found in LTR retrotransposons. To that aim, we use the RPS-Blast program [START_REF] Camacho | BLAST+ Release Notes[END_REF]. RPS-Blast (Reverse Position-Specific Blast) uses a query sequence to search a database of precalculated position specific score matrices (PSSMs), built from multiple sequence alignments. In order to accelerate the search process, we constructed a database that only contains the PSSMs related to conserved protein domains found in LTR retrotransposons. More precisely, our database consists of 26 PSSMs, corresponding to the list of conserved subdomains shown in Table 1.

The resulting domain hits are used to provide an estimation of the boundaries of the TE candidate. In particular, we considered that each candidate starts at the beginning of the first predicted domain and ends at the final position of the last predicted domain. In case a domain is found in the overlap region between two sequences, we merged these sequences before determining their boundaries. In case a candidate contains predicted domains in both DNA strands (direct and reverse), we split it in two or more new candidates. Note that our boundaries do not include the LTR regions.

To summarize, for each candidate TE sequence returned by Step 1, this second step produces a list of conserved protein domains, which is used to delineate the candidate. Finally, the candidates with at least one domain hit are passed on to Step 3.

Step 3: Predicting superfamily membership. In the final step of our method, we take the candidate TE sequences, represented by their predicted protein domains, and apply a random forest model on them. We use a first-order logic based format to represent candidate sequences. This format is more expressive than traditional tabular representations and can be handled by so-called relational learning systems (see [START_REF] Raedt | Logical and Relational Learning[END_REF] for more background on relational In our approach, the learning process involves learning, for each LTR retrotransposon superfamily, a separate model that maps a sequence, represented as in Fig 3 , to the probability that it belongs to that superfamily. The model is constructed using the FORF (first-order random forests) approach [START_REF] Van Assche | First Order Random Forests: Learning Relational Classifiers with Complex Aggregates[END_REF], which is implemented in the relational data mining system ACE (http://dtai.cs.kuleuven.be/ACE/). A random forest [START_REF] Breiman | Random forests[END_REF] is an ensemble of decision trees that classifies examples (candidate sequences in this case) by combining the predictions of each individual tree. The decision trees in a forest are constructed by resampling the training set and the feature set in the different tree induction processes. In this paper we use relational decision trees [START_REF] Blockeel | Top-down induction of first-order logical decision trees[END_REF] to construct the random forest. The difference with the standard approach [START_REF] Breiman | Classification and Regression Trees[END_REF][START_REF] Quinlan | C4.5: Programs for Machine Learning[END_REF] is that it allows to use the data representation defined above (instead of a traditional tabular format), to express background knowledge, and to use user-defined tests in the nodes of the trees. Apart from that, a standard top-down decision tree induction algorithm is applied.

The tests that we allow include the following: (1) the occurrence of a particular protein domain, (2) the occurrence of a particular protein domain with a certain minimum or maximum length limit (the same domain can be predicted with different lengths), and (3) the spatial relationship between two domains. The protein domains that are considered in the tree nodes are the subdomains listed in Table 1, including the more general domains of the left column. Moreover, each test includes a minimum e-value for the domain hits, which is chosen from the following list: [1e -50, 1e -40, 1e -30, 1e -20, 1e -10, 1e -05, 1e -02, 1e -01]. The list of length limits employed is [START_REF] Xu | LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons[END_REF]50,100,200,500,700]. Note that limiting the tests to individual domains or pairs of domains (rather than focusing on the complete TE domain structure) enables the classification of incomplete TEs as well. A final type of tests that are allowed in the tree nodes is (4) whether the number of occurrences of general domains (left column of Table 1) exceeds a number between 1 and 5. As such, the total number of tests considered at each tree node is 10,947.

As an illustration, Fig 4 shows the (partial) relational decision tree obtained when constructing a single tree instead of a random forest for the Gypsy superfamily. Each internal (oval) node of the tree contains one of the tests defined above, and is used to route down the sequence according to the outcome of the test. The root node, in Fig 4, tests whether the sequence has a domain RNase_HI_RT_Ty3 with an e-value of less than 0.01. If it does, the sequence moves down to the "yes" branch of the tree; otherwise it moves to the "no" branch. This procedure is repeated until the sequence arrives in a leaf (rectangle) node, which provides the probability for belonging to Gypsy. The tree shows that not all Gypsy instances are discovered by merely checking the occurrence of the RNase_HI_RT_Ty3 key domain (consider the leaf node with 42 training sequences: 97.6% of them are Gypsy elements, although they followed the "no" branch at the root).

As we can see in Table 1, domains have subdomains; thus, we provide the hierarchical "is a subdomain" relationship as background knowledge to the system. This allows the tests to check for the occurrence of both domains or subdomains, even though the data representation only contains subdomains. For example, in the tree in Fig 4 there is a test that checks for the occurrence of the domain RT and tests that check for the occurrence of its subtypes (RT_LTR, for example).

TE-LEARNER is available for download at http://dtai.cs.kuleuven.be/software/te-learner.

Methodology and parameter settings

We evaluate the predictive performance of our framework on the genomes of D. melanogaster and A. thaliana. We use version 6-15 of the annotated genome from Flybase (http://flybase. org/), as the official annotation for D. melanogaster, which was made publicly available in April 2017. We use the Flybase annotations "Transposable Elements" and "Repeat Regions" to constitute the golden standard in our experiments. For A. thaliana, we used the Araport11 annotation, released in June 2016, for genome TAIR10, from The Arabidopsis Information Resource (TAIR) (http://www.arabidopsis.org). We will compare our results for the Copia, Gypsy and Bel-Pao superfamilies (Bel-Pao only for D. melanogaster because there is no annotation for it in A. thaliana) with those of three methods for TE identification or classification that can make predictions at the superfamily level: REPEATMASKER, CENSOR and LTRDIGEST. For each superfamily, we also compare the results to those of a baseline model. We now discuss the specific parameter settings for each of the tools used in our framework, as well as for the methods we compare to.

Baseline: The baseline model starts from the TE candidates obtained in Step 2 of our framework and makes predictions solely based on the presence of one key protein domain (as predicted by the RPS-Blast program): RNase_HI_RT_Ty1, RNase_HI_RT_Ty3, and RT_pep_A17 for Copia, Gypsy, and Bel-Pao, respectively. As such, it evaluates the impact of the machine learning aspect (step 3) in our framework.

RPS-Blast:

We constructed the database used by RPS-Blast by taking for each domain of interest the set of sequences from the Conserved Domain Database (CDD) (http://www.ncbi. nlm.nih.gov/Structure/cdd/cdd.shtml) [START_REF] Marchler-Bauer | CDD: a Conserved Domain Database for protein classification[END_REF], used to generate the original multiple sequence alignment, except that we excluded the sequences of the organisms used in our tests (D. melanogaster and A. thaliana, respectively). The reason to exclude the D. melanogaster or A. thaliana sequences is that we want to provide an evaluation as blindly as possible, without using any known information from the target organism. The PSI-Blast (Position-Specific Iterative Blast) program was used to obtain the new PSSMs and Makeprofiledb application for creating the database for RPS-Blast.

FORF: The relational trees were built with default parameters, except for the minimum number of examples in a leaf, which was set to 5. No pruning was used. The forests consist of 100 trees, with a feature sample size at each node equal to the square root of the number of possible features.

For the training of FORF we used sequences from Repbase (http://www.girinst.org/server/ RepBase/), volume 17-issue 3, one set for each superfamily of interest here. In order to provide a fair evaluation, as before, we excluded from these sets the sequences of the target organism. Note that each analysis was performed twice: each time leaving out one target organism. The resulting sets were also used as the databases for RepeatMasker and Censor applicationsas described further. We ran the RPS-Blast program, with the PSSM database created in Step 2 of our framework, to search these sequences for regions related to the conserved domains of interest (Table 1), retrieving the same types of information obtained from the screening of the candidates (Step 2). We observed that the longest predicted domain region in these sequences has a length smaller than 800 nucleotides, which indicates that the overlap size of 1,000 nucleotides we used in Step 1 of our framework, is sufficient. We removed training sequences without domain hits and those that contained domain hits in both strands of their genomic sequence. The resulting number of sequences is 3188 for Copia, 4718 for Gypsy, and 891 for Bel-Pao when leaving out D. melanogaster sequences. Leaving out A. thaliana, the numbers become 3077 for Copia and 4728 for Gypsy. These sequences constitute the positive training set. For each superfamily, we also constructed a negative set, by sampling without replacement from the other superfamilies. For Copia and Bel-Pao the negative set has an equal size as the positive set; however for Gypsy, given the size of its positive set, the negative set contains less sequences (all Copia and Bel-Pao sequences), which still yields a balanced classification task.

RepeatMasker and Censor: These systems were run using their standard parameter settings. For a fair evaluation, we used as reference library the same training sets as for FORF as described above. As each of the training sets belongs to a particular superfamily, we can label hits with the corresponding superfamily. Both applications were run on the complete genomes.

LtrDigest: This method was also run with its standard parameter settings on the complete genomes. We only retained predictions with an assigned DNA strand and used the authors' guidelines to assign a particular superfamily to each prediction as follows. Every predicted sequence is annotated with protein domain hits. If the sequence has a "Peptidase_A17" hit it is classified as BelPao; otherwise, if the sequence has a "Gypsy" hit, it is classified as Gypsy; otherwise, following [START_REF] Steinbiss | Fine-grained annotation and classification of de novo predicted LTR retrotransposons[END_REF], if the sequence has an "INT" followed by an "RT" (there may be other hits in between), it is classified as Copia and if an "RT" is followed by an "INT", it is classified as Gypsy. The remaining sequences are not classified.

Evaluation methodology

We report the predictive performance of the different methods with precision-recall (PR) curves [START_REF] Davis | The relationship between Precision-Recall and ROC curves[END_REF]. The motivation for preferring PR curves over the more popular ROC curves is as follows. Only a small fraction of the genome contains TE sequences of a specific superfamily, thus we are more interested in recognizing the positive cases, i.e. the candidate sequences that actually belong to the superfamily, than in correctly predicting the negatives. Precision is the percentage of predictions that are correct and recall is the percentage of annotations that were predicted. A PR curve plots the precision of a model as a function of its recall. Assume the model predicts the probability that a new example is positive, and that we threshold this probability with a threshold t to obtain the predicted class (positive or negative). A given threshold corresponds to a single point in PR space, and by varying the threshold we obtain a PR curve: while decreasing t from 1.0 to 0.0, an increasing number of examples is predicted positive, causing the recall to increase whereas precision may increase or decrease (with normally a tendency to decrease). A domain expert can choose the threshold corresponding to the point on the curve that looks most interesting.

To consider a prediction as a true positive, we do not require it to match the exact same boundaries of the corresponding annotation of the genome, as this would be an overly strict criterion. Instead, we allow some tolerance by defining a true positive as a prediction which has a minimum overlap of 100 nucleotides with an annotation, or a prediction which overlaps a complete annotation and vice versa. Our motivation for this evaluation is that a domain expert can inspect each prediction and determine the exact boundaries of the complete TE.

Combined predictions

The random forests in TE-LEARNER LTR only make predictions w.r.t. the superfamily for which they were built. For example, one forest outputs the probability whether a sequence belongs to Copia or not. However, one might be interested in having a model that can make predictions w.r.t. many superfamilies at the same time. An advantage of such a model is that the user does not need to combine the results of individual models, avoiding conflicting predictions.

As our models output probabilities, one straightforward idea to obtain this more general model consists of selecting the superfamily with the highest probability. To avoid that a superfamily with a very low probability is predicted, we include the category None (i.e., the sequence does not belong to any of the considered superfamilies), which is predicted when none of the probabilities exceeds a certain threshold.

In this setting, we construct a single average PR curve for all superfamilies together as follows. When a sequence is predicted to have a certain superfamily, we consider it correct if the sequence indeed belongs to that superfamily. The definition of precision and recall is then as before. Thus, for precision, the denominator contains all candidate predictions, minus those predicted as None; for recall, the denominator contains all annotations (for all considered superfamilies).

We compare our results to those of LTRDIGEST, which is also able to make predictions w.r.t. different superfamilies at the same time.

Results

In Step 1 and 2 of our framework, we generated 3372 possible candidates for D. melanogaster and 2141 candidates for A. thaliana. In the third step, we constructed a random forest for each organism (i.e., leaving out sequences of this genome from the training data) and each superfamily. Table 2 shows the average number of nodes per tree, the training set accuracy for the forest, and the induction time. The latter shows the number of seconds needed to construct the entire forest on a MacBook Pro, 2.8 GHz Intel Core i5. With respect to the tests in the nodes of the trees, for all generated forests, we observed around 65% of the nodes testing for the occurrence of a domain with a certain length condition, around 20% without length condition, and around 15% checking the occurrence of one domain before an other. Interestingly, the number of occurrences of a particular domain did not show up in the forests.

We first evaluate the ability to predict each superfamily independently and then we evaluate the combined predictions, i.e., how well our framework performs in classifying a TE as Copia, Gypsy or Bel-Pao.

Predicting each superfamily independently

Before discussing each superfamily in detail, we first show for both genomes the number of predictions that were made and the average prediction length of each method and for each superfamily (Table 3). Note that TE-LEARNER LTR presents the same numbers of candidates and average length of candidates for the three superfamilies. This happens because in our implementation Steps 1 and 2 output one common candidate set for the three superfamilies. From the table it is clear that REPEATMASKERand CENSOR make a lot of predictions, which are on average much smaller than the predictions of TE-LEARNER LTR . LTRDIGEST on the other hand, makes much less predictions, which are considerably longer.

Copia. Let us first look at the curves for D. melanogaster (Fig 5). First, we observe that the curve for TE-LEARNER LTR has a maximal recall of 0.81. This is due to the candidate set returned by Step 2 of our method: the candidate set contains 171 of the 210 Copia TE sequences described in the annotations of D. melanogaster. Second, as expected, we observe that the precision of TE-LEARNER LTR 's curve is high for high thresholds and goes down as the thresholds are lowered. Third, REPEATMASKER, LTRDIGEST and the baseline model only output (positive) predictions with 100% confidence rather than giving probabilities. Therefore, they correspond to a single point in PR space. While the point of the baseline model is clearly below the curve of TE-LEARNER LTR , the points of REPEATMASKER and LTRDIGEST have a slightly higher value in either recall or precision, respectively, compared to TE-LEARNER LTR . Fourth, the curve of CENSOR is below the one of TE-LEARNER LTR , obtaining a maximal precision of 0.43, which quickly drops down close to zero. It does, however, obtain a higher recall than the other methods, with a maximum value of 0.86 (at a precision of 0.01). The curve for A. thaliana can be found in Fig 6 . The recall values for A. thaliana are lower than the ones for D. melanogaster, with TE-LEARNER LTR obtaining a maximal recall of 0.36, while CENSOR obtains a maximal recall of 0.54, and REPEATMASKER reaches a recall of 0.53. Both LTRDIGEST and TE-LEARNER LTR are able to reach a precision of 1, at a recall of 0.09 and 0.23, respectively. As before, CENSOR obtains a very low precision, and the baseline point is below TE-LEARNER's curve.

Gypsy. Fig 7 shows the results for D. melanogaster. As for Copia, we find that the candidate sequences returned by TE-LEARNER LTR only contain a subset of the known Gypsy sequences in D. melanogaster (1707 of the 3353 known Gypsy sequences). This results in a curve with a maximal recall of 0.51. The points of CENSOR and REPEATMASKER have a lower precision than the points on TE-LEARNER LTR 's curve. However, both are able to obtain a higher recall (0.80). The point of LTRDIGEST obtains a slightly higher precision (0.83) than TE-LEARNER LTR (0.81) at a recall of 0.14. However, TE-LEARNER LTR is able to increase the recall considerably without giving up on precision. Interestingly, the baseline model, which scans the candidate sequences for the presence of the RNase_HI_RT_Ty3 protein domain, is situated exactly on the curve of TE-LEARNER LTR . This shows that this protein domain is highly predictive for Gypsy.

The curve for A. thaliana can be found in point of REPEATMASKER and the curve of CENSOR are able to reach a higher recall (0.88 and 0.85, resp.) than TE-LEARNER LTR , while the latter obtains a much precision at lower recall values. LTRDIGEST's point is slightly above the curve of TE-LEARNER LTR , obtaining a precision of 0.87 at recall 0.23, while TE-LEARNER LTR obtains a precision of 0.79 at the same recall. TE-LEARNER LTR does however obtain a precision of 0.89 at a recall of 0.15.

Combined predictions

Figs 10 and 11 report the combined PR curves of TE-LEARNER LTR and the point of LTRDIGEST. For D. melanogaster, the point of LTRDIGEST obtains a slightly higher precision (0.80) than TE-LEARNER LTR (0.77) at a recall of 0.15. For A. thaliana, LTRDIGEST and TE-LEARNER LTR both obtain a precision of 1, however, the latter obtains a higher recall. Moreover, our combined model has the advantage of allowing the user to choose an appropriate threshold.

Discussion

To summarize the above analyses, first, TE-LEARNER LTR outperformed LTRDIGEST in terms of recall, while it obtains a similar precision. Second, TE-LEARNER LTR outperforms REPEATMASKER and CENSOR in terms of precision, while the latter are able to obtain a higher recall. As we can see in Table 3, REPEATMASKER and CENSOR make on average 3.7, resp. 14.2 times as many predictions as TE-LEARNER LTR in order to obtain this recall, which comes at the cost of precision. Additionally, TE-LEARNER LTR 's lower recall can be explained by the fact that a considerable amount of transposable elements either do not have the protein domains considered in this study (see Table 1) or have domains that were not detected by RPS-Blast and are for these reasons not found by TE-LEARNER LTR . Third, the fact that TE-LEARNER LTR and CENSOR output probabilities does give the advantage for a domain expert to choose a threshold based on the performance he or she prefers (either precision or recall). precision and recall, and measures the effectiveness of identifying TEs when giving equal importance to precision and recall. For the methods that output a PR (i.e., TE-LEARNER and CENSOR), we list the maximal F1 that can be obtained along the curve. The results, shown in Table 5, show that, on average, TE-LEARNER LTR obtains the best F1-score. Only in 2 cases it was outperformed by REPEATMASKER only; in all other cases it outperformed the other methods. CENSOR and LTRDIGEST are not competitive w.r.t. this evaluation measure. The F1 scores for the combined predictions confirm these results: on D. melanogaster we obtain 0.56 for TE-LEARNER LTR versus 0.26 for LTRDIGEST, on A. thaliana the values are 0.43 and 0.06, respectively.

On a MacBook Pro 2.8 GHz Intel Core i5, it took TE-LEARNER LTR less than 10 minutes per superfamily to produce predictions from the D. melanogaster and A. thaliana genomes, including the annotation of the protein domains. LTRDIGEST runs as fast as TE-LEARNER LTR , but REPEATMASKER and CENSOR have a much longer runtime: they are respectively 232 and 7 times slower than TE-LEARNER LTR (averaged over the two target organisms and the three superfamilies).

Finally, since our method outputs probabilities and achieves high precision, it can be used to discover missing annotations: if a prediction receives a high score but is not in Flybase or Araport11, it may indicate a missing or incorrect annotation. To verify this, we have performed a BLAST search against the Nucleotide Collection (nt) of NCBI, for the top five false positive predictions for each superfamily and each target organism. Default search parameters were used. The results, shown in Table 6, confirm that for D. melanogaster all these predictions match with a high confidence to a retrotransposon hit, belonging to the respective superfamily. The only exception is the top two false positive predictions for Gypsy, which show a very low query coverage. Further inspection learned that all these annotations are indeed missing in Flybase, except the first false positive for Copia, which is in Flybase, but is annotated on the wrong strand. In general, we observed that 14 annotations for Copia in Flybase describe strands different of the strands retrieved by RPS-Blast for their conserved domains; for Gypsy and Bel-Pao, these numbers are 34 and 5, respectively. For A. thaliana, only three false positives lead to a TE hit, which belongs to the respective superfamily.

Conclusion

In this paper we have proposed a framework based on machine learning to identify and classify TEs in a genome. We evaluated our approach on three Class I TE superfamilies in D. melanogaster, and two Class I TE superfamilies in A. thaliana, using a relational random forest model. We found a better predictive performance (w.r.t. F1 measure) and runtime compared to three widely used methods for TE identification and classification. In terms of F1, the performance of REPEATMASKER comes close to TE-LEARNER LTR , however, it obtains a higher recall, because it is able to recover TEs that have no conserved protein domains. The fact that we rely on these protein domains is a clear limitation of our method, yet, we are able to find TEs that other methods did not find. This suggests that TE-LEARNER LTR presents a viable alternative to the state-of-the-art methods, in case one prefers predictions with a very high precision, or as complement to the other methods when one is interested in finding more TEs. Furthermore, for our top predictions not confirmed by the official annotations, we validated their homology to known TEs of the corresponding superfamilies, showing that our method could be useful to detect missing annotations. While our implementation has been focusing on LTR retrotransposons, it is possible to train it on other TE orders with superfamilies that have recognizable protein domains. Alternatively, one could change the implementation of any of the steps of the framework: the machine learning model, the features used, and the candidate generation procedure. For instance, to identify TEs from the TIR order (a Class II order with Terminal Inverted Repeats), the first step could use software tools to identify a candidate set of sequences surrounded by TIRs (such as [START_REF] Warburton | Inverted repeat structure of the human genome: the X-chromosome contains a preponderance of large, highly homologous inverted repeats that contain testes genes[END_REF][START_REF] Gambin | TIRfinder: A Web Tool for Mining Class II Transposons Carrying Terminal Inverted Repeats[END_REF]).

A possible direction for further work is to explore hierarchical classification methods in the machine learning step of the framework. This would allow to exploit the underlying structure of the TE classification scheme. Additionally, one could try to still boost the performance of the different steps of the framework, e.g., by improving protein domain detection, or by including additional features (including features not related to protein domains) in the decision trees.

Fig 1 .

 1 Fig 1. Copia, Gypsy and Bel-Pao superfamilies positioned in Wicker's taxonomy [2]. https://doi.org/10.1371/journal.pcbi.1006097.g001

  Fig 2 shows a schematic representation of our framework based on this implementation.

Fig 2 .

 2 Fig 2. A schematic representation of our framework, applied to LTR retrotransposons. https://doi.org/10.1371/journal.pcbi.1006097.g002

Table 1 .Fig 3 .

 13 Fig 3. Logical representation. (a) Illustration of the typical structure of a TE from the Gypsy superfamily, delimited by LTRs and with protein domains identified. (b) An example of a candidate sequence, annotated with protein domain predictions. Each domain prediction lists the candidate ID, the domain, the predicted start and end positions in the sequence, and the e-value for the PSSM match. Note that domains may have subtypes. For example, RT_LTR is a subtype of domain RT. https://doi.org/10.1371/journal.pcbi.1006097.g003

Fig 4 .

 4 Fig 4. Decision tree. (Partial) decision tree predicting the probability that a given sequence belongs to the Gypsy superfamily. The abbreviation e in the nodes stands for e-value, length for the number of base pairs, seq for the number of training sequences reaching the leaf. https://doi.org/10.1371/journal.pcbi.1006097.g004

Fig 8 .

 8 Apart from the lower recall for all methods, and the baseline point being below the TE-LEARNER curve, the same conclusions as for D. melanogaster hold. Bel-Pao. Fig 9 shows the results for Bel-Pao. Here, our candidate sequences contain 349 of the 599 known Bel-Pao sequences in D. melanogaster, so the curve of TE-LEARNER LTR has a maximal recall of 0.58. Its curve is above the point of the baseline model. As for Gypsy, the

Fig 5 .Fig 7 .

 57 Fig 5. Precision-recall curves for the Copia superfamily (D. melanogaster). https://doi.org/10.1371/journal.pcbi.1006097.g005

Fig 6 .Fig 8 .Fig 9 .

 689 Fig 6. Precision-recall curves for the Copia superfamily (A. thaliana). https://doi.org/10.1371/journal.pcbi.1006097.g006

Fig 10 .

 10 Fig 10. PR curve of the combined predictions in D. melanogaster for TE-Learner LTR and LtrDigest. https://doi.org/10.1371/journal.pcbi.1006097.g010

Table 2 . Average number of nodes per tree, training set accuracy, and induction times for the random forests.

 2

	Target organism		D. melanogaster		A. thaliana
	Superfamily	Copia	Gypsy	Bel-Pao	Copia	Gypsy
	nodes	14.63	22.97	8.99	13.45	20.91
	training acc.	99.65%	99.59%	99.61%	99.71%	99.59%
	induction time	248.98s	354.23s	72.55s	186.33s	320.16s
	https://doi.org/10.1371/journal.pcbi.1006097.t002

Table 3 . Number of predictions and average prediction length for each organism, method and superfamily. D. melanogaster A. thaliana Superfamily # pred. avg. length (bp) # pred. avg. length (bp)

 3

	TE-LEARNER LTR
	Copia	3372	2387.25	2141	1779.67
	Gypsy	3372	2387.25	2141	1779.67
	Bel-Pao	3372	2387.25	-	-
	REPEATMASKER
	Copia	3662	222.25	5317	406.25
	Gypsy	27188	415.63	11987	479.63
	Bel-Pao	5176	498.78	-	-
	CENSOR
	Copia	30702	77.60	48055	99.80
	Gypsy	55645	182.35	61414	125.77
	Bel-Pao	8068	300.44	-	-
	LTRDIGEST
	Copia	43	4984.53	117	5196.29
	Gypsy	441	6860.89	25	6433.72
	Bel-Pao	152	8120.15	-	-

https://doi.org/10.1371/journal.pcbi.1006097.t003

Table 5 . F1-score of each method for each target organism and superfamily. Target organism D. melanogaster A. thaliana Avg.

 5

	Superfamily	Copia	Gypsy	Bel-Pao	Copia	Gypsy
	TE-LEARNER LTR	0.76	0.57	0.56	0.50	0.42	0.56
	REPEATMASKER	0.20	0.47	0.58	0.48	0.59	0.46
	CENSOR	0.12	0.19	0.28	0.06	0.14	0.16
	LTRDIGEST	0.34	0.24	0.36	0.16	0.02	0.22
	https://doi.org/10.1371/journal.pcbi.1006097.t005

Table 6 . BLAST analysis of the false positive (FP) predictions with the highest score.

 6 The first four columns show the location of the FP sequence, the last four columns show the details of the BLAST TE hit.

	chrom.	start pos.	end pos.	strand	hit	FP cover	e-value	identity
	D. melanogaster
	Copia
	3R	1354188	1356929	+	AF492763.1	65%	0.0	94%
	3R	2252810	2259304	-	X02599.1	70%	0.0	99%
	3R	1364299	1367899	-	X07656.1	100%	0.0	96%
	2R	3956904	3960831	-	FJ238509.1	100%	0.0	99%
	2L	23105894	23109919	-	X02599.1	100%	0.0	99%
	Gypsy
	X	21872221	21878159	-	S68526.1	1%	3e-22	96%
	2R	3790612	3798451	+	AY048125.1	0%	7e-10	89%
	X	23476916	23480018	-	AF541949.1	100%	0.0	97%
	X	23340454	23343541	-	X59545.1	100%	0.0	97%
	3R	3675376	3680906	+	X14037.1	84%	0.0	99%
	Bel-Pao
	3R	31890608	31894664	-	AY180917.1	100%	0.0	99%
	3R	3035596	3042707	-	AY180917.1	68%	0.0	99%
	3R	1657218	1661273	+	AY180917.1	100%	0.0	99%
	3L	27623316	27627575	-	AJ487856.1	100%	0.0	99%
	2R	4107747	4111802	+	AY180917.1	100%	0.0	99%
	A. thaliana
	Copia
	5	19473171	19473572	-	NM_124179.2	100%	0.0	100%
	3	15404097	15404507	-	XM_018601062.1	99%	4e-100	83%
	Gypsy
	4	5080042	5082690	+	FJ197993.1	11%	5e-37	77%

https://doi.org/10.1371/journal.pcbi.1006097.t006

PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1006097April 23, 2018

This work was supported by the Explorative Scientific Co-operation Programme between KU Leuven and São Paulo State University (UNESP), the Research Foundation Flanders (FWO-Vlaanderen) [project G.0413.09 to EC, postdoctoral grant to CV, grant GA.001.15N (Chist-ERANET call 2013 Adalab project) to JR], the Research Fund KU Leuven, ERC Starting Grant 240186 and IWT-SBO Nemoa to LS, the São Paulo Research Foundation (FAPESP -Brazil) [project 2015/14300-1 to RC, project 2012/24774-2 to CNF, project 2013/15070-4 to CMAC], the National Council for Scientific and Technological Development (CNPq-Brazil) [project 306493/2013-6 to CMAC], and Coordination for the Improvement of Higher Education Personnel (CAPES-Brazil) to EC. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data and software are available at https://dtai.cs.kuleuven.be/ software/te-learner.

In order to have an overview about how many TEs are found exclusively by each method, we created Table 4. It shows the number of annotations by every method at the highest recall point, regardless of the precision. The table confirms that REPEATMASKER and CENSOR find indeed more TEs than TE-LEARNER LTR , however, the latter is still able to find a significant number of TEs that the other methods were not able to find.

To conclude the above analysis, we computed the widely used F1 measure, which summarizes precision and recall in a single measurement. F1 corresponds to the harmonic mean of


Robrecht Cannoodt

Joeri Ruyssinck

Jan Ramon

Katleen De Preter

Yvan Saeys
email: yvan.saeys@ugent.be

IncGraph: Incremental graphlet counting for topology optimisation

Motivation

Graphlets are small network patterns that can be counted in order to characterise the structure of a network (topology). As part of a topology optimisation process, one could use graphlet counts to iteratively modify a network and keep track of the graphlet counts, in order to achieve certain topological properties. Up until now, however, graphlets were not suited as a metric for performing topology optimisation; when millions of minor changes are made to the network structure it becomes computationally intractable to recalculate all the graphlet counts for each of the edge modifications.

Results

IncGraph is a method for calculating the differences in graphlet counts with respect to the network in its previous state, which is much more efficient than calculating the graphlet occurrences from scratch at every edge modification made. In comparison to static counting approaches, our findings show IncGraph reduces the execution time by several orders of magnitude. The usefulness of this approach was demonstrated by developing a graphletbased metric to optimise gene regulatory networks. IncGraph is able to quickly quantify the topological impact of small changes to a network, which opens novel research opportunities to study changes in topologies in evolving or online networks, or develop graphlet-based criteria for topology optimisation.

Introduction

Even the simplest of living organisms already consist of complex biochemical networks which must be able to respond to a variety of stressful conditions in order to survive. An organism can be characterised using numerous interaction networks, including gene regulation, metabolic, signalling, and protein-protein interaction. The advent of high-throughput profiling methods (e.g. microarrays and RNA sequencing) have allowed to observe the molecular contents of a cell, which in turn have enabled computational network inference methods to reverse engineer the biochemical interaction networks [START_REF] Albert | Network inference, analysis, and modeling in systems biology[END_REF]. Improving the accuracy of inferred networks has been a long-standing challenge, but the development of ever more sophisticated algorithms and community-wide benchmarking studies have resulted in significant process [START_REF] Marbach | Revealing strengths and weaknesses of methods for gene network inference[END_REF][START_REF] Narendra | A comprehensive assessment of methods for de-novo reverse-engineering of genome-scale regulatory networks[END_REF][START_REF] Marbach | Wisdom of crowds for robust gene network inference[END_REF][START_REF] Bonneau | Biophysically motivated regulatory network inference: Progress and prospects[END_REF].

Several recent developments involve incorporating topological priors, either to guide the inference process [START_REF] Lopes | A feature selection technique for inference of graphs from their known topological properties: Revealing scale-free gene regulatory networks[END_REF] or post-process the network [START_REF] Ruyssinck | Netter: re-ranking gene network inference predictions using structural network properties[END_REF]. A common prior is that biological networks are highly modular [START_REF] Rives | Modular organization of cellular networks[END_REF], as they consist of multiple collections of functionally or physically linked molecules [START_REF] Hartwell | From molecular to modular cell biology[END_REF][START_REF] Baraba ´si | Network biology: understanding the cell's functional organization[END_REF]. At the lowest level, each module is made up out of biochemical interactions arranged in small topological patterns, which act as fundamental building blocks [START_REF] Milo | Network motifs: simple building blocks of complex networks[END_REF].

Graphlets [START_REF] Prz ˇulj | Modeling interactome: scale-free or geometric?[END_REF] are one of the tools which could be used to add a topological prior to a biological network, Graphlets are small connected subnetworks which can be counted to identify which low-level topological patterns are present in a network. By comparing how topologically similar a predicted network is to what would be expected of a true biological network, a predicted network can be optimised in order to obtain a better topology.

By counting the number of occurrences of each of the different graphlets (Fig 1A) touching a specific node, one can characterise the topology surrounding it. The graphlet counts of a network can be represented as a matrix with one row for each of the nodes and one column for each of the graphets (Fig 1B). An orbit represents a class of isomorphic (i.e. resulting in the same structure) positions of nodes within a graphlet (Fig 1A, coloured in red). Both graphlets and orbits have been used extensively for predicting the properties of nodes such as protein functionality [START_REF] Milenković | Uncovering biological network function via graphlet degree signatures[END_REF][START_REF] Guerrero | Characterization of the proteasome interaction network using a QTAX-based tag-team strategy and protein interaction network analysis[END_REF][START_REF] Singh | Graphlet signature-based scoring method to estimate protein-ligand binding affinity[END_REF] and gene oncogenicity [START_REF] Milenković | Systems-level cancer gene identification from protein interaction network topology applied to melanogenesis-related functional genomics data[END_REF], by performing network alignment [START_REF] Kuchaiev | Topological network alignment uncovers biological function and phylogeny[END_REF][START_REF] Milenković | Global network alignment in the context of aging[END_REF] or using them as a similarity measure in machine learning tasks [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF][START_REF] Vacic | Graphlet kernels for prediction of functional residues in protein structures[END_REF]. In this work, we focus on optimising gene regulatory networks by incorporating a topological prior as part of a topology optimisation process. We seek to optimise a predicted network by iteratively modifying the network and accepting modifications that lead to topological properties that resemble better those of true biological networks.

However, using graphlets to perform topology optimisation was hitherto not possible. Calculating the graphlet counts using the most state-of-the-art graphlet counting of a moderately sized gene regulatory network already has an execution time of about five seconds (E. coli, * 3000 genes, * 10000 interactions, up to graphlets up to 5 nodes). While this computational time poses no issue for regular static networks, recalculating all graphlet counts for every network modification made as part of a topology optimisation is computationally intractable. For example, performing 100'000 iterations of topology optimisation on a similarly sized network and calculating the topological impact of 10 possible edge modification at each iteration, already results in a computational time of about 12 days. Graphlet-based topology optimisation thus quickly becomes infeasible for larger networks.

When adding or removing an edge to a large network, the number of altered graphlets is much smaller than the total number of graphlets in the network. It could therefore be much more efficient to iterate over and count all the graphlets that have been added or removed as a result of the edge modification, than it is to recalculate the graphlet counts from scratch.

Eppstein et al. introduced data structures and algorithms for updating the counts of sizethree [START_REF] Eppstein | The h-index of a graph and its application to dynamic subgraph statistics[END_REF] and size-four [START_REF] Eppstein | Extended dynamic subgraph statistics using h-index parameterized data structures[END_REF] subgraphs in a dynamic setting. The data structures were determined such that the amortised time is O(h) and O(h 2 ), respectively, where h is the h-index of the network [START_REF] Hirsch | An index to quantify an individual's scientific research output[END_REF].

We developed IncGraph, an alternative algorithm and implementation for performing incremental counting of graphlets up to size five. We show empirically that IncGraph is several orders of magnitude faster at calculating the differences in graphlet counts in comparison to non-incremental counting approaches. In addition, we demonstrate the practical applicability by developing a graphlet-based optimisation criterion for biological networks.

Materials and methods

Assume a network G of which the graphlet counts C G are known. C G is an n-by-m matrix, with n the number of vertices in the network, m = 73 is the number of different orbits, and where C G [i, j] is the number of times node i is part of a graphlet at orbit O j . Further assume that one edge has either been added or removed from G, resulting in G 0 , at which point the counts C G 0 need to be observed. If multiple edges have been modified, the method described below can be repeated for each edge individually.

Incremental graphlet counting

As stated earlier, recalculating the graphlet counts for each modification made to the network quickly becomes computationally intractable for larger network sizes. However, as the differences in topology between G and G 0 are small, it is instead possible to calculate the differences in graphlet counts Δ G, G 0 instead. This is much more efficient to calculate, as only the neighbourhood of the modified edges needs to be explored. C G 0 can thus be calculated as

C G 0 = C G + Δ G, G 0 .
The differences in graphlet counts Δ G, G 0 are calculated by iterating recursively over the neighbours surrounding each of the modified edges (See S1 Pseudocode). Several strategies are used in order to calculate Δ G, G 0 as efficiently as possible (Fig 2). (A) The delta matrix is calculated separately for each modified edge. Since the edge already contains two out of five nodes and any modified graphlet is a connected subgraph, the neighbourhood of this edge only needs to be explored up to depth 3 in order to iterate over all modified graphlets. (B) We propose a lookup table to look up the graphlet index of each node of a given subgraph. By weighting each possible edge in a 5-node graphlet, the sum of the edges of a subgraph can be used to easily look up the corresponding graphlet index. (C) During the recursive iteration of the neighbourhood, the same combination of nodes is never visited twice. (D) As the network can be relatively large, the adjacency matrix is binary compressed in order to save memory. One integer requires 4 bytes and contains the adjacency boolean values of 32 edges, whereas otherwise 32 booleans would require 32 bytes. For example, 1GB of memory is large enough to store a compressed adjacency matrix of 92681 nodes. If the network contains too many nodes to fit a compressed adjacency matrix into the memory, a list of sets containing each node's neighbours is used instead.

IncGraph supports counting graphlets and orbits of subgraphs up to five nodes in undirected networks. By modifying the lookup table, the method can be easily extended to directed graphlets or larger-node graphlets, or limited to only a selection of graphlets. This allows for variations of the typical graphlets to be studied in an incremental setting.

Timing experiments

We compared the execution time needed to calculate the graphlet counts in iteratively modified networks between our method and a state-of-the-art non-incremental algorithm, Orca [START_REF] Hočevar | A combinatorial approach to graphlet counting[END_REF]. Orca is a heavily optimised algorithm for counting 5-node graphlets in static networks, and outperforms all other static graphlet counting algorithms by an order of magnitude [START_REF] Hočevar | A combinatorial approach to graphlet counting[END_REF].

The timing experiments were performed by generating networks from 3 different network models, making modifications to those networks while still adhering to the network model, and measuring the execution times taken for both approaches to calculate the new graphlet counts (Fig 3).

The network models were based on three static network models: Baraba ´si-Albert [START_REF] Albert | Statistical mechanics of complex networks[END_REF], Erdős-Re ´nyi [START_REF] Erdős | On random graphs[END_REF], and Geometric [START_REF] Appel | The minimum vertex degree of a graph on uniform points in [0, 1]d[END_REF]. Pseudo code for these random evolving network models can be found in S2 Pseudocode, S3 Pseudocode, and S4 Pseudocode respectively. Each model generates an initial network according to the static network model it is based on, and a list of network modifications (removing an edge from or adding an edge to the network). These network modifications are made such that at any given time point in the evolving network, it is likely that the network at its current state could have been generated by the original static network model.

Networks were generated with varying network models, between 1000 and 16000 nodes, node degrees between 2 and 64, and 10000 time points. We measured the time needed to calculate the delta matrix at random time points until 1 hour has passed. All timings experiments were carried out on Intel(R) Xeon(R) CPU E5-2665 @ 2.40GHz processors, with one thread per processor. The generation of networks with higher node counts or degrees was constrained by the execution time of the network generators, not by IncGraph. All data and scripts are made available at github.com/rcannood/incgraph-scripts.

Gene regulatory network optimisation experiments

We demonstrate the usefulness of IncGraph by using a simple graphlet-based metric to optimise gene regulatory networks. One of the striking differences between real and predicted gene regulatory networks is that the predicted networks contain highly connected subnetworks, which contain high amounts of false positives. We determine a penalty score such that predicted networks containing graphlets with many redundant edges will be penalised in comparison to very sparse networks.

The redundancy penalty (Fig 4A ) of a network is defined as the sum of occurrences of each graphlet multiplied by the redundancy associated with each individual graphlet. The redundancy of a graphlet is the number of edges that can be removed without disconnecting the nodes from one another. By using the redundancy penalty score, we aim to improve the gene regulatory network (Fig 4B).

The topology optimisation procedure uses an empty network as initialisation and grows the network by selecting interactions iteratively. Each iteration, the top k = 100 highest ranked The time for the non-incremental approach to calculate the graphlet counts of the new network is measured (T NI ).

Steps 1 to 3 are repeated until all modifications generated at step 0 are processed, or until an execution time threshold has been reached. https://doi.org/10.1371/journal.pone.0195997.g003 interactions that are not currently part of the network are evaluated, and the highest ranked interaction passing the redundancy criterion is selected (Fig 4C). This procedure is repeated until a predefined amount of time has passed. As the aim of this experiment is not to obtain the highest performing topology optimisation method, parameter optimisation of k has not been performed and is considered to be outside the scope of this work.

We optimised gene regulatory networks of E. coli and S. cerevisiae. The predicted networks were generated using the network inference method GENIE3 [START_REF] Huynh-Thu | Inferring Regulatory Networks from Expression Data Using Tree-Based Methods[END_REF] with default parameters. Gene expression data was obtained from COLOMBOS [START_REF] Moretto | COLOMBOS v3.0: Leveraging gene expression compendia for cross-species analyses[END_REF] and GEO [START_REF] Edgar | Gene Expression Omnibus: NCBI gene expression and hybridization array data repository[END_REF], respectively. The predicted networks and the optimised versions thereof were compared against respective lists of known gene regulatory interactions [START_REF] Gama-Castro | RegulonDB version 9.0: High-level integration of gene regulation, coexpression, motif clustering and beyond[END_REF][START_REF] Ma | De-novo learning of genome-scale regulatory networks in S. cerevisiae[END_REF].

Results and discussion

The contributions of this work are twofold. Firstly, we propose a new method for incrementally calculating the differences in graphlet counts in changing graphs, and show that it is orders of magnitude faster than non-incremental approaches. Secondly, we demonstrate its A) The number of redundant edges in each graphlet are counted. B) The network is optimised in order to obtain a lower redundancy over time. Two networks are shown, one before and one after the optimisation procedure. Edges coloured in red have been removed from the network after optimisation, green edges have been added. C) Starting from an empty network, the interactions are modified by iteratively evaluating the increase in redundancy of the next k interactions, and adding the first edge for which its redundancy is less than the 90 th percentile redundancy. https://doi.org/10.1371/journal.pone.0195997.g004 applicability by optimising a predicted gene regulatory network in order to reduce the false positive rate therein.

Execution time is reduced by orders of magnitude

Timing experiments show that IncGraph is significantly faster in calculating the delta matrix in comparison to calculating the graphlet counts from scratch at each iteration (Fig 5). The observed speedup ratios between IncGraph and the non-incremental approach Orca ranges from about 50× to 10000×. The speedup ratio increases as the network size increases. For larger networks, IncGraph can thus calculate the delta matrices of 10000 edge modifications while the non-incremental approach calculates one graphlet count matrix.

Surprisingly, IncGraph obtains higher execution times for networks with 5657 nodes than for networks with 8000 nodes. One possible explanation is that the size of the data structures containing those networks are particularly favourable in avoiding cache misses. Confirmation of this explanation, however, would require further investigation.

Comparing the execution time of IncGraph to the h-index of the networks indicates that the amortised time of IncGraph could be O(h 3 ) (S1 Fig). This is in line with the amortised times O(h) and O(h 2 ) of the algorithm described by Eppstein et al. [START_REF] Eppstein | Extended dynamic subgraph statistics using h-index parameterized data structures[END_REF] for counting three-size and four-size subgraphs respectively.

IncGraph allows for better regulatory network optimisation

We implemented a graphlet-based optimisation algorithm for improving the false positive rate of the predicted gene regulatory networks of E. coli and S. cerevisiae. After reranking the regulatory interactions, the F1 score of the first 1000 interactions had increased by 7.6% and 2.2% respectively (Fig 6A). The obtained speedup of about 15-30× (Fig 6B) is in line with the experiments on in silico networks. Namely, for the E. coli network at 9618 interactions and 889 nodes (average degree = 10.8), a speedup of about 30× was obtained. Similarly, for the S. cerevisiae network at 8013 interactions and 1254 nodes (average degree = 6.4), a speedup of about 15× was obtained. These speedups are in the same order of magnitude of similarly sized networks (1000 nodes and 8000 interactions) generated with a Baraba ´si-Albert model, with a speedup of 65×. This is to be expected, as such networks share the same scale-free property that gene regulatory networks have. IncGraph is significantly faster than the non-incremental approach. Note that for each interaction added to the network, the graphlet counts of 100 putative interactions were evaluated. https://doi.org/10.1371/journal.pone.0195997.g006

Conclusion

Many improvements over the past few years have resulted in efficient graphlet counting algorithms, even for large static networks. However, needing to perform the simplest of tasks tens of thousands of times quickly becomes computationally intractable. As such, recalculating the graphlet counts of a network after each of the many network modification is infeasible.

This study introduces a method for calculating the differences in graphlet (and orbit) counts, which we call incremental graphlet counting or IncGraph for short. We show that IncGraph is at least 10-100 times faster than non-incremental methods for networks of moderate size, and that the speedup ratio increases even further for larger networks. To demonstrate the applicability of IncGraph, we optimised a predicted gene regulatory network by enumerating over the ranked edges and observing the graphlet counts of several candidate edges before deciding which edge to add to the network.

IncGraph enables graphlet-based metrics to characterize online networks, e.g. in order to track certain network patterns over time, as a similarity measure in a machine learning task, or as a criterion in a topology optimisation.

Fig 1 .

 1 Fig 1. Graphlet counting in a network characterises its local topologies. (A) In total, there are 30 different graphlets containing 2 to 5 nodes, ranging from G 0 to G 29 . Orbits are an extension of graphlets which also take into account the position of a node within a graphlet. The 73 different orbits are coloured in red. (B) By counting the occurrences of these graphlets in the network, the local topology surrounding a node can be quantified. https://doi.org/10.1371/journal.pone.0195997.g001

Fig 2 .

 2 Fig 2. Several strategies are employed in order to reduce time and memory consumption. (A) Only the depth 3 neighbourhood of each modified edge needs to be explored in order to have visited all modified five-node graphlets. (B) A lookup table can be used to easily look up the graphlet index of a subgraph, by weighing each edge in a 5-node subgraph by a power of 2. (C) The same combination of five nodes is never repeated, as to avoid counting the same graphlet multiple times. (D) The adjacency matrix of the network is compressed in order to reduce memory usage. https://doi.org/10.1371/journal.pone.0195997.g002

Fig 3 .

 3 Fig 3. Static network model generators were modified to generate dynamic networks. Three network models were used: Baraba ´si-Albert, Erdős-Re ´nyi, and Geometric. Step 0: a network is generated according to the network model and the given parameters. Step 1: the network is modified such that the result is as likely to have been generated by the network model. Step 2: The time for IncGraph to calculate the differences in graphlet counts is measured (T IG ). Step 3:The time for the non-incremental approach to calculate the graphlet counts of the new network is measured (T NI ).Steps 1 to 3 are repeated until all modifications generated at step 0 are processed, or until an execution time threshold has been reached.

Fig 4 .

 4 Fig 4. Predicted gene regulatory networks of model organisms are optimised to reduce the false positive rate.A) The number of redundant edges in each graphlet are counted. B) The network is optimised in order to obtain a lower redundancy over time. Two networks are shown, one before and one after the optimisation procedure. Edges coloured in red have been removed from the network after optimisation, green edges have been added. C) Starting from an empty network, the interactions are modified by iteratively evaluating the increase in redundancy of the next k interactions, and adding the first edge for which its redundancy is less than the 90 th percentile redundancy.

Fig 5 .

 5 Fig 5.IncGraph is significantly faster than non-incremental approaches. For small networks, the execution time of IncGraph T IG is already 50 times less than that of non-incremental approaches T NI . This ratio increases even further for networks with higher numbers of nodes or higher average degrees. https://doi.org/10.1371/journal.pone.0195997.g005

PLOS ONE | https://doi.org/10.1371/journal.pone.0195997April 26, 2018

Acknowledgments

R.C. is supported by Fund for Scientific Research FWO Flanders. Ja.R. is supported by ERC StG 240186 "MiGraNT". Y.S. is an ISAC Marylou Ingram scholar.

Funding: R.C. is supported by Fund for Scientific Research FWO Flanders (www.fwo.be). Ja.R. is supported by European Research Council StG 240186 "Mi-GraNT" (erc.europa.eu).

All scripts and data in order to fully reproduce the results and figures in the paper have been made available at https:// github.com/rcannood/incgraph-scripts.

Availability

IncGraph is freely available as an open-source R package on CRAN (incgraph). The development version is also available on GitHub (rcannood/incgraph).

Author Contributions

Conceptualization: Robrecht Cannoodt, Joeri Ruyssinck, Yvan Saeys.

Investigation: Robrecht Cannoodt, Joeri Ruyssinck, Jan Ramon, Katleen De Preter, Yvan Saeys.

Methodology: Robrecht Cannoodt, Joeri Ruyssinck, Jan Ramon, Yvan Saeys.

Software: Robrecht Cannoodt.

Supervision: Joeri Ruyssinck, Katleen De Preter, Yvan Saeys.

Writing -original draft: Robrecht Cannoodt.

Writing -review & editing: Robrecht Cannoodt, Joeri Ruyssinck, Jan Ramon, Katleen De Preter, Yvan Saeys.


J F Deü

L Rouleau
email: lucie.rouleau@cnam.fr

D Matignon
email: denis.matignon@isae.fr

TIME DOMAIN FINITE ELEMENT ANALYSIS OF STRUCTURES WITH FRACTIONAL VISCOELASTIC DAMPING USING TIME-DIFFUSIVE SCHEME

This work aims at simulating the time response of a structure damped by viscoelastic materials. The structure is discretised by finite elements and a 4-parameter fractional derivative model is used to describe the frequency-dependency of the mechanical properties of the viscoelastic material. The proposed approach combines a classical Newmark time-integration scheme to solve the semi-discretised equation of motion with a diffusive representation of fractional derivatives. This approach is applied to a finite element model, and validated on a single degree-of-freedom system for which an analytical solution can be derived.

INTRODUCTION

The importance of fractional calculus for modeling viscoelastic material behavior has been recognized by the mechanical scientific community since the pioneering work of Bagley and Torvik [START_REF] Bagley | Fractional calculus -a different approach to the analysis of viscoelastically damped structures[END_REF]. The merits of using fractional differential operator lie in the fact that few parameters are needed to accurately describe the constitutive law of damping materials and the resulting model can be easily fitted to experimental data over a broad range of frequencies. While the use of such models is quite straightforward in the frequency domain, some difficulties arise from their application in the time-domain, due to the presence of fractional derivatives. The resolution methods are classically either based on time discretization of the fractional dynamics (see e.g. [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]), or on diffusive representations (cf. [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF]). For large scale systems, the first method proves memory consuming because it is necessary to store the whole displacement history of the system due to the non-local character of the fractional derivatives. The second method, based on diffusive realizations of fractional derivatives, is numerically more efficient because it has no hereditary behavior, thus avoiding the storage of the solution from all past time steps. The diffusive representation, coupled with a Newmark integration scheme, has already been developed and validated for a fractionally damped single-of-freedom system [START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF]. In [START_REF] Deü | A numerical scheme of time-domain FE analysis of viscoelastic structures with fractional derivative contitutive equations[END_REF], an extension of this approach to viscoelastic structures using FE modeling and a fractional derivative model has been presented but not tested. The purpose of this work is to implement the method described in [START_REF] Deü | A numerical scheme of time-domain FE analysis of viscoelastic structures with fractional derivative contitutive equations[END_REF] and to apply it to a structure with viscoelastic damping. The approach is validated on a single degree-of-freedom system for which an analytical solution can be derived.

FINITE ELEMENT VISCOELASTIC PROBLEM

We consider a structure composed of elastic and viscoelastic materials. A fractional derivative model is identified to described the frequency-dependency of the complex shear and the bulk moduli (resp. Ĝ and K) of the viscoelastic material:

Ĝ(ω) = G 0 + (G ∞ -G 0 )(iωτ G ) α G 1 + (iωτ G ) α G and K(ω) = K 0 + (K ∞ -K 0 )(iωτ K ) α K 1 + (iωτ K ) α K (1)
where G 0 and K 0 are relaxed moduli, G ∞ and K ∞ are unrelaxed moduli satisfying G ∞ > G 0 and K ∞ > K 0 , τ G > 0 and τ K > 0 are relaxation times and α G and α K are fractional coefficients comprised bewteen 0 and 1. Figure 1 shows that the fractional derivative model enables a good representation of the frequency-dependency both the shear and the bulk moduli over a wide frequency range with few parameters. The finite element discretization of the equation of motion leads to the following matrix system:

K e + iω ĥG (ω)K G v + iω ĥK (ω)K K v -ω 2 M Û = F (2)
where

K e = K ep + G 0 K G v 0 + K 0 K K v 0 , K G v = (G ∞ -G 0 ) K G v 0 , and K K v = (K ∞ - K 0 ) K K v 0 .
The matrix K ep is the stiffness matrix associated to the volume of elastic part of the model, M is the mass matrix of the whole system, and the stiffness matrix associated to the volume of viscoelastic material (computed with unitary moduli), is separated into a spheric part K K v 0 and a deviatoric part K G v 0 . According to Equation (1), the functions ĥG (ω) and ĥK (ω) are expressed as:

ĥG (ω) = τ α G G (iω) 1-α G [1 + (iωτ G ) α G ] and ĥK (ω) = τ α K K (iω) 1-α K [1 + (iωτ K ) α K ] (3)

COUPLED NEWMARK-DIFFUSIVE SCHEME

Equation ( 2) can be rewritten in the time domain as follows:

M Ü + h G (t) K G v + h K (t) K K v U + K e U = F(t) (4)
where the symbol represents a convolution product. Following e.g. [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF][START_REF] Deü | A numerical scheme of time-domain FE analysis of viscoelastic structures with fractional derivative contitutive equations[END_REF], letting V = U, the functions h G (t) and h K (t) will be realized by a standard diffusive representation of the form:

∂ϕ(ξ, t) ∂t = -ξϕ(ξ, t) + V(t), with ϕ(ξ, 0) = 0 (5)
observed through the continuous superposition:

(h j V)(t) = ∞ 0 µ j (ξ)ϕ(ξ, t)dξ with j = G, K (6)
This exact diffusive representation can be approximated as follows:

∞ 0 µ j (ξ)ϕ(ξ, t)dξ ≈ N n=1 µ j n ϕ(ξ n , t) with j = G, K (7)
where N is the number of approximation nodes, ξ n a sequence of angular frequencies in the frequency range of interest and µ G n and µ K n are the corresponding optimal weights computed by minimising the respective functions C G (µ G ) and C K (µ K ) defined as [START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF]:

C j (µ j ) = L l=1 N n=1 µ j n iω l + ξ n - 1 (iω l ) 1-α j 2 with j = G, K (8)
MEDYNA 2017

25-27 Apr 2017, Sevilla (Spain)

where ω l are angular frequencies and L >> N . This diffusive representation is integrated into a Newmark integration scheme, by considering the functions ϕ n (ξ n ) := ϕ(ξ n , t) as internal variables updated at each time steps. More details on the time integration scheme have been presented in [START_REF] Deü | A numerical scheme of time-domain FE analysis of viscoelastic structures with fractional derivative contitutive equations[END_REF].

ANALYTICAL SOLUTION FOR A SINGLE-DEGREE-OF-FREEDOM SYSTEM

The equation of motion for a single-degree-of-freedom system is:

k e + (iωτ ) α 1 + (iωτ ) α k v -ω 2 m û = f , (9)
and can be rewritten in the time domain as:

mτ α (D t ) 2+α + m(D t ) 2 + (k e + k v )τ α (D t ) α + k e u(t) = [1 + τ α (D t ) α ] f (t) (10)
where (D t ) β represents the time derivative of order β (integer or fractional), and α = p/q, with p and q integers satisfying p/q ∈ [0, 1].

To analytically solve this equation, the exact solution is expressed in terms of fractional power series:

u(t) = ∞ n=0 u n t n q (11)
where the coefficients u n are calculated from initial conditions and recurrence relationships.

CONCLUSION

The coupled Newmark-diffusive scheme described in this paper will be used to compute the time response of a viscoelastically damped structure, modelled by 3D finite elements. Validation of the proposed approach will be carried out on the single-degree-of-freedom system, by comparing the analytical solution with that obtained by the proposed approach.

Figure 1 :

 1 Figure 1: Master curves of the complex shear and bulk moduli of Deltane 350 (Paulstra R ) measured (points) by DMA and fitted (lines) by a fractional derivative model (left). Real part of the complex Poisson ratio from identified models (right).