A tandem repeat (TR) in genomic sequence is a subsequent recurrence of a single sequence motif. TRs are described by the length of the minimal repeating motif (unit), the number of units, and the similarity among its units. The similarity of initially identical TR units fades with time through point mutations and indels, masking their shared ancestry. Diverged TR units, even when unrecognizable by eye, can maintain structural similarity over long evolutionary times [e.g., **Figure 1** in Kajava (2012)]. While the mechanisms shaping TRs are poorly understood, they can evolve by duplication/loss of TR units, recombination, and gene conversion (Pearson et al., 2005; Richard et al., 2008). TRs can mutate by replication slippage (Levinson and Gutman, 1987; Ellegren, 2000), whereby the mispairing of a slipping-strand during the DNA synthesis causes a loss or gain of units as loops of TR units form hairpin structures (Mirkin, 2006).

Tandem repeats are frequent in coding and non-coding DNA in species throughout the kingdoms of life. Genomic TRs are a rich source for genetic variability, providing a wide range of possible genotypes at a given locus (Nithiananthrajah and Hannan, 2007) and apt opportunity for selection, not only on long evolutionary scales but also during somatic cellular processes. Particularly staggering

variation in TR unit lengths and numbers is characteristic to genomic TRs, such as ribosomal DNA arrays crucial for the translation machinery, and satellite DNA comprising the main component of functional centromeres (Richard et al., 2008). In protein-coding genes, mutations in TRs are likely to alter the structure/function of the protein product. Even in non-coding TRs, mutations can have serious fitness consequences by affecting gene regulation, transcription, and translation (Usdin, 2008). Crucially, TRs have attracted attention because of their medical relevance: many human proteins with TRs have been linked to monogenic disorders, typically affecting the nervous system (Siwach and Ganesh, 2008; Hannan, 2010). There is a high incidence of TRs in virulence factors of pathogenic agents, toxins, and allergens (Jorda et al., 2010).

### **The Methodological Challenges for Tandem Repeats Detection**

Systematic analyses of genomic TRs will help to better understand the biological processes governing and governed by TRs and their functional relevance. Such studies rely on the large-scale TRs detection (TRD). Numerous methods for TRD have been developed. Yet, since they are based on different algorithmic paradigms and heuristics, there is a large discrepancy between TR annotations produced by different algorithms for the same sequence (Leclercq et al., 2007; Merkel and Gemmell, 2008; Mudunuri et al., 2010; Schaper et al., 2012). For example, with four TRD methods applied to the human proteome, the majority of TRs were annotated by a single detector, only 9.8% were annotated by two and a meager 1.1% by at least three (Schaper et al., 2012).

The TR heterogeneity contributes to the large variability among TRD methods. The TRD is relatively simple for identical units. If the TR motif is unknown, this task is computationally expensive for long sequences, requiring an exhaustive search with no information on TR unit length, number, or position in the sequence (search space in O(*N* 3 ) for sequence length *N*). Substitutions and indels in the TR region cause major challenges to TRD: with decreasing unit similarity, TR regions become hard to discern. Indels introduce length variability between individual TR units, increasing the TR search space to O(2*NN* 3 ).

Furthermore, the original TR unit boundaries can be shifted due to new unit duplications (**Figure 1A**, Benson and Dong, 1999; Rivals, 2004). Clear boundaries are preserved only in some cases, for example, when protein TRs are confined by the exon structure. Therefore, unambiguously dividing a TR region into units of similar lengths may not accurately reflect the TR duplication history. Occasionally, the TR history is described by different phylogenies for different parts of the TR motif (**Figure 1B**). Thus, defining the consensus TR motif is problematic, and TRD methods typically differ in the predicted unit lengths and boundaries.

Ultimately, TRD methods differ by TR definitions. One TR definition borrows from string matching in computer science, whereby TRs are defined by a repetitive regular expression, allowing for a fixed proportion of dissimilar characters among TR units. This viewpoint enables straightforward exhaustive TRD algorithms, but lacks biological interpretation. Alternatively, from a structure perspective, protein TRs may be defined by structural repetitions, which allows TR detection for structurally conserved TR units, even with low sequence similarity [see, e.g., the structural TR database Repeats DB, Di Domenico et al. (2014)]. Yet, defining structural repeats is in itself problematic. Finally, the evolutionary definition states that a TR stems from ancestral unit duplications. This viewpoint has a direct biological interpretation reflecting the TR generating mechanism. Most TRD methods, however, lack an explicit TR definition, which obscures the search objective and TRs are typically detected by empirical properties of unit similarity. This further impedes the ability to evaluate the statistical properties of a method. Improving TRD requires a rigorous statistical framework based on a clear TR definition described as a biologically meaningful mathematical model. Then, a genuine TR can be distinguished from a non-TR sequence by comparison with a model describing random sequences. Relying on a mathematical TR model means that the TRD method's behavior can be predicted for different scenarios, including the evaluation of false predictions.

One possibility is to define TR units as related by a common ancestral unit under a Markov substitution model and a standard phylogeny model reflecting the unit duplication history (Schaper et al., 2012). The evolutionary distance *t* from currently observed TR units to their common ancestor can be estimated by maximum likelihood. For any tentative TR region with predicted units, this conveniently allows for statistical hypothesis testing using a likelihood ratio test (LRT; Schaper et al., 2012). The null hypothesis "*t* = *∞*" represents that the estimated evolutionary distance is so large that TR units have no common origin and could have appeared by chance. Rejecting the null suggests that the given units are related (with finite *t*) and therefore are assumed to be TRs by definition. Such approach provides a statistical framework to validate predicted TRs, filtering out potential false positive predictions.

The variability in TR annotations produced by different TRD methods warns against relying on one specific algorithm. Different methods not only achieve optimal power for different combinations of TR divergence and unit length, but also vary in their accuracy across the TR space. Therefore, to obtain the most complete and accurate set of TR annotations for a given sequence, we suggest that multiple TR detectors should be used [maximizing the number of true positives, e.g., as in Pellegrini et al. (2012)] followed by validation with an LRT (controlling false positives at a fixed significance level).

### **Annotating TRs with Sequence Profile Models**

Many common protein domains found in tandem are listed in sequence profile databases, such as Pfam (Punta et al., 2011), PROSITE (Sigrist et al., 2010), SMART (Letunic et al., 2012), Repbase (Jurka et al., 2005), and Dfam (Travis et al., 2013). For example, of all *de novo* annotated TRs with unit length *≥*15, we found that only few had not been described in Pfam (2.1% in *Arabidopsis thaliana* and 11.5% in human). Thus, TR annotation can profit strongly from the existing databases.

Profile-based annotation typically relies on sequence profile hidden Markov models (HMMs). Circular connections in an HMM allow the annotation of full TR units (Bucher et al., 1996; Schaper et al., 2014), as implemented in pftools (Sigrist et al., 2013) and in our Python TR Annotation Library TRAL (http://elkeschaper.github.io/tral/). General profile HMM annotation can be used to detect TR units [e.g., HMMER; Eddy (2011)], but a subsequent analysis is required to annotate the whole TR region, potentially including diverged TR units that, without considering the whole TR region, lack statistical significance.

Importantly, sequence profile HMMs can be used to refine *de novo* annotations. *De novo* detected TR units provide seed motifs that can be converted to sequence profile models (e.g., with *hmmbuild* from the HMMER package). These models can then be used to re-annotate sequences. The advantage is twofold: first, the quality of annotation is homogenized among TRs from different TRD methods; second, annotations on homolog sequences become comparable.

## **An Example Pipeline for Meta-Prediction of Genomic TRs**

A plausible TR annotation pipeline in three steps could be (1) identify a putative TR unit seed motif; (2) detect all tandem occurrences of this motif in the sequence, forming a putative TR; and (3) for each putative TR validate its statistical significance and filter out redundant predictions. We describe this TR annotation workflow in **Figure 2**; all functionalities are implemented in TRAL (http://elkeschaper.github.io/tral/).

Tandem repeat seed motifs are obtained from sequence profile databases and multiple *de novo* TRD algorithms. Circular profile HMMs are built from these TR seeds and consequently used to annotate TR regions in a sequence. All annotated TRs must be statistically validated, for example, using an LRT. A multiple testing correction may be required dependent on the application [e.g., Saville (1990)]. TRs that fail the test are assumed to be false positive predictions and are discarded. Combining several methods leads to redundant predictions, which is not limited to the rare case where two TR predictions fully coincide. Due to differences in predicted unit boundaries or unit numbers, it is often difficult to decide whether two overlapping TR annotations describe the same TR or, rather, nested or neighboring TRs (**Figure 1C**). To filter redundant predictions, several *ad hoc* criteria may be used. For example, overlapping TRs may be seen as redundant (the

required degree of overlap demands another *ad hoc* decision; **Figure 1C**, b–d). Using the evolutionary TR definition, we propose one possible criterion based on the representation of a TR region as a multiple sequence alignment (MSA) of its TR units. Characters grouped in one MSA column are assumed to derive from a common ancestral character. Given this, one of two TR annotations may be seen as redundant if the characters in two annotations are grouped into columns similarly by the two corresponding MSAs (**Figure 1C**, b). Model-based test statistics may be defined to recognize redundant TRs, deciding if two independent TR models provide a better description of these TRs compared to one common model.

We used *de novo* and profile-based methods to annotate TRs in the entire UniprotKB/Swiss-Prot (UniProt Consortium, 2014, v2013-08) with the proposed pipeline (**Figure 2**). A considerable proportion of these sequences contain TRs: 46% in Eukaryotes, 29% in Archaea, and 27% in Bacteria. These TRs, across all kingdoms of life, are dominated by microsatellite (typically *<* 10 bp) or short minisatellite TRs (~10–100 bp) with few units.

## **Probabilistic MSA with TRs**

When aligning homologous sequences, often we are not aware of the presence of TRs. Yet, due to uneven TR unit gain/loss among the homologs, TR-containing MSAs are error-prone, since standard indel penalty schemes do not account for the potential variation in TR region length. Some MSA methods accounting for sequence repeats produce local alignments (Raphael, 2004; Phuong et al., 2006; Treangen et al., 2009). However, a global alignment is required for evolutionary inferences, such as for the estimation of TR unit history. Sammeth and Heringa (2006) proposed a global MSA method with fixed TR unit boundaries. Yet, unit boundaries may be distorted by indels, slippage, and recombination. Modeling TRs explicitly in the MSA graph representation allows TR units to start at any position, adequately penalizing indels corresponding to unit gains/losses, and to reconstruct the evolutionary history of TR unit events. The implementation Pro-GraphMSA + TR (Szalkowski and Anisimova, 2013) uses a probabilistic phylogeny-aware approach similar to PRANK (Löytynoja and Goldman, 2005), achieving not only improved alignment quality, but also *a posteriori* estimation of rates of evolutionary events, such as TR unit indels. For example, ProGraphMSA + TR was applied to leucine-rich repeats (LRRs) in a gene family of type III effectors determining the pathogenicity in agriculturally important bacteria *Ralstonia solanacearum*. The estimates of TR indel frequencies in different clades of a gene phylogeny suggested that TR indel rate variation contributes to the diversification of this protein family [Figure 9 of Szalkowski and Anisimova (2013)]. Variation in LRR unit numbers might contribute to adaptive processes in this gene and to pathogenesis on different plant hosts.

### **Phylogenetic Approach to Study the Evolution of TRs**

Tandem repeat unit phylogenies reconstructed from homologous TRs carry much evolutionary signal, even with short units (Schaper and Anisimova, 2014; Schaper et al., 2014). These phylogenies inform about unit duplication histories and TR unit gain/loss rates, allowing to study selection on TRs and their functional relevance. Clustering patterns in TR unit phylogenies describe the unit conservation between species. If the TR unit number and order in orthologs regions from different species are perfectly conserved throughout the evolution, then the phylogeny of all TR units consists of clades formed by orthologous unit copies, each reflecting the phylogeny of the whole region (or species) [Figure 1B of Schaper et al. (2014)]. In contrast, if TR regions are fully separated, then a speciation event is followed by a series of TR unit gain/loss events, and the TR unit phylogeny consists of species-specific monophyletic clades of TR units [Figure 1A of Schaper et al. (2014)].

Tandem repeat unit phylogenies from pairs of orthologs can be used to backtrack the evolution of TRs from a single species (Schaper et al., 2014) or across an entire species tree – using the all-against-all pair-wise approach (Schaper and Anisimova, 2014). In contrast to multispecies TR unit phylogenies, for pairwise TR unit phylogenies, the statistical significance of observing perfect conservation and separation patterns is computed exactly (Schaper et al., 2014). On the other hand, multispecies TR unit phylogenies suffer fewer reconstruction errors due to the additional information on the unit evolution from additional orthologs.

Our large-scale analyses of eukaryotic proteomes revealed an extremely deep conservation of some protein domain TRs (dTRs), many dating to hundreds million years ago and some even to the times of separation between human and yeast (0.6–1.6 billion years ago) or red algae and green plants (~1.6 billion years ago). Conserved dTRs span much of the TR diversity of proteomes. For example, in human 81% of detected distinct dTR types have been conserved at least to the ancestor of mammals at least in one protein (Schaper et al., 2014). The distribution of conserved domain types is highly heterogeneous: 68% of all conserved dTRs are described by only 5% of all TR types detected in human. Yet, many conserved TRs are rare and occur only in a single protein. Similar numbers were observed in plants (Schaper and Anisimova, 2014).

In contrast, very few dTRs have separated between closely related species. In human, dTRs separated within mammals are dominated by zinc finger repeats (~50%), followed by DUF1220 (~8%; Schaper et al., 2014). In *A. thaliana*, dTRs separated within magnoliophytes are dominated by LRRs (~40%), followed by ankyrin (~12%) (Schaper and Anisimova, 2014). For both species, separated dTRs are enriched in *de novo* annotations, i.e., rare and presumably recent dTRs, which are perhaps more prone to unit gains/losses due to relaxed selection or due to higher mutation rates as a result of a high among-unit sequence similarity.

## **Simulating Genomic TRs**

Benchmarking and hypothesis testing in bioinformatics must often rely on simulated data since the truth is rarely known. For example, only the observation of identical TR units qualifies them as a true TR with certainty. Yet this is only the simplest scenario, which is of little practical relevance. With diverged or shifted TR units, unbiased benchmarks are challenging to construct from real data. Model-based simulation offers a powerful means to benchmarking and hypothesis testing in bioinformatics studies of TRs. Simulations enable comparisons of competing hypotheses and help to reveal methodological weaknesses or detect and estimate important factors. Simulations are not only crucial for benchmarking new TRD methods, but also to study the underlying evolutionary mechanisms of genomic TRs. For example, to describe the biological mechanisms for specific TR types, patterns or parameter estimates observed in these data can be compared with those obtained from sequences simulated with alternative models of TR evolution. Sequences with TRs may be generated not only with the dedicated models of evolution by unit gain/loss with fuzzy units boundaries, e.g., SlippageSim (Szalkowski and Anisimova, 2013), but also with other general sequence simulators that allow gene family evolution by mutations, indels, gene gain/loss, recombination, etc. [e.g., Dalquen et al. (2012)]. For example, sequences with TRs of different divergence and unit lengths were used to benchmark the MSA method Pro-GraphMSA + TR that accounts for sequence TRs (see above). The evaluation of the power of TRD methods also relies on simulated sequences with TRs (the alternative hypothesis). Yet the high power is irrelevant without the evaluation of false positive TRD rates, which must be done on TR-free data (the null hypothesis). Furthermore, simulated TR-free data helps to validate TR-specific findings: the comparison of patterns found in simulated TR-free sequences with those observed in sequences with TRs serves to disentangle TR-specific findings from those that may occur in

### **References**

genomic sequences in general. A simple approach is to simulate TR-free data by drawing *k*-mers from a (*k* - 1)-th order Markov model based on empirical frequencies (Robin et al., 2007). In contrast to drawing single characters from their frequency distribution, simulating *k*-mers mimics natural local correlations while choosing small *k* minimizes the chance of hidden TRs within a *k*-mer.

### **Conclusion and Perspectives**

Tandem repeats are diverse in their size, type, unit similarity, and distribution across genomes. Methods discussed above enable accurate detection of TR orthologs with strongly conserved unit configurations, or on the contrary, with highly changing unit numbers. Due to the heterogeneity of TRs, large-scale studies should be followed up by studies that focus on specific TR types and their effects on molecular processes. Our TR scans of eukaryotic proteomes provide a plentitude of cases to investigate with respect to their functional roles. While many TRs were linked to key functions, phenotypic changes, or disease predisposition, the biological mechanisms generating and preserving TRs in genomes are poorly understood. New studies of genomic TRs only fuel our fascination with these genomic features, calling for further research and for the development of dedicated methods. Further development of rigorous statistical models of TR generating mechanisms will help to improve TRD methods, and to shed some light on biological forces shaping these sequences.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Anisimova, Pe*ˇ*cerska and Schaper. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Accurate Prediction of the Statistics of Repetitions in Random Sequences: A Case Study in Archaea Genomes**

*Mireille Régnier 1,2 \* and Philippe Chassignet 2*

*1 Inria, Palaiseau, France, 2 LIX, Ecole Polytechnique, Palaiseau, France*

Repetitive patterns in genomic sequences have a great biological significance and also algorithmic implications. Analytic combinatorics allow to derive formula for the expected length of repetitions in a random sequence. Asymptotic results, which generalize previous works on a binary alphabet, are easily computable. Simulations on random sequences show their accuracy. As an application, the sample case of Archaea genomes illustrates how biological sequences may differ from random sequences.

**Keywords:** *K***-mers, combinatorics, probability**

## **1. INTRODUCTION**

### *Edited by:*

*Marco Pellegrini, Consiglio Nazionale delle Ricerche, Italy*

### *Reviewed by:*

*Travis Gagie, University of Helsinki, Finland Solon P. Pissis, King's College London, UK*

> *\*Correspondence: Mireille Régnier<EMAIL_ADDRESS>
### *Specialty section:*

*This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology*

> *Received: 03 December 2015 Accepted: 08 April 2016 Published: 08 June 2016*

### *Citation:*

*Régnier M and Chassignet P (2016) Accurate Prediction of the Statistics of Repetitions in Random Sequences: A Case Study in Archaea Genomes. Front. Bioeng. Biotechnol. 4:35. doi: 10.3389/fbioe.2016.00035* This paper provides combinatorial tools to distinguish biologically significant events from random repetitions in sequences. This is a key issue in several genomic problems as many repetitive structures can be found in genomes. One may cite microsatellites, retrotransposons, DNA transposons, long terminal repeats (LTR), long interspersed nuclear elements (LINE), ribosomal DNA, and short interspersed nuclear elements (SINE). In Treangen and Salzberg (2012), it is claimed that half of the genome consists of different types of repeats. Knowledge about the length of a maximal repeat is a key issue for assembly, notably the design of algorithms that rely upon de Bruijn graphs. In resequencing, it is a common assumption for aligners that any sequenced "read" should map to a single position in a genome: in the ideal case where no sequencing error occurs, this implies that the length of the reads is larger than the length of the maximal repetition. Average lengths of the repeats are given in Gu et al. (2000). Recently, heuristics have been proposed and implemented (Devillers and Schbath, 2012; Rizk et al., 2013; Chikhi and Medvedev, 2014).

A similar problem has been extensively studied: the prediction of the length of maximal common prefixes for words in a random set. Typical parameters are the background probability model, the size *V* of the alphabet, the length *n* of the sequence, and so on. Deviation from uniformity was studied for a uniform model as early as 1988 (Flajolet et al., 1988). A complexity index that captures the richness of the language is addressed in Janson et al. (2004). A distribution model, valid for binary alphabets and biased distributions, was introduced in Park et al. (2009), the so-called *trie profile* and extended to Patricia tries in Magner et al. (2014). The authors pointed out different "regimes" of randomness and a phase transition, by means of analytic combinatorics (Sedgewick and Flajolet, 2009). It was observed in Jacquet and Szpankowski (1994) that the average length of maximal common prefixes in a random set of *n* words is asymptotically equivalent to the average length of maximal repetitions in a random sequence of length *n*. Sets of words are considered below in the theoretical analysis. A comparison with the distribution of maximal repetitions in random sequences or real Archaea genomic sequences is presented in Section 3.

Our first goal is to extend results of Park et al. (2009) to the case of a general *V*-alphabet, including the special case {*A*, *C*, *G*, *T*} where *V* is 4. A second goal is to compare the results consistency with random data and real genomic data in the finite range.

To achieve the first goal, we rely on an alternative, and simpler, probabilistic and combinatorial approach that is interesting *per se*. It avoids generating functions and the Poissonization–dePoissonization cycle that is used in Park et al. (2009) and it extends to non-binary alphabets. In that case, there is no closed formula for the asymptotic behavior. Nevertheless, the Lagrange multipliers allow to derive it as the solution of an equation that can be computed numerically.

Explicit and computable bounds for the profile of a random set of *n* words are provided. Three domains can be observed. A first domain is defined by a threshold *k* for the length, called the *completion length*: any prefix with a length smaller than this threshold occurs at least twice. This threshold is extremely stable over the data sets and it is highly predictable. A similar phenomenon was observed for a uniform model in Fagin et al.(1979a) and a biased model (Mahmoud, 1992; Park et al., 2009). For larger lengths, some prefixes occur only once. In a second domain, called the *transition phase*, the number of maximal common prefixes is sublinear in the size *n* of the sequence: increasing first, then decreasing slowly, and, finally, dropping rapidly. In the third domain, for a length larger than some *extinction length*, almost no common prefix of that length occurs. Despite the fact that these bounds are asymptotic, a good convergence is shown in practice for random texts when a second-order term is known.

Differences between the model and the observation are studied on the special case of Archaea genomes. A dependency to the GCcontent, which is a characteristic of each genome, is exhibited. Regimes and transitions are studied on these genomic data and theoretical results are confirmed, with a drift in the values of transition thresholds. Notably, the length of the largest repetitions is much larger than expected. This difference between the model and the observation arises from the occurrences of long repeated regions.

Section 2 is devoted to Main Results, to be proved in Section 4. First, some notations are introduced; then, an algebraic expression for the expectation of the number of maximal common prefixes in a sequence is derived in Theorem 2.1. Second, this expression is split between two sums that are computable in practical ranges. Then, it is shown that a Large Deviation principle applies. It yields first and second order asymptotic terms, and oscillations, that are provided in Theorem 2.2. A comparison between exact, approximate, and asymptotic expressions is presented in Section 3.

### **2. MAIN RESULTS**

It is assumed throughout this study that sequences and words are randomly generated according to a biased Bernoulli model on an alphabet of size *V*. Let *p*1, *· · ·* , *pV* denote the probabilities of the *V* characters *χ*1, *· · ·* ,*χV*.

**Definition 2.1.** For any *i* in {1, *· · ·* ,*V*}, one notes

$$
\beta\_i = \log \frac{1}{p\_i}.
$$

Additionally,

$$p\_{\min} = \min\{p\_i; 1 \le i \le V\} \quad \text{and} \quad \alpha\_{\min} = \frac{1}{\log \frac{1}{p\_{\min}}} = \frac{1}{\max(\beta\_i)};\tag{1}$$

$$\mathfrak{p}\_{\max} = \max \{ \mathfrak{p}\_i; 1 \le i \le V \} \quad \text{and} \quad \alpha\_{\max} = \frac{1}{\log \frac{1}{\mathfrak{p}\_{\max}}} = \frac{1}{\min(\beta\_i)}. \tag{2}$$

The two values min(*βi*) and max(*βi*) are different when the Bernoulli model is non-uniform.

### **2.1. Enumeration**

**Definition 2.2.** Given *U* a set of words and an integer *k*, *k ≥* 2, a unique *k*-mer in *U* is a word *wχi* of length *k* such that

By convention, a unique 1-mer is a character *χi* that is a prefix of a single word.

**Definition 2.3.** Let *U* be a set of *n* words.

For *k ≥* 1, one denotes *B*(*n, k*) the number of unique *k*-mers in *U*.

One denotes *µ*(*n*, *k −* 1) the expectation of *B*(*n, k*) over all sets of *n* words.

**Remark:** It follows from Definition 2.2 that quantity *B*(*n, k*) is upper bounded by *n*. Observe that, for each random set *U*, it is the sum of a large number – *V k* – of correlated random variables. Expectation *µ*(*n, k*) is studied below and compared in Section 3 with *B*(*n, k* + 1).

Profiles of repetitions can be expressed as a combinatorial sum. **Theorem 2.1.** Given a length *k*, the expectation *µ*(*n, k*) satisfies:

$$\mu(n,k) = n \sum\_{k\_1+\cdots+k\_V=k} \binom{k}{k\_1,\cdots,k\_V} \phi(k\_1,\cdots,k\_V) \psi\_n(k\_1,\cdots,k\_V) \tag{3}$$

where

$$\phi(k\_1, \cdots, k\_V) = p\_1^{k\_1} \cdots p\_V^{k\_V} \tag{4}$$

$$\psi\_n(k\_1, \cdots, k\_V) = \sum\_{i=1}^V p\_i[(1 - \phi(k\_1, \cdots, k\_V)p\_i)^{n-1}]$$

$$\begin{aligned} \quad & = 1 \\ & \quad - \left( 1 - \phi(k\_1, \dots, k\_V) \right)^{n-1} \text{.} \end{aligned} \tag{5}$$

**Proof.** A word *wχi* is a unique (*k* + 1)-mer iff (i) *w* has length *k* and is the prefix of at least two words, including *wχi*; (ii) *wχi* is not repeated.

Event (i) has probability

$$n\phi(k\_1, \cdots, k\_V)p\_i[1 - (1 - \phi(k\_1, \cdots, k\_V))^{n-1}].$$

$$\text{Event (ii), which is a sub-event of (i), has probability }$$

$$n\phi(k\_1, \cdots, k\_V)p\_i[1 - (1 - \phi(k\_1, \cdots, k\_V)p\_i)^{n-1}].$$

**2.2. A Combinatorial Expression Definition 2.4.** Given a *k*-mer *w*, let *α* denote *k* log *n* and *ki* denote the number of occurrences of character *χi* in *w*. The *objective function* is

$$\rho(k\_1, \cdots, k\_V) = \sum\_{i=1}^{V} \frac{k\_i}{k} \beta\_i - \frac{1}{\alpha}. \tag{6}$$

The character distribution (*k*1, *· · ·* , *kV*) of a *k*-mer may be viewed as *barycentric coordinates* for a point *β*(*k*1*, · · · , kV*) = ∑*V i*=1 *ki k βi* that lies in [min(*βi*); max(*βi*)] = [ 1 *α*max ; 1 *α*min ] . The order of *β* points on that interval allows for a classification of *k*-mers that is a key to this study.

**Definition 2.5.** A *k*-mer *w* is said

**Remark:** If *ρ*(*k*1, *· · ·* , *kV*) = 0, the condition on the ancestor is trivially satisfied.

**Definition 2.6.** Given a set *U* of *n* words and an integer *k*, let *Dk(n)* denote the set of character distributions (*k*1, *· · ·* , *kV*) for rare and transition *k*-mers. Let *Ek*(*n*) denote the set of character distributions for common *k*-mers.

The set *Dk*(*n*) is the empty set if *k < αmin* log *n* and is the set of character distributions (*k*1, *· · ·* , *kV*) if *k > αmax* log *n*. Computation of (3) is split among the two sets *Dk*(*n*) and *Ek*(*n*). Computations show that the main contribution arises from transition *k*-mers. A probabilistic interpretation will be discussed in 2.4.

**Notation:** Let *S*(*k*) and *T*(*k*) be

$$S(k) = n \sum\_{D\_k(n)} \binom{k}{k\_1 \cdots k\_V} \phi(k\_1, \cdots, k\_V) \psi\_n(k\_1, \cdots, k\_V); \tag{7}$$

$$T(k) = n \sum\_{k\_k(n)} \binom{k}{k\_1 \cdots k\_V} \phi(k\_1, \cdots, k\_V) \psi\_n(k\_1, \cdots, k\_V). \tag{8}$$

So *µ*(*n, k*) rewrites

$$
\mu(n,k) = \mathcal{S}(k) + T(k). \tag{9}
$$

These sums *S*(*k*) and *T*(*k*) can be efficiently computed for moderate *k*, up to a few hundred, approximately. In practice, *αmax* log *n* is below this threshold for the sizes of actual genomes and for their ordinary GC content value. The simulations in Section 3 show that this estimation is rather tight. Behavior and asymptotic estimates are derived and discussed in the next section.

### **2.3. Asymptotic Estimates**

In this section, asymptotic estimates for (3) are derived. First, some characteristic functions are introduced. Then, it is observed that a Large Deviation Principle applies for the combinatorial sums to be computed and asymptotics for the dominating term follow. Amortized terms are also computed. It is shown in Section 3 that this second-order term cannot be neglected in the finite range.

### 2.3.1. Notations

For general alphabets, asymptotic behavior is a function of the solution of an equation and depends on domains whose bounds are defined below.

**Definition 2.7.** Let (*pi*)1*≤i≤V* be a Bernoulli probability distribution. Let *σ*2 denote ∑*V i*=1 *p* 2 *i* .

The *fundamental ratio*, noted *α*˜, is (∑ *i pi* log 1 *pi* )*−*1 .

The transition ratio, noted  $\bar{\alpha}$ , is  $\sigma\_2 \left(\sum\_i p\_i^2 \log \frac{1}{p\_i}\right)^{-1}$ . The extinction ratio, noted  $\alpha\_{ext}$ , is  $\frac{2}{\log \frac{1}{\sigma\_2}}$ .

**Definition 2.8.** Let *α* be a real value in [*αmin*, *αmax*]. Let *τ α* be the unique real root of the equation

$$\frac{1}{\alpha} = \frac{\sum\_{i=1}^{V} \beta\_i e^{-\beta\_i \tau}}{\sum\_{i=1}^{V} e^{-\beta\_i \tau}}\tag{10}$$

Let *ψ* be the function defined in [*αmin*, *αext*] as

$$\begin{aligned} \alpha\_{\min} \le \alpha \le \bar{\alpha} : \psi(\alpha) &= \tau\_{\alpha} + \alpha \log \left( \sum\_{i=1}^{V} e^{-\beta\_i \tau\_{\alpha}} \right); \\\bar{\alpha} \le \alpha : \psi(\alpha) &= 2 - \alpha \log \frac{1}{\sigma\_2}. \end{aligned}$$

**Proposition 2.1.** The following property holds

$$
\alpha\_{\min} \le \tilde{\alpha} \le \bar{\alpha} \le \alpha\_{\max} \le \alpha\_{\text{ext}}.
$$

Function *ψ* increases on [*αmin, α*˜] and decreases on [˜*α, ∞*]. It satisfies

$$
\psi(\alpha\_{\min}) = \psi(\alpha\_{\text{ext}}) = \mathbf{0} \text{ and } \psi(\tilde{\alpha}) = 1. \tag{11}
$$

**Remark:** Uniqueness of *τ α* is shown in Section 4.2. As *τα*¯ = 2, *ψ* is continuous at *α* = ¯*α*, with *ψ*(¯*α*) = 2 *− α*¯ log 1 *σ*2 .

#### 2.3.2. Asymptotic Results

**Theorem 2.2.** Given a length *α* log *n*, when *n* tends to *∞* the ratio log *µ*(*n, α* log *n*) log *n* satisfies:

$$0 \le \alpha \le \alpha\_{\min} \text{ or } \alpha\_{\text{ext}} \le \alpha : \frac{\log \mu(n, \alpha \log n)}{\log n} \le 0;\tag{12}$$

$$\alpha\_{\min} \le \alpha \le \alpha\_{\text{ext}} : \frac{\log \mu(n, \alpha \log n)}{\log n} \sim \psi(\alpha). \tag{13}$$

Moreover, let *ξ* be the function defined in [*αmin*,*αext*] as *ξ*(*α*) = *µ*(*n, α* log *n*) log *n − ψ*(*α*). It satisfies

$$\alpha\_{\min} \le \alpha \le \bar{\alpha} : \xi(\alpha) \sim -\frac{V-1}{2} \frac{\log(\alpha \log n)}{\log n}; \qquad \text{(14)}$$

$$\bar{\alpha} \le \alpha \le \alpha\_{\text{ext}} : \xi(\alpha) \sim \frac{\log(1 - \sigma\_2)}{\log n}. \tag{15}$$

**Proof**. The key to the proof when*α*ranges in [*αmin*,*αmax*] is that *ψn*(*k*1, *· · · kV*) is maximal when *ρ*(*k*1, *· · · kV*) is close to 0. Sum *T*(*k*) satisfies a Large Deviation Principle.

$$\frac{\log T(\tilde{k})}{k} \sim \max \left\{ -\sum\_{i=1}^{V} \frac{k\_i}{k} \log \frac{k\_i}{k}; \rho(k\_1, \dots, k\_V) = 0 \right\}. \tag{16}$$

The maximization problem rewrites as

$$\max \left\{ \sum\_{i=1}^{V} \theta\_i \log \frac{1}{\theta\_i}; \sum\_{i=1}^{V} \theta\_i = 1; \sum\_{i=1}^{V} \beta\_i \theta\_i = \frac{1}{\alpha}; 0 \le \theta\_i \le 1 \right\} \tag{17}$$

The maximum value is *τα* + *α* log (∑*V i*=1 *e −βiτα* ) that is reached for the *V*-tuple ( *θi* = *e −βiτα* ∑*V i*=1 *e−βiτα* ) 1*≤i≤V* .

*S*(*k*) satisfies again a Large Deviation Principle when *α < α*¯, which yields the asymptotic result in this range. For larger *α*, *S*(*k*) is approximately (1 *− σ*2)*n* 1*−α* log 1 *σ*2 that dominates *T*(*k*).

Details for the proof, including the short and long lengths, are provided in Section 4.

**Remark:** The discussion will depend of the ratio *α* = *k* log *n* . Possible values for *α* range over a *discrete* set as they are constrained to be the ratio of an integer by the log of an integer. An interesting property is that, for any real *α*, the set *T* = {*n ∈ N*; *α* log *n ∈ N*} is either empty or infinite. Indeed, when *T* is non-empty, it contains all values *n*(*α*) *p* where *n*(*α*) denotes the minimum value of *T*. It is beyond the scope of this paper to establish the number of other possible solutions.

### 2.3.3. Domains

Different domains arise from this Theorem, which were observed in Park et al. (2009). Equalities *ψ*(*αmin*) = 0 and *ψ*(¯*α*) = 2 *− α*¯*log* 1 *σ*2 show that there is a continuity between domains.

When *α* lies inside the domain [*αmin*, *αext*], the ratio log *µ*(*n,α* log *n*) log *n* is positive and parameters *µ*(*n*, *α* log *n*) are *sublinear* in the size *n* of the text: some *k*-mers – mostly transition *k*-mers – are unique *k*-mers. Observe that the maximum value for *ψ*(*α*) is 1. When the Bernoulli model is uniform, this central domain is empty.

When the length is smaller than the *completion length αmin* log *n* or greater than the *extinction length αext* log *n*, the ratio log *µ*(*n, α* log *n*) log *n* is negative.

### 2.3.4. Oscillations

Parameters (*k*1, *· · ·* , *kV*) in the combinatorial sums are integers. As the optimum values (*kθi*)1*≤i≤V* may not be integers, the practical maximum is a close point on the lattice (*k*1, *· · ·* , *kV*). The difference introduces a multiplicative factor that ranges in [ *−* log *pmax pmin ,* log *pmax pmin* ] . This leads to a small *oscillation* of log *µ*(*n, k*). For large *n*, this contribution to log *µ*(*n,k*) log *n* becomes negligible. As mentioned above, the set of lengths *n* that are *admissible* for a given *α* is very sparse. Nevertheless, an approximate value may be used: for instance, for an integer *k ′* , 1 *k ′* log ⌈ *n*(*α*) *k ′ k* ⌉ is very close to *α*. This oscillation phenomenon was first observed in Nicodème (2005).

### 2.3.5. Binary Alphabets

Results for binary alphabets in Park et al. (2009) steadily follow from Theorem 2.2. A rewriting of *ψ* leads to alternative expression (18). This *explicit* expression points out the dependency to the distances to *αmin* and *αmax*, and the behavior around these points.

**Corollary 2.1.** Assume that the alphabet is binary. Then

$$\psi(\alpha) = \frac{\alpha}{\log \frac{p\_{\text{max}}}{p\_{\text{min}}}} \log \left[ s\_{\alpha}^{-\frac{1}{\alpha} - \frac{1}{\alpha\_{\text{min}}}} + s\_{\alpha}^{-\frac{1}{\alpha} - \frac{1}{\alpha\_{\text{max}}}} \right] \tag{18}$$

where

$$s\_{\alpha} = \frac{\alpha\_{\min}}{\alpha\_{\max}} \cdot \frac{\alpha - \alpha\_{\min}}{\alpha\_{\max} - \alpha}. \tag{19}$$

A similar result holds for DNA sequences when the alphabet is 4-ary and the probability distribution satisfies *pA* = *pT* and *pC* = *pG*. Such a distribution is defined by its GC-content *pG* + *pC*.

### **2.4. A Probabilistic Interpretation**

The main contribution to *µ*(*n, k*) arises from *k*-mers with an objective function close to 0, i.e., transition *k*-mers. Such *k*-mers exist in the *transition phase* [*αmin* log *n*, *αmax* log *n*] where they coexist with rare or common *k*-mers. Observe that this phase is *shrinked* when the Bernoulli model is uniform, as *pmin* = *pmax* and *αmin* = *αmax*. Therefore, most unique *k*-mers are concentrated on the two lengths *⌊αmin* log *n⌋* and *⌈αmin* log *n⌉*, as observed initially in Fagin et al. (1979b).

Let *k* be some integer in the transition phase. First, the relative contribution of *S*(*k*) and *T*(*k*) to *µ*(*n, k*) varies with the length *k*. For lengths close to *αmin* log *n*, most words are common and *T*(*k*) dominates *S*(*k*). When *k* increases, the proportion of common words decreases and the relative contribution of *T*(*k*) decreases.

Second, the dominating term in *µ*(*n, k*) arises from transition *k*-mers. Let *w* be a word of length *k*, the character distribution in *w* be (*k*1, *· · ·* , *kV*) and *χi* be some character. The number of words that admit *w* or *wχi* as a prefix fluctuates around the expectations *nϕ*(*k*1, *· · ·* , *kV*) and *nϕ*(*k*1, *· · ·* , *kV*)*pi*, respectively. On the one hand, when word *wχi* is a rare word, *nϕ*(*k*1, *· · ·* , *kV*) is less than 1. The smallest *nϕ*(*k*1, *· · ·* , *kV*) is, the less likely the actual number of occurrences of *w* is greater than 2 and the smallest the contribution of *wχi* to *S*(*k*), and *µ*(*n, k*), is. On the other hand, let *wχi* be a common *k* + 1-mer; *w* is a common *k*-mer and then *nϕ*(*k*1, *· · ·* , *kV*) is greater than 1. The largest *nϕ*(*k*1, *· · ·* , *kV*) is, the more likely the word *wχi* is repeated and the smallest the contribution to *T*(*k*), and *µ*(*n, k*), is.

For a short length, i.e., *k* smaller than the completion length *kmin*, all words are common. In a given sequence, most *k*-mers are repeated at least twice and there is (almost) no unique *k*-mers.

For a large length *k*, i.e., *k* greater than *kmax*, all words are rare. Nevertheless the number of unique *k*-mers remains sublinear in *n* in the range [*αmax* log *n*, *αext* log *n*]: the sum of small contributions arising from a large number of possible words is significant.

A folk theorem (Szpankowski, 2001; Jacquet and Szpankowski, 2015) claims that the objective function is concentrated around 1 *α*˜ *−* 1 *α* . Consequently, when *α* = ˜*α*, most *k*-mers are transition *k*-mers and the exponent, the *ψ* function, is maximal.

### **3. EXPERIMENTS AND ANALYSIS**

Simulations are presented for random and real data. For each simulation, a suffix tree (Ukkonen, 1995) is built, where each leaf represents a unique *k*-mer. For random cases, the Ukkonen's insertion step is iterated until a tree with exactly *n* leaves is build. This requires *n* + *kins* insertions of symbols, where *kins >* 0 is relatively small (there is a value of a few dozen in practice for considered *n*). One can observe that the event of having *n* leaves after *n* + *k −* 1 insertions corresponds to the fact that the trailing *k*-mer is unique in the sequence of length *n* + *k −* 1.

Even if a statistical bias exists, with respect to the case of a set of N random words analyzed in previous sections, this bias for respective values on *k* and *n* is below the numeric precision used for tables below.

Then, one simulation that is related to the case of a set of *n* random words, requires the generation of the order of N random symbols from a small alphabet, following a Bernoulli scheme. For this range of *n*, and even in the case of a hundred consecutive simulations, this corresponds to a regular use of a common random number generator (Knuth, 1998).

A first set of simulation deals with the case of random sequences over a binary alphabet, since the results can be compared with previous work. A second set addresses the case of random sequences over a quaternary alphabet {*A*, *C*, *G*, *T*} with a constrained distribution such that probabilities *pA ≈ pT* and *pC ≈ pG* as it is the case for DNA sequences (where the sum *pC* + *pG* is also known as the GC-content). Results on such random sequences are then compared with the sample biological sequence of an Archaea (*Haloferax volcanii*).

An implementation with a suffix array (Manber and Myers, 1993) allows for a compact representation and an efficient counting (Beller et al., 2013).

### **3.1. Random data**

A hundred binary sequences were randomly generated. The number of leaves in each tree was fixed to *n* = 5000000 and the Bernoulli parameter was *pmax* = 0.7000. Therefore, *pmin* = 0.3000, ˜*p* = 0.5429, and log *n* = 15.4249. The thresholds for *α* and the corresponding lengths *α* log *n* are:

### 3.1.1. Statistical Behavior on Random Sets

Throughout experiments, every sample profile for a given sequence fluctuates very little around the expectation.

**Table 1** provides experimental results averaged over a hundred binary sequences. Short length with no observed unique *k*-mer is removed. Column 2 gives the mean of *B*(*k* + 1), i.e., the mean number of observed leaves at depth *k* + 1, over the set of a hundred simulations. Columns 3 to 5 give the computed values for *S*(*k*), *T*(*k*), and *µ*(*k*), using the expressions, equations (7–9).

The actual number of leaves *B*(*n*, *k* + 1) is very close to the average value *µ*(*n, k*), and simulations show that this is the general case when (only) a hundred simulations are performed: *µ*(*n, k*) is a very good prediction.

Observed lengths of extinction also show very little variations. In array below, each column gives *nk*, the number of sequences out of the one hundred sample set for which the longest repetition had length *k*.

#### **Distribution of the extinction level for 100 random binary sequences.** *pmax* **is 0.7.**

In the binary case, the predicted extinction length is between 56 and 57. It is noticeable that, in most cases, the observed depth is slightly smaller than this value. In **Table 1**, value 0.04 for *µ*(*n*, 61) means that one expects a total of four leaves at depth 60 over one hundred sequences. In that run, exists a total amount of 8.

### 3.1.2. Quality of Estimates

### 3.1.3. Dependency to Probability Bias

Thresholds were computed for a given sequence length *n* and various probabilities. The more *pmax* departs from 0.5, the value for the uniform model, the largest the extinction length is. The completion length, *kmin*, slightly decreases, while the extinction length significantly increases. Nevertheless, this effect is limited when the largest probability *pmax* remains in the range [0.5;0.7].

### **3.2. Long Repetitions in Archaea Genomes**

The experimental data set is the sequence from *Haloferax volcanii DS2 chromosome, complete genome* (Hartman et al., 2010). The alphabet is quaternary. Profile results are shown in **Table 2**.

Sequence length is *n* = 2847757. The observed symbol frequencies are *pA* = 0.1655; *pC* = 0.3334; *pG* = 0.3330; *pT* = 0.1681. Therefore, observed *GC-content* is 0.6664. Parameters for an approximate degenerated quaternary model are *pA* = *pT* = *pmin* = 0.1668; *pC* = *pG* = *pmax* = 0.3332; ˜*p* = 0*.*2645; and log *n* = 14.8620. The thresholds for the domain are

Statistics on one hundred random sequences with same parameters are shown in **Table 3**. GC-content is 0.6664. Extinction level is provided in **Table 4**. Observe first a good match between the observed values, the predicted values for *µ*(*n, k*), and the asymptotic values for random data. As shown for binary alphabets,

#### **TABLE 1 | Mean profile for 100 random binary sequences**.

*(pmax; pmin)* = *(0.7; 0.3).*

**TABLE 2 | Profile for the sequence from** *Haloferax volcanii* **DS2 chromosome, complete genome**.

the observed extinction level for random sequences departs very little from the predicted *kext* level.

Numerous differences with random data can be observed on real genomes.

Interestingly, the behavior for short lengths and in the transition phase is similar to the random behavior. Observation and prediction have the same order of magnitude. In particular, the number of unique *k*-mers is maximum for length ˜*k* where observation and prediction coincide. For a real genome and a length *k* smaller than *kmin*, observed *B*(*n*, *k* + 1) is larger than predicted *µ*(*n, k*). This indicates, at a level *k* + 1 where completion is expected, more leaves in the real trie, more missing words at level *k* + 2. Simultaneously, less internal nodes occur at level *k* + 1 because the total sum is constant and equal to *V k*+1 .

The effect of (non-random) repetitions is more sensible in the decreasing domain. First, the number of unique *k*-mers decreases much more slowly than expected for lengths larger than *kmax*. A significant gap can be observed around extinction level *kext*. The decrease rate, which was around 0.02–0.04 drops to 0.007 and then becomes even lower. Finally, the extinction level is much larger than the predicted value 23: the largest repetition is 1395 bp long.

To evaluate the contribution of long repetitions, one may erase the longest ones. When a word *w* is repeated, any proper suffix of *w* is also repeated. Consequently, once the longest repeated word is erased, one unique *k*-mer (only) disappears for each length larger than the length of the second largest subsequence (here, 935). The profile remains far from the random profile. This observation is still true if the 10 longest subsequences are erased.

### **4. COMBINATORIAL AND ANALYTIC DERIVATION**

### **4.1. Lagrange Multipliers**

Lagrange multipliers method allows to maximize an expression under constraints. To compute (17), one sets

$$F = \sum\_{i=1}^{V} \theta\_i \log \theta\_i;\tag{20}$$

$$G = \sum\_{i=1}^{V} \theta\_i;\tag{21}$$

$$H = \sum\_{i=1}^{V} \theta\_i \beta\_i. \tag{22}$$

Two constraints are given:

$$G = 1 \text{ and } H = \frac{1}{\alpha}.$$

An intermediary function *ϕα*(*τ* 1, *· · · τV*) is defined

$$
\phi\_{\alpha} = F + \lambda\_{\alpha} G + \tau\_{\alpha} H \tag{23}
$$

In order to maximize *ϕ* under these two constraints, *ϕ* function is derived with respect to each random variable *τ i*. This yields *V* equations

$$1 + \log \theta\_i + \lambda\_\alpha + \tau\_\alpha \beta\_i = 0. \tag{24}$$

Two indices *imin* and *imax* are chosen that satisfy *βimin ̸*= *βimax* . For instance

$$\begin{aligned} \beta\_{i\_{\min}} &= \min \left( \beta\_i \right)\_{1 \le i \le V} = \log \frac{1}{p\_{\max}}; \\ \beta\_{i\_{\max}} &= \max \left( \beta\_i \right)\_{1 \le i \le V} = \log \frac{1}{p\_{\min}}. \end{aligned}$$

Solving equation (24) with indices *imin* and *imax* yields

$$\begin{aligned} \tau\_{\alpha} &= \frac{\log \theta\_{i\_{\rm{min}}} - \log \theta\_{i\_{\rm{max}}}}{\beta\_{i\_{\rm{max}}} - \beta\_{i\_{\rm{min}}}} = \log \frac{\theta\_{i\_{\rm{min}}} \frac{1}{\beta\_{i\_{\rm{max}}} - \beta\_{i\_{\rm{min}}}}}{\theta\_{i\_{\rm{max}}}};\\ 1 + \lambda\_{\alpha} &= \frac{\beta\_{i\_{\rm{min}}} \log \theta\_{i\_{\rm{max}}} - \beta\_{i\_{\rm{max}}} \log \theta\_{i\_{\rm{min}}}}{\beta\_{i\_{\rm{max}}} - \beta\_{i\_{\rm{min}}}}. \end{aligned}$$

Remaining equations rewrite:

$$\log \theta\_{\dot{\imath}} = \log \theta\_{\dot{\imath}\_{\text{min}}} + \tau\_{\alpha} (\beta\_{\dot{\imath}\_{\text{min}}} - \beta\_{\dot{\imath}}).\tag{25}$$

Using the constraint ∑*V i*=1 *θi* = 1 that yields

$$
\theta\_{i\_{\min}} e^{\beta\_{i\_{\min}} \tau\_{\alpha}} \sum\_{i=1}^{V} e^{-\beta\_{i} \tau\_{\alpha}} = 1,
$$

*GC-content is 0.6664.*

#### **TABLE 4 | Distribution of the extinction level for 100 random degenerated quaternary sequences**.

*GC-content is 0.6664.*

and an expression for *θimin* follows. Therefore Equation 25 rewrites:

$$\theta\_i = \frac{e^{-\beta\_i \tau\_\alpha}}{\sum\_{i=1}^{V} \beta\_i e^{-\beta\_i \tau\_\alpha}}. \tag{26}$$

Finally, Equation ∑*V i*=1 *θiβi* = 1 *α* yields equation (10).

$$\frac{1}{\alpha} = \frac{\sum\_{i=1}^{V} \beta\_i e^{-\beta\_i \tau\_\alpha}}{\sum\_{i=1}^{V} e^{-\beta\_i \tau\_\alpha}}.$$

For this *V*-tuple

$$\begin{split} \sum\_{i=1}^{V} \theta\_i \log \theta\_i &= -\left(\sum\_{i=1}^{V} \theta\_i \beta\_i\right) \tau\_\alpha - \left(\sum\_{i=1}^{V} \theta\_i\right) \log \left(\sum\_{i=1}^{V} e^{-\beta\_i \tau\_\alpha}\right) \\ &= -\frac{\tau\_\alpha}{\alpha} - \log \left(\sum\_{i=1}^{V} e^{-\beta\_i \tau\_\alpha}\right). \end{split}$$

### **4.2. Approximation Orders**

Derivating the RHS of (10) yields ∑ *i̸*=*j* (*βi*+*βj*) 2 *e −*(*βi*+*βj* )*τ* ( ∑ *i e−βiτ* ) 2 that is positive. Therefore, for any *α*, the solution to (10) is unique. Moreover, *τ α* increases with *α*. Let

$$\psi\_1(\alpha) = \tau\_\alpha + \alpha \log \left( \sum\_{i=1}^V e^{-\beta\_i \tau\_\alpha} \right);\tag{27}$$

$$
\psi\_2(\alpha) = 2 - \alpha \log \frac{1}{\sigma\_2}.\tag{28}
$$

Notably, the solutions *τ α* of (10) associated with the four increasing values of *α*: (*αmin, α,* ˜ *α, α* ¯ *max*) are (–*∞*, 1 + 2, +*∞*). Computing *ψ* for these values yields (11) and Equality *ψ*1(˜*α*) = *ψ*2(˜*α*).

Derivating both expressions yields

$$\frac{\partial \psi\_1}{\partial \alpha}(\alpha) = \log \left( \sum\_{i=1}^{V} e^{-\beta\_i \tau\_\alpha} \right);\tag{29}$$

$$\frac{\partial \psi\_1}{\partial \alpha}(\alpha) - \frac{\partial \psi\_2}{\partial \alpha}(\alpha) = \log \left( \frac{1}{\sigma\_2} \sum\_{i=1}^{V} e^{-\beta\_i \tau\_\alpha} \right). \tag{30}$$

Both derivatives are monotone functions of *τ α*. In equation (30), derivative is 0 when *α* = ¯*α*. Therefore, *ψ* is the maximum of the two values *ψ*1 and *ψ*2 over the interval [*αmin*, *αmax*]. The former equation is 0 if *α* = ˜*α*. Therefore, *ψ* is maximum when *α* = ˜*α*.

### **4.3. Approximations** 4.3.1. Short Lengths

Assume that *k ≤ αmin* log *n*. Each term *ϕ*(*k*1, *· · ·* , *kV*) is lower bounded by *p k min* = *n αlogpmin* = *n − α αmin* . Each term *ψn*(*k*1, *· · ·* , *kV*) is trivially bounded by *e −n* 1*− α αmin* that is upper bounded by 1 and *nψn*(*k*1, *· · ·* , *kV*) tends to 0 when *n* goes to *∞*. As ∑( *k k*1*···kV* ) *ϕ*(*k*1*, · · · , kV*) = 1, the ratio log *µ*(*n, k*) log *n* is negative.

### 4.3.2. Moderate and Large Lengths

For a length *k* in the transition domain [*αmin* log *n*, *αmax* log *n*], the objective function may be either positive or negative. When *k > αmax* log *n*, set *Ek*(*n*) is empty and *µ*(*n, k*) reduces to *S*(*k*).

The maximum *M* among the terms *e k* ( *−* ∑ *i ki k* log *ki k −* 1 *k* log *nϕ*(*k*1*,··· ,kV*) ) in *T*(*k*) is reached when *ρ*(*k*1, *· · ·* , *kV*) is 0. Due to the exponential decrease of *e −nϕ*(*k*1*,··· ,kV*) when *nϕ*(*k*1, *· · ·* , *kV*) *≥* 1, *T*(*k*) *k* is upper bounded. Computation of log *M* is done with Lagrange multipliers, as explained above.

Computation of *S*(*k*) relies on the local development of *ψn*(*k*1, *· · ·* , *kV*), that is *n*(1–*σ*2)*ϕ* (*k*1, *· · ·* , *kV*). *S*(*k*) rewrites *σ*2 *k* ˜*S*(*k*) + (*S*(*k*) *− σ*2 *k* ˜*S*(*k*)) where ˜*S*(*k*) = ∑ ( *ρ*(*k*1*,··· ,kV*)*≤*0 *k ki* ) ( *p* 2 1 *σ*2 )*k*1 *· · ·* ( *p* 2 *V σ*2 )*kV* . This sum satisfies a Large Deviation Principle when *ρ*(*k*1*, · · · , kV*) + 1 *α ≥* 1 *α*˜ , or *α < α*˜. In this range, ˜*S*(*k*) *k ∼* max { *−* ∑*V i*=1 *ki k* log *ki k* } , which was shown to be *ψ*(*α*).

When *α > α*˜, sum ˜*S*(*k*) rewrites 1 *−* ¯*S*(*k*) where

$$\bar{S}(k) = \sum\_{\rho(k\_1,\dots,k\_V) + \frac{1}{\alpha} < \frac{1}{\alpha}} \binom{k}{k\_l} \left(\frac{p\_1^2}{\sigma\_2}\right)^{k\_1} \cdots \left(\frac{p\_V^2}{\sigma\_2}\right)^{j\_V}.$$

This sum satisfies a Large Deviation Principle and

$$\frac{\bar{\mathbf{S}}(k)}{k} \sim \max \left\{ -\sum\_{i=1}^{V} \frac{k\_i}{k} \log \frac{k\_i}{k} + \sum\_{i=1}^{V} \frac{k\_i}{k} \log \frac{p\_i^2}{\sigma\_2} \right\}.$$

$$\text{As } \sum\_{i=1}^{V} \frac{k\_i}{k} \log \frac{p\_i^2}{\sigma\_2} = -\frac{2}{\alpha} + \log \frac{1}{\sigma\_2}, \text{ this maximum is}$$

$$-\frac{1}{\alpha} [2 - \alpha \log \frac{1}{\sigma\_2} - \psi(\alpha)]$$

that is negative.

### **4.4. Binary Case**

Barycentric coordinates of *α* are unique. Indeed, equation (10) reduces to a linear equation on the variable *e −*(*β*2*−β*1)*τ*

$$\frac{1}{\alpha} = \frac{\beta\_1 + \beta\_2 e^{-(\beta\_2 - \beta\_1)\tau}}{1 + e^{-(\beta\_2 - \beta\_1)\tau}}$$

where *β*2*−β*1 = *βmin−βmax* = log *pmax pmin* . Therefore, *e −*(*β*2*−β*1)*τ* = 1*−αβ*1 *αβ*2*−*1 . Finally

### **REFERENCES**

$$\tau\_{\alpha} = \frac{1}{\log \frac{p\_{\text{max}}}{p\_{\text{min}}}} \log \frac{\alpha \beta\_2 - 1}{1 - \alpha \beta\_1} = \frac{1}{\log \frac{p\_{\text{max}}}{p\_{\text{min}}}} \log \frac{\frac{1}{\alpha\_{\text{min}}} - \frac{1}{\alpha}}{\frac{1}{\alpha} - \frac{1}{\alpha\_{\text{max}}}}.$$

Function *ψ* rewrites, in the binary case:

$$
\psi\_{\alpha} = \tau\_{\alpha} = \alpha \log e^{-\frac{1}{\alpha}\tau\_{\alpha}} \left( e^{-\left(\beta\_1 - \frac{1}{\alpha}\right)\tau\_{\alpha}} + e^{-\left(\beta\_2 - \frac{1}{\alpha}\right)\tau\_{\alpha}} \right) .
$$

Observing that *e −*(*β*2*−β*2)*τα* = *sα* and changing variable *τ α* into (*β*2 *− β*1) yields *e −* ( *β*1*−* 1 *α* ) *τα* = *sα −* ( 1 *αmin −* 1 *α* ) and *e −* ( *β*2*−* 1 *α* ) *τα* = *sα −* ( 1 *αmax −* 1 *α* ) .

### **5. CONCLUSION**

This paper describes the behavior of the number of unique or repeated *k*-mers in a random sequence, on a general alphabet. Derivation relies on a combination of analytic combinatorics and on Lagrange multipliers. It simplifies an approach provided for binary alphabets and allows to address larger alphabets, including the quaternary alphabets, such as DNA alphabet. Precise asymptotic estimates are provided and a probabilistic interpretation is given. They are validated on random simulated data and shown to be valid in the finite range. Therefore, they provide a valuable tool to estimate a suitable read length for assembly purposes and tune parameters for assembly algorithms. Real genomes significantly depart from the random behavior for long repetitions. The general shape of the trie profile is observed, with a maximum of the number of unique *k*-mers at the expected length. However, for real genomes, a number of very short *k*-mers are missing and, on the contrary, one observes a number of very long repetitions. Besides these events, the behaviors are rather similar.

In the future, it is worth extending the method to generalized Patricia tries, Markov models and approximate repetitions.

### **AUTHOR CONTRIBUTIONS**

Both authors contributed equally.

### **FUNDING**

Inria-Cnrs-Poncelet grant Carnage.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2016 Régnier and Chassignet. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Transposable element insertions in long intergenic non-coding RNA genes

*Sivakumar Kannan1†, Diana Chernikova2†, Igor B. Rogozin1†, Eugenia Poliakov 3 , David Managadze1 , Eugene V. Koonin1 and Luciano Milanesi 4 \**

*1National Center for Biotechnology Information, National Library of Medicine, National Institutes of Health, Bethesda, MD, USA, 2Department of Genetics, Institute for Quantitative Biomedical Sciences, Geisel School of Medicine, Dartmouth College, Hanover, NH, USA, 3 Laboratory of Retinal Cell and Molecular Biology, National Eye Institute, National Institutes of Health, Bethesda, MD, USA, 4 Institute for Biomedical Technologies, National Research Council, Segrate, Italy*

### *Edited by:*

*Marco Pellegrini, Consiglio Nazionale delle Ricerche, Italy*

*Reviewed by:*

*Roderic Guigo, Center for Genomic Regulation, Spain Justin Blumenstiel, University of Kansas, USA*

#### *\*Correspondence:*

 *Luciano Milanesi, Institute for Biomedical Technologies, National Research Council, Via Fratelli Cervi 93, 20090 Segrate, Italy<EMAIL_ADDRESS>
*† Sivakumar Kannan, Diana Chernikova and Igor B. Rogozin have contributed equally to this work.*

#### *Specialty section:*

*This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology*

*Received: 10 December 2014 Accepted: 06 May 2015 Published: 09 June 2015*

#### *Citation:*

*Kannan S, Chernikova D, Rogozin IB, Poliakov E, Managadze D, Koonin EV and Milanesi L (2015) Transposable element insertions in long intergenic non-coding RNA genes. Front. Bioeng. Biotechnol. 3:71. doi: 10.3389/fbioe.2015.00071*

Transposable elements (TEs) are abundant in mammalian genomes and appear to have contributed to the evolution of their hosts by providing novel regulatory or coding sequences. We analyzed different regions of long intergenic non-coding RNA (lincRNA) genes in human and mouse genomes to systematically assess the potential contribution of TEs to the evolution of the structure and regulation of expression of lincRNA genes. Introns of lincRNA genes contain the highest percentage of TE-derived sequences (TES), followed by exons and then promoter regions although the density of TEs is not significantly different between exons and promoters. Higher frequencies of ancient TEs in promoters and exons compared to introns implies that many lincRNA genes emerged before the split of primates and rodents. The content of TES in lincRNA genes is substantially higher than that in protein-coding genes, especially in exons and promoter regions. A significant positive correlation was detected between the content of TEs and evolutionary rate of lincRNAs indicating that inserted TEs are preferentially fixed in fast-evolving lincRNA genes. These results are consistent with the repeat insertion domains of LncRNAs hypothesis under which TEs have substantially contributed to the origin, evolution, and, in particular, fast functional diversification, of lincRNA genes.

Keywords: mobile elements, molecular domestication, exaptation, junk DNA, long non-coding RNA, repetitive elements

## Introduction

Traditionally, genomes have been perceived mostly as repositories of protein-coding genes. Although this might be largely true in the case of viruses, prokaryotes, and unicellular eukaryotes, numerous recent studies on the genomes of multicellular eukaryotes, particularly animals, have revealed a vast non-coding RNome, i.e., numerous genes encoding various classes of non-coding RNAs (ncRNAs) (Carninci et al., 2005; Mattick and Makunin, 2006; Ponting et al., 2009; Derrien et al., 2012; Amaral et al., 2013). Strikingly, the total number of genes for ncRNAs that are expressed from a mammalian genome seems to exceed the number of protein-coding genes several fold (Mattick and Makunin, 2006; Amaral et al., 2013). The classification of ncRNAs and validation of their functionality remain matters of intensive investigation and debate (Van Bakel and Hughes, 2009; Ponting and Belgard, 2010; Graur et al., 2013). Among many distinct classes of ncRNAs, the long non-coding RNA (lncRNA) is probably the most enigmatic group. The definition of a lncRNA is based solely on the transcript size: lncRNAs are defined as ncRNAs longer than 200 nt (Mattick and Makunin, 2006; Ponting et al., 2009). Many lncRNAs are spliced, 5′capped, and polyadenylated (Okazaki et al., 2002; Carninci et al., 2005; Kapranov et al., 2007; Ponjavic et al., 2007). Based on the localization in the genome, lncRNAs can be divided into two distinct classes: (i) transcripts that overlap protein-coding genes, many of which are likely to be involved in sense–antisense regulation (Chen et al., 2005; Ponting and Belgard, 2010; Rinn and Chang, 2012) and (ii) long intergenic non-coding (linc)RNAs that are transcribed from genome regions separating protein-coding genes (Ponjavic et al., 2007; Mercer et al., 2008; Ponting et al., 2009).

The current knowledge on the functions of long intergenic non-coding RNAs (lincRNAs) is scarce because very few of the lincRNAs have been experimentally characterized. Nevertheless, the functional range of this class of ncRNA is believed to be broad on the basis of indirect evidence (Bertone et al., 2004; Ponjavic et al., 2007; Mercer et al., 2008; Ponting and Belgard, 2010; Ulitsky et al., 2011; Glazko et al., 2012; Ng et al., 2013). It has been proposed that lincRNAs could be involved in the regulation of many cellular processes (Mattick and Makunin, 2006; Loewer et al., 2010; Wang et al., 2011; Rinn and Chang, 2012). For example, they can affect transcription locally on the gene level (Martens et al., 2004; Martianov et al., 2007; Osato et al., 2007; Hirota et al., 2008) as well as target transcription regulators and thus affect transcription of many genes (Feng et al., 2006; Goodrich and Kugel, 2006). They can also target RNA polymerase II in human and mouse (Espinoza et al., 2007; Mariner et al., 2008) and thus affect the expression of an even broader range of genes. Furthermore, lincRNAs participate in the regulation of splicing (Munroe and Lazar, 1991; Beltran et al., 2008) and translation (Wang et al., 2005; Centonze et al., 2007). Well-characterized examples of lincRNAs involved in epigenetic processes are *Xist* (Brockdorff et al., 1992; Elisaphenko et al., 2008), *Kcnq1ot1* (Umlauf et al., 2004; Pandey et al., 2008), and *Air* (Nagano et al., 2008).

It is well established that, compared to protein-coding sequences and structural RNAs, lincRNAs are weakly conserved in evolution. Many early studies, therefore, branded the lincRNAs "transcriptional dark matter" and considered them to be generally non-functional (Van Bakel and Hughes, 2009; Robinson, 2010). However, low level or lack of detectable conservation does not necessarily imply that these molecules have no function (Pang et al., 2006). A case in point is the best-characterized, functionally important lincRNA gene, *Xist*, which is weakly conserved although it does contain evolutionary constrained regions (Elisaphenko et al., 2008). In general, lincRNAs show reduced substitution and insertion–deletion rates, which has been attributed to purifying selection (Ponjavic et al., 2007; Managadze et al., 2011). Taking into account that some lincRNA genes originated from proteincoding genes [for example, *Xist* (Duret et al., 2006; Elisaphenko et al., 2008)], it appears likely that many properties of lincRNAs would generally resemble those of protein-coding genes, despite the typically lower level of constraint. In particular, protein-coding genes that are highly expressed in many tissues typically evolve slower than genes with lower expression level and breadth (Duret and Mouchiroud, 2000; Krylov et al., 2003; Drummond and Wilke, 2008), and a similar dependence has been observed for lincRNA genes (Managadze et al., 2011). Taken together, these findings imply that an unknown but substantial fraction of lincRNAs are functional molecules rather than transcriptional noise and have evolutionary properties similar to those of protein-coding genes. However, the number of functionally characterized lincRNAs remains scarce (Amaral et al., 2013).

The origin of lincRNA genes generally remains enigmatic. However, analysis of the well-characterized *Xist* lincRNA has revealed fragmentary homology to a protein-coding gene *Lnx3* suggesting that the *Xist* genes emerged in early eutherians via integration of transposable elements (TEs) into the *Lnx3* gene, which gave rise to simple tandem repeats (Duret et al., 2006; Elisaphenko et al., 2008). The *Xist* gene promoter region and 4 of its 10 exons retain homology to exons of the *Lnx3* gene. The remaining six *Xist* exons including those containing simple tandem repeats show similarity to different TEs (Elisaphenko et al., 2008). Integration of TEs into the *Xist* gene apparently had been occurring throughout the course of evolution of this gene and most likely continues in contemporary eutherian species. Additionally, it has been shown that the combination of remnants of protein-coding sequences and TEs is not unique to the *Xist* gene but is also found in neighboring genes that encode non-coding nuclear RNAs (Elisaphenko et al., 2008; Kolesnikov and Elisafenko, 2010).

The discovery of the pivotal contribution of TEs to the evolution of the *Xist* gene prompts the question on a possible general role of TEs in the evolution of lincRNAs. Diverse TEs are widespread and abundant in the genomes of most eukaryotes (Smit, 1996; Brosius, 1999; Kidwell and Lisch, 2001; Deininger and Batzer, 2002). Different classes of TEs include mobile retrovirus-like elements, or retrotransposons, which transpose within the genome via RNA intermediates, and DNA transposons, which can relocate directly. Retrotransposons including long interspersed repetitive elements (LINEs), short interspersed repetitive elements (SINEs), and long terminal repeat (LTR) retrotransposons are widely represented in mammals (Smit, 1996; Deininger and Batzer, 2002). The LINEs are transcribed by RNA polymerase II and contain open reading frames (ORFs) (Temin, 1985). A complete and transpositionally active L1 element (the most common variety of LINEs) is ~7 kb long and contains a 5′-untranslated region (UTR) with an internal promoter, two ORFs (ORF1 and ORF2) and a 3′-UTR terminated by a polyadenylate-rich tail (Smit, 1996; Deininger and Batzer, 2002). The ORF1 encodes a putative RNA-binding protein ~40 kDa in size (Martin, 2006) whereas ORF2 encodes a protein with endonuclease and reverse transcriptase (RT) activities that generates cDNAs from RNA transcripts of the element (Loeb et al., 1986). The mobility of the LINE elements had been demonstrated in mouse and human genomes (Kazazian et al., 1988; Boccaccio et al., 1990). The SINEs are characterized by the presence of a split intragenic RNA polymerase III promoter and a 3′A-rich region often followed by an oligo(A) tail (Smit, 1996; Rogozin et al., 2000; Kapitonov and Jurka, 2003). The SINEs do not contain long ORFs and do not encode enzymes for transposition. Instead, transposition of SINEs apparently requires RT encoded by other TEs, in particular, LINEs (Smit, 1996; Deininger and Batzer, 2002). The LTR retrotransposons have LTRs that range from ~100 bp to over 5000 bp in size (Smit, 1996; Deininger and Batzer, 2002). The LTR retrotransposons are similar to retroviruses in organization, with transcriptional regulatory sequences located in the flanking LTRs, a RT priming site that is typically located immediately downstream of an first LTR, and several ORFs encoding proteins involved in retrotransposition, in particular, RT and integrase (Smit, 1996; Deininger and Batzer, 2002).

The TEs are the primary contributors to the bulk of the genomic DNA in many eukaryotes, in particular mammals, and have the potential to contribute to the evolution of the hosts by providing novel regulatory or coding sequences (Makalowski, 2000). Different classes of regulatory regions in the human genome have been surveyed for the presence of TE-derived sequences (TES) to systematically assess the potential contribution of TEs to the regulation of human genes, and almost 25% of the analyzed promoter regions have been found to contain TES (Jordan et al., 2003; Feschotte, 2008; Bourque, 2009). In addition, numerous examples where experimentally characterized *cis*-regulatory elements are derived from TE sequences have been identified (Jordan et al., 2003; Bourque et al., 2008; Faulkner et al., 2009). Thus, thousands of human (and other mammalian) genes appear to be regulated, at least in part, by sequences derived from TEs (Jordan et al., 2003; Feschotte, 2008; Bourque, 2009). The TES are likely to have substantially contributed to evolutionary change in both gene specific and global patterns of mammalian protein-coding gene regulation (Makalowski, 2000; Jordan et al., 2003).

In light of the regulatory and structural effects that some TEs exert on host protein-coding and lncRNA genes (Makalowski, 2000; Jordan et al., 2003; Elisaphenko et al., 2008; Mattick et al., 2010; Wang et al., 2011; Kapusta et al., 2013; Johnson and Guigo, 2014), we sought to examine the contribution and conservation of TES to regulatory regions, exons and introns of human and mouse lincRNA genes. We found that introns of lincRNA genes contain the highest fraction of TES, followed by exons. The promoters of the lincRNAs contain the lowest fraction of TES but the largest fraction of ancient TES that are conserved between primates and rodents. The content of TES in lincRNA genes is substantially greater than in protein-coding genes, particularly in exons and promoter regions. These results are compatible with the view that TEs are major contributors to the origin and evolution of lincRNAs. We further sought to assess the potential utility of TES as an "evolutionary variable" by analyzing the correlations between the TES content, lincRNA expression, and sequence conservation.

### Materials and Methods

Human and mouse lincRNA genes, the corresponding genomic alignments and expression data were taken from our previous work (Managadze et al., 2011) where the procedures of data processing are described in full details. Briefly, complete mouse and human probe sets were downloaded from the NRED database (Dinger et al., 2009) in the tab delimited and browser extensible data (BED, containing genomic coordinates) formats. The probe sets from platform GNF Atlas 2 (Mouse and Human), with the target classification "Non-coding Only," were used for further analysis. This protocol yielded 917 human and 5444 mouse probe sets. Only the probe sets that mapped to intergenic regions of the human and mouse genomes (i.e., between two adjacent protein-coding genes) were used for further analysis. The resulting list of lincRNAs was further filtered: sequences shorter than 200 nt were removed. This procedure yielded the final set of NCBI GenBank Accession IDs of 2390 mouse and 589 human lincRNAs and their corresponding microarray expression probe sets. The genomic coordinates and sequences of exons and introns of lincRNA and protein-coding genes were downloaded from the UCSC Table Browser (Karolchik et al., 2004), specifically, from "all\_mrna" tables of mouse mm8 and human hg18 assemblies. Multiple alignments of these regions were fetched from Galaxy (Goecks et al., 2010). For the detection of TES, lincRNA and protein-coding genes were analyzed using RepeatMasker version open-3.1.31 with the following parameters: -w -s -no\_is -cutoff 255 -frag 20000 -gff -species mouse/human. A TE insert was considered ancient if the pairwise alignment between human and mouse orthologous TE sequences was longer than 100 bp and contained <5% insertions/deletions (stringent definition) or 25% insertions/deletions (relaxed definition). Microarray data for normal (non-cancerous) tissues (73 human and 61 mouse tissues) were used to analyze the lincRNA expression. Log2-normalized median values of expression for each probe set across the tissues were calculated (Managadze et al., 2011). As an alternative method of measuring expression levels, the mouse RNA-seq data for eight tissues (the ENCODE project; modENCODE Consortium) were downloaded from the UCSC genome browser Web site2 and pooled together. The RPKM value was calculated for each mouse lincRNA (Managadze et al., 2011). Pairwise evolutionary distances for human–macaque and mouse– rat lincRNA alignments were calculated using the DNADIST program from the PHYLIP package (Felsenstein, 1996), with the Kimura nucleotide substitution model. The lists of lincRNA genes and expression data are available at ftp://ftp.ncbi.nlm.nih.gov/pub/ managdav/paper\_suppl/TEs\_lincRNA/.

One of the problems in the analysis of lincRNAs is that there is little overlap between lincRNA sets produced in different studies (Ulitsky et al., 2011; Chew et al., 2013; Managadze et al., 2013; Schuler et al., 2014). We used human and mouse datasets because these curated lincRNA sets have known evolutionary and gene expression properties (Managadze et al., 2011, 2013). Another reason for this choice is that we sought to analyze lincRNA datasets as different as possible from those used in previous studies (Kapusta et al., 2013; Johnson and Guigo, 2014) and to check how small sample size of human lincRNA set influences results. As shown below, the sample size does not perceptibly affect the conclusions of this study.

### Results

### Transposable Elements in Human and Mouse lincRNA Genes

Transposable element-derived sequences (TES for short) comprise at least half of the mammalian genomes, and in particular, are found in most lincRNAs. We identified TES in 69% of the human lincRNAs and 51% of the mouse lincRNAs. These values are somewhat lower than the previously reported 83% of TES in human lincRNAs (Kelley and Rinn, 2012) but nevertheless clearly show the importance of TE for lincRNA evolution. The

1http://www.repeatmasker.org

2http://hgdownload.cse.ucsc.edu/goldenPath/mm9/encodeDCC/ wgEncodeLicrRnaSeq/

distribution of TES in 5′ flanking regions (putative core promoter regions), lincRNA exons, and introns is shown in the **Figure 1**. The lowest fraction of TES was found in the predicted core promoter regions (100 bp upstream regions), and the highest fraction of TES was observed in introns, whereas exonic sequences showed intermediate densities of TES (**Figure 1A**). This distribution of TES is compatible with the previously described general tendency of TES avoidance in functionally important regions of proteincoding genes (Jordan et al., 2003). In particular, similar to the protein-coding genes, the TES density in extended promoter regions has been found to be significantly greater than that in core promoter regions (Jordan et al., 2003). Notably, the fractions of TES in introns of lincRNAs and protein-coding genes are nearly identical, suggestive of comparable (weak) functional constraints. By contrast, in the exons and the core promoter regions of lincRNA genes, the fractions of TES are substantially and statistically significantly (*P* < 10−5 according to the Fisher exact test) higher than in the respective regions of protein-coding genes (compare **Figures 1A,B**). These findings are consistent with the results of a previous study that employed different datasets of lincRNA genes (Kapusta et al., 2013), indicating that the distribution of TES in lincRNA genes is a robust feature. A more detailed analysis of the distribution of TES across lincRNAs is shown in Figure S1 in Supplementary Material. The most prominent

FIGURE 1 | Fractions (proportions) of lincRNA gene regions (concatenated promoters, exons, and introns) (A) and protein-coding gene regions (B) occupied by TE-derived sequences. The differences for pairwise comparisons "promoters vs. introns" and "exons vs. introns" are statistically significant for both classes of genes (*P* < 10−5 according to the Fisher exact test; the raw counts of nucleotides in TES vs. the raw counts of nucleotides in TE-free regions was used as the input for 2 × 2 contingency tables).

feature of this distribution is the high fraction of lincRNAs with a low TES content: in 66% of human lincRNAs and 78% of mouse lincRNAs, the fraction of TES is <20% (Figure S1 in Supplementary Material).

The avoidance of TES in lincRNAs is consistent with purifying selection, which is an important feature of lincRNA evolution (Ponjavic et al., 2007; Managadze et al., 2011). The significant positive correlation between the evolutionary rate and the content of TES was observed for both human and mouse lincRNA sets (for human and mouse, respectively, the Pearson correlation coefficient are 0.183 and 0.337, *P* < 10−5 for both comparisons) (**Figure 2**). We also tested the correlation between the expression level and the content of TES (Figure S2 in Supplementary Material). In many independent previous studies, it has been shown that proteincoding genes that are highly expressed in many tissues typically evolve slower than genes with lower expression level and breadth (Duret and Mouchiroud, 2000; Krylov et al., 2003; Drummond and Wilke, 2008), and a similar dependence has been observed for lincRNA genes (Managadze et al., 2011). Consistent with these observations, here we found a significant negative correlation between the content of TES and the expression level of lincRNAs (Figure S2 in Supplementary Material; for mouse RNA-seq data, Pearson correlation coefficient is −0.158, *P* < 10−5 ; for mouse microarray data, Pearson correlation coefficient is −0.07, *P*< 10−5 ; for human microarray data, Pearson correlation coefficient is −0.253, *P* < 10−5 ).

### Different Classes of Transposable Elements in lincRNA Genes

Analysis of different classes of TEs indicates that the fractions of each class are similar for introns of lincRNA and protein-coding

FIGURE 2 | Correlation between the fraction (proportion) of TEs in concatenated exons and evolutionary rate for human (A) and mouse (B) lincRNAs. Pearson correlation coefficient is 0.183 for human and 0.337 for mouse (*P* < 10−5 for both comparisons).

genes and whole genomes (**Figures 3** and **4**). In each case, the fraction of LINEs is substantially greater than those of SINEs and LTR elements (**Figures 3A** and **4A**). However, there is a significant suppression of LINEs in exonic and promoter regions, in both human and mouse (**Figures 3A** and **4A**). This effect cannot be explained by fluctuations of the base composition in different gene regions because there are no significant compositional differences between exons, introns, and promoter regions for human and mouse lincRNA genes (results not shown). The same trend was observed for different lincRNA sets (Kapusta et al., 2013) suggesting that re-distribution of TEs is a general property of mammalian lincRNA genes. Furthermore, similar tendency is observed in promoter sequences of protein-coding genes (**Figures 3B** and **4B**), the overall lower abundance of TEs notwithstanding. Conceivably, when the smaller SINEs are inserted into functionally important parts of genes, they typically exert a milder deleterious effect than the larger LINEs and LTR elements and accordingly, are more often fixed in the course of evolution.

### Higher Frequency of Ancient Transposable Element-Derived Sequences in Promoters and Exons Compared to Introns

Evolutionary conservation of TEs is likely to reflect molecular domestication of the respective elements (Jordan et al., 2003; Feschotte, 2008; Jurka, 2008; Bourque, 2009; Sinzelle et al., 2009). We analyzed the fraction of ancient mobile elements in different regions of lincRNA genes (**Figure 5**). A significantly higher abundance of ancient TEs (*P* < 10−5 according to the Fisher exact test) was detected in exons and especially in promoter regions compared to introns (**Figure 5**). This finding is consistent with the hypothesis that TEs, in some cases, may perform novel functions in the host organisms (Makalowski, 2000; Jordan et al., 2003). The excess of ancient TEs was more pronounced in human compared to mouse lincRNA genes (**Figure 5**), possibly reflecting differences in evolutionary processes in rodents and primates although a bias caused by technical problems with the detection of 5′-ends of human lincRNA sequences cannot be ruled out (Kutter et al., 2012). We searched the putative promoter regions of lincRNA genes for the presence of TATA boxes and found a substantially elevated frequency of TATA-like sequences in the region −25 to −35 (Figure S3 in Supplementary Material). Given that a similar distribution is observed in many well-annotated human protein-coding genes (Yang et al., 2007), these observations suggest acceptable accuracy of 5′-end identification in lincRNA genes. The fractions of TATA-containing promoters are similar for protein-coding genes (10–25%) (Yang et al., 2007; Anish et al., 2009) and the analyzed sets of lincRNA genes (19–30%; Table S1 in Supplementary Material). The higher frequency of ancient TEs in promoters and exons compared to introns (**Figure 5**; Table S2 in Supplementary Material) suggests that many lincRNA genes emerged before the split of primates and rodents, and that TEs contributed to the origin of these ancient lincRNAs. This finding is consistent with recent observations that 60–70% of the lincRNAs genes are conserved between human and mouse (Kutter et al., 2012; Managadze et al., 2013), and with the observed

sequence [the data for the whole genome sequence were from Waterston et al. (2002)] (A) and protein-coding gene regions (B) occupied by sequences derived from different types of TEs. Differences Fisher exact test; the raw counts of nucleotides in TES vs. the raw counts of nucleotides in TE-free regions was used as the input for 2 × 2 contingency tables).

FIGURE 5 | Ancient transposable elements (TEs) in putative promoter regions, exons, and introns of lincRNA genes. A TE was considered ancient if the alignment between human–mouse orthologous TE sequences was longer than 100 bp and contained <5% insertions/deletions (the stringent threshold). "% TEs" stands for the fraction (proportion) of ancient TEs. Results for the relaxed threshold (the alignment between human–mouse orthologous TE sequences was longer than 100 bp and contained no more than 25% insertions/deletions) are shown in the Table S2 in Supplementary Material. Differences between pairwise comparisons "promoters vs. introns" and "exons vs. introns" are statistically significant (*P* < 10−5 according to the Fisher exact test; the raw counts of ancient TES vs. the raw counts of lineage-specific TES was used as the input for 2 × 2 contingency tables).

higher conservation of lincRNA promoter regions compared to exons (Elisaphenko et al., 2008; Kapusta et al., 2013; Johnson and Guigo, 2014).

## Discussion

The staggering evolutionary success of TEs in eukaryotes is often attributed to their ability to out replicate the host genomes in which they reside, as opposed to any selective advantage that they might provide to their hosts. Indeed, it has been shown that TEs can spread within and among genomes even in the face of a selective cost to the host (Hickey, 1982). Hence, the selfish DNA concept of TEs focuses on the parasitic nature of these elements and emphasizes the deleterious effects of transposition as well as the negligible evolutionary benefit that TEs provide to their hosts (Orgel et al., 1980; Gould and Vrba, 1982). However, the sheer abundance of TEs in the genome, as well as the variety of mutation effects induced by their mobility, suggest that they might, in some cases, be exapted (Gould and Vrba, 1982) or domesticated (Miller et al., 1999), to serve the evolutionary interests of the host (Makalowski, 2000; Jordan et al., 2003). Indeed, multiple lines of evidence indicate that the presence of TEs can result in host adaptation by shaping and reshaping the genome in many different ways (Smit, 1996; Makalowski, 2000; Rogozin et al., 2000; Kidwell and Lisch, 2001; Deininger and Batzer, 2002; Jordan et al., 2003).

The TEs comprise at least half of the mammalian genomes, and in particular, are found in most lincRNAs [this study and Kelley and Rinn (2012)]. Here, we demonstrate that TEs substantially contribute to the evolution of lincRNAs and their promoter regions. Although the densities of TES in these regions are much lower than those in introns, ostensibly, due to the purifying selection that affects functional regions, the contributions of TEs to the evolution of these regions is substantially greater than in the respective regions of protein-coding genes. The higher density of TES in the exons of lincRNAs compared to protein-coding exons appears to reflect the much lower level of functional constraint characteristic of the former (Ponjavic et al., 2007; Managadze et al., 2011). The promoters of lincRNA appear to similarly enjoy greater plasticity and flexibility compared to the promoters of protein-coding genes.

Thus, TE insertion is an important factor that affects lincRNA evolution and biological function. An analysis of TEs in human lincRNAs revealed that the TES composition in lincRNA genes significantly differs from genomic averages: LINEs and SINEs are depleted whereas LTR retrotransposons are enriched (Kelley and Rinn, 2012). The TES occur in biased positions and orientations at lincRNA transcription start sites suggesting a functional role in lincRNA transcriptional regulation (Kelley and Rinn, 2012). In many cases, lincRNAs devoid of TES are expressed at higher levels than lincRNAs containing TES in all tested tissues and cell lines (Kelley and Rinn, 2012). Thus, it has been suggested that TES divide lincRNAs into classes and have contributed to lincRNA evolution and function by conferring tissue-specific expression from extant transcriptional regulatory signals (Kelley and Rinn, 2012). Here, we add another facet to these observations by showing that the promoter regions of lincRNAs are specifically enriched for ancient TES. This finding indicates that not only have many lincRNA genes evolved before the radiation of primates and rodents but also that at least some features of their regulation were already established at that time through TE insertion.

The possibility that some lincRNA genes encode short peptides that are translated, perhaps in a tissue-specific manner, is the subject of an ongoing debate (Brosius and Tiedge, 2004; Mattick and Makunin, 2006; Dinger et al., 2008; Makalowska et al., 2010; Carvunis et al., 2012; Chew et al., 2013). It is extremely hard to rule out such a role for a fraction of purported lincRNAs. A recent peptidomics study demonstrated that most annotated lincRNAs do not generate stable protein (peptide) products (Banfai et al., 2012). Furthermore, ribosomal profiling of lincRNAs suggests that ribosomal engagement with lincRNAs is likely to be regulatory (Chew et al., 2013). The presence of ORFs in the analyzed lincRNA data sets had been analyzed before using different approaches (Managadze et al., 2011, 2013). Importantly, removal of ORFcontaining lincRNAs did not affect the conclusions of both studies (Managadze et al., 2011, 2013). The much higher abundance of TES in lincRNA compared to 5′UTR and protein-coding regions of mRNAs is consistent with the low frequency of protein-coding regions in the analyzed data sets.

It has been proposed that lncRNAs are organized into combinations of discrete functional domains, but the nature of these domains and their identification remain elusive (Guttman and Rinn, 2012). Insertion of TEs and exaptation of TES could represent an important route of evolution of the domain structure of lncRNAs. More specifically, Johnson and Guigo (2014) have proposed that exonic TES comprise functional domains of lncRNAs that they dubbed repeat insertion domains of LncRNAs (RIDLs). A growing number of RIDLs have been experimentally identified whereby lncRNA TES function as RNA-, DNA-, and protein-binding domains/motifs (Elisaphenko et al., 2008; Kelley and Rinn, 2012; Grote and Herrmann, 2013; Holdt et al., 2013; Johnson and Guigo, 2014). These examples are likely to reflect a more general phenomenon of exaptation and/or domestication during lncRNA evolution whereby TES are employed as DNA-, RNA-, and protein-binding domain/motifs (Johnson and Guigo, 2014). The RIDL hypothesis has the potential to explain how functional evolution can keep pace with the fast evolution observed in many lncRNA genes (Johnson and Guigo, 2014). The findings on the distribution of TES across different regions of lincRNA genes, the higher occurrence of TES in lincRNA promoters and exons compared to introns, and significant correlations between the content of TES and evolutionary rate presented here appear to be compatible with the RIDL hypothesis. More specifically, even if a substantial fraction of TES are not fixed in lincRNA exons and promoter regions, those TES that are fixed tend to persist in the genome longer than intronic TES. Moreover, given the near ubiquity of recognizable TES in lincRNA genes, TE mapping can be a useful approach for characterization of lincRNAs and possibly even prediction of their functions. The correlations between the content of TES and various features of lincRNA genes described here could be useful for the characterization of lincRNA functions.

## Conclusion

The results of the present analysis, along with several previous studies, indicate that TEs have contributed to the evolution of many if not most mammalian lincRNAs. Whereas the density of TES in the introns of lincRNA genes is about the same as in introns of protein-coding genes exons, and promoters of lincRNAs are markedly enriched in TES compared to the counterparts in protein-coding genes. This high prevalence of TES reflects the relatively weak evolutionary constraints on lincRNA genes and itself appears to contribute to the plasticity and functional diversification of lincRNAs. Furthermore, the distribution of TE types in the functional regions of lincRNA genes significantly differs from that in introns (or whole genomes), conceivably, because the smaller SINEs that encode no proteins are more suitable for exaptation than the larger, protein-coding LINEs. The prodigious exaptation of TE could account, at least in part, for the functionality of many lincRNAs despite their rapid evolution.

## Acknowledgments

We thank Joshua Cherry, Jean Thierri-Mieg, Kira Makarova, and Yuri Wolf for useful discussions. *Funding*: This research was supported in part by the Intramural Research Program of the National Library of Medicine at the National Institutes of Health and the National Institutes of Health Intramural Research Program of the National Eye Institute (US Department of Health and Human Services), Italian Ministry of Education, University and Research (MIUR) Flagship "InterOmics" (PB05), HIRMA (RBAP11YS7K) and by the "MIMOMICS" European projects.

## Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/article/10.3389/fbioe.2015.00071/ abstract

### References

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Kannan, Chernikova, Rogozin, Poliakov, Managadze, Koonin and Milanesi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **SPECTRA: an integrated knowledge base for comparing tissue and tumor-specific PPI networks in human**

#### *Giovanni Micale 1 , Alfredo Ferro 2 , Alfredo Pulvirenti 2 \* † and Rosalba Giugno 2 \* †*

*1 Department of Computer Science, University of Pisa, Pisa, Italy, 2 Department of Clinical and Molecular Biomedicine, University of Catania, Catania, Italy*

### *Edited by:*

*Marco Pellegrini, Consiglio Nazionale delle Ricerche, Italy*

#### *Reviewed by:*

*Arsen Arakelyan, Institute of Molecular Biology, Armenia Mohammed El-Kebir, Brown University, USA*

#### *\*Correspondence:*

*Alfredo Pulvirenti and Rosalba Giugno, Department of Clinical and Molecular Biomedicine, University of Catania, Via Andrea Doria 6, Catania 95037, Italy<EMAIL_ADDRESS><EMAIL_ADDRESS>
*†Alfredo Pulvirenti and Rosalba Giugno have contributed equally to this work.*

#### *Specialty section:*

*This article was submitted to Bioinformatics and Computational Biology, a section of the journal Frontiers in Bioengineering and Biotechnology*

> *Received: 31 October 2014 Accepted: 17 April 2015 Published: 08 May 2015*

#### *Citation:*

*Micale G, Ferro A, Pulvirenti A and Giugno R (2015) SPECTRA: an integrated knowledge base for comparing tissue and tumor-specific PPI networks in human. Front. Bioeng. Biotechnol. 3:58. doi: 10.3389/fbioe.2015.00058* Protein–protein interaction (PPI) networks available in public repositories usually represent relationships between proteins within the cell. They ignore the specific set of tissues or tumors where the interactions take place. Indeed, proteins can form tissue-selective complexes, while they remain inactive in other tissues. For these reasons, a great attention has been recently paid to tissue-specific PPI networks, in which nodes are proteins of the global PPI network whose corresponding genes are preferentially expressed in specific tissues. In this paper, we present SPECTRA, a knowledge base to build and compare tissue or tumor-specific PPI networks. SPECTRA integrates gene expression and protein interaction data from the most authoritative online repositories. We also provide tools for visualizing and comparing such networks, in order to identify the expression and interaction changes of proteins across tissues, or between the normal and pathological states of the same tissue. SPECTRA is available as a web server at http://alpha.dmi.unict.it/spectra.

**Keywords: tissue, tumor, database, proteins, interactions**

## **1. Introduction**

In the last 10 years, there has been a rapid growth of available protein–protein interaction (PPI) data. They represent all known physical interactions between proteins within a cell. The collection of PPI data yields a network depicting a global overview of the relationships between proteins.

Nowadays, PPIs of species and associated data are stored in many databases, which usually are weekly or monthly updated. Primary sources of PPI data include BioGRID (Stark et al., 2006), DIP (Xenarios et al., 2000), HPRD (Peri et al., 2004), IntAct teorchard2013mintact, and MINT (Licata et al., 2012).

DIP (Xenarios et al., 2000) was the first database, which combined information from multiple observations and experimental techniques into networks of interacting proteins for different species. HPRD (Peri et al., 2004) contains manually curated proteomic information regarding human proteins, which are annotated and linked to OMIM database (Hamosh et al., 2005). BioGRID (Stark et al., 2006) collects protein–protein and genetic interactions for all major model organisms trying to remove redundancy and create a single mapping of interactions. The IntAct database (Orchard et al., 2013) provides tools for both textual and graphical representations of protein interactions. Interacting proteins can be annotated with GO terms for functional analysis. MINT (Licata et al., 2012), which is based on the IntAct database infrastructure, collects experimentally verified PPIs by extracting experimental evidences from the scientific literature.

Some databases integrate PPIs data of human and other organisms from primary sources, by removing redundancies and assigning a unique reliability score. These include STRING (Franceschini et al., 2013), IRefIndex (Razick et al., 2008), ConsensusPathDB (Kamburov et al., 2013), and HitPredict (Patil et al., 2011).

STRING (Franceschini et al., 2013) combines physical interaction data and curated pathways of different organisms with predicted interactions from text mining, genomic features and interactions transferred from model organisms based on orthology. IRefIndex (Razick et al., 2008) is a set of tools to index and retrieve proteins and interactions from major public databases. Indexes are built according to protein sequences and taxonomy identifiers and mapping scores evaluate the quality of the mapping. ConsensusPathDB (Kamburov et al., 2013) integrates human protein–protein interactions, biochemical pathways, gene regulatory, and drug–target interactions into a global network, containing genes, proteins, and metabolites, which can be visualized, analyzed, and annotated. HitPredict (Patil et al., 2011) combines PPI data from IntAct (Orchard et al., 2013), BIOGRID (Stark et al., 2006), and HPRD (Peri et al., 2004), by assigning a confidence score based on sequence, structure, and functional annotations of the interacting proteins. The reliability score is calculated using the Bayesian networks.

The analysis of PPI networks has provided novel biological insights on the function of many previously uncharacterized proteins in *Human* through module identification (Bader and Hogue, 2003; Adamcsek et al., 2006; Mete et al., 2008; Rhrissorrakrai and Gunsalus, 2011), network querying (Ferro et al., 2007; Banks et al., 2008; Bruckner et al., 2010), and network alignment (Flannick et al., 2006; Kalaev et al., 2009; Liao et al., 2009; Sahraeian and Yoon, 2013; Micale et al., 2014a) algorithms. Furthermore, the annotation of PPI networks with external data (i.e., diseases, expression data, phenotypes) has helped to classify genes according to the expression profiles (Dao et al., 2011), predict new gene–disease associations (Huang et al., 2012; Zhao et al., 2012), and discover new drugs (Huang et al., 2012; Alaimo et al., 2013; Csermerly et al., 2013).

These tasks have been accomplished thanks to the availability of authoritative repositories of gene expression data in normal/cancer tissues and at different diseases stages (Uhlen et al., 2010; Barrett et al., 2013; Rustici et al., 2013). For example, ArrayExpress (Rustici et al., 2013) and GEO (Barrett et al., 2013) include gene expression data from microarray and highthroughput sequencing experiments, which can be easily queried or downloaded. Users can also submit data directly by using the standard MIAME format. More recently, new projects have started with the aim of cataloging tissue or tumor sequencing data. The Cancer Genome Atlas (TCGA)1 collects complete highthroughput genome data (clinical information, expressions data, methylations, mutations) for specific cancer tissues, with the purpose of helping the diagnosis and the treatment of cancers. The Human Protein Atlas (Uhlen et al., 2010) is a database with histological images showing the spatial distribution of proteins in normal and cancer tissues. Protein Atlas contains also transcription expression levels, protein expression profiles, and subcellular localization data.

All above PPI networks data are constructed by ignoring the role of proteins in human tissues. On the other hand, human diseases often occur in specific tissues (Lage et al., 2008). Some genes can be predominantly expressed in one or few tissues and can control the formation of protein complexes (Emig and Albrecht, 2011). Furthermore, genes can use alternative splicing as a powerful mechanism to enlarge the number of their interactors and perform distinct functions in different tissues (Emig and Albrecht, 2011). Therefore, the integration of PPI networks with tissuespecific gene expression data can help to highlight the role of some genes in specific disease or tumors. The result of such integration gives the so called Tissue-Specific PPI (TS-PPI) network (Bossi and Lehner, 2009), which is a subgraph of a PPI network where the genes corresponding to both interacting proteins are expressed in one or more selected tissues.

Some studies focus on the analysis of global and local properties of TS-PPI networks. In Bossi and Lehner (2009), authors prove that most housekeeping proteins form highly tissue-specific protein interactions, suggesting a key role of those proteins in tissuespecific biological processes. Emig and Albrecht (2011) show that the number of tissue-specific proteins is very low and the receptoractivated signaling processes and the transcriptional regulation are two key factors for tissue specificity. In Souiai et al. (2011), a gradient model is used to describe the structure of TS-PPI networks, containing interactions of regulatory and developmental functions at the core of the TS-PPI network and physiological functions at the periphery.

Several recent works highlight the advantages of using TS-PPI networks. In Lopes et al. (2011), a set of proteins related to the response of viral infection in a TS-PPI network lead to a more reliable functional enrichment. Magger et al. (2012) use TS-PPI networks to improve the prioritization of candidate diseasecausing genes with respect to a generic PPI network. In Chen and Wang (2012), authors identify functional modules in TS-PPI networks using CFinder (Adamcsek et al., 2006) and show that they exhibit more biological meaning than modules in a PPI network. Xiao et al. (2014) propose a new method for the identification of multi-tissue gene co-expression networks associated with specific functional processes relevant for phenotype variation and disease in humans. Barshir et al.(2014)show that genes causing hereditary diseases tend to have higher transcript levels and more interacting partners in the TS-PPI network of disease tissues than in the TS-PPI network of unaffected tissues.

To the best of our knowledge, few tools are available for querying and analyzing TS-PPI networks (Barshir et al., 2013; Nersisyan et al., 2014). CyKeggParser (Nersisyan et al., 2014) is a Cytoscape app for generating and analyzing tissue-specific KEGG pathways. Pathways can be checked for inconsistencies and modified based on gene expression data from normal and cancer tissues. TissueNet (Barshir et al., 2013) is a dataset of TS-PPIs in humans, which integrates a collection of four PPInetworks (BioGRID, DIP, IntAct, and MINT) with three expression datasets (GEO, Human Protein Atlas, and Illumina Body Map 2.0). The database provides

1 http://cancergenome.nih.gov

Micale et al. SPECTRA database

a web interface for retrieving tissue-specific interactions of a query protein. However, it handles only 16 normal tissues and does not provide any tool for the analyses of TS-PPI networks.

In this paper, we propose SPECTRA (SPECific Tissue/Tumor Related PPI networks Analyzer), a framework to build and analyze TS-PPI networks. SPECTRA integrates tissue and tumorspecific gene expression data from the most authoritative online repositories such as Protein Atlas, ArrayExpress, GEO, and TCGA. Expression data are then integrated with high-quality protein–protein interactions, taken from HPRD, BioGRID, MIPS, IntAct, and the work of Havugimana et al. (2012). We provide a web interface for constructing, visualizing, and comparing TS-PPI networks, with the aim of identifying differential interaction/expression patterns in TS-PPI networks (i.e., distinct tissues, or normal and pathological states of the same tissue). The TS-PPI networks together with the results of differential analysis can be easily visualized by using Cytoscape facilities (Shannon et al., 2003) and downloaded as text files for further investigations. SPECTRA is free for all users and available at http://alpha.dmi. unict.it/spectra.

### **2. Materials and Methods**

SPECTRA combines protein–protein interactions in human with gene expressions, by integrating 13 authoritative resources. The final integrated SPECTRA database contains 16,435 protein coding genes and 175,841 gene interactions (GIs), 1,350,637 tissuespecific gene expression data entries covering 107 normal tissues, and 2,171,808 tumor-specific expression data entries covering 160 different tumors.

### **2.1. Interaction Datasets**

Human protein interaction data were taken from BioGRID2 , DIP3 , a recent work by Havugimana et al. (2012), HPRD4 , IntAct5 , and MINT6 .

**Table 1** describes the features of the PPI networks integrated in SPECTRA. Networks taken from the work of Havugimana et al. (2012), IntAct and MINT are weighted with edge weights ranging in [0,1], while the other PPI networks are unweighted. Proteins of the considered PPI networks, including splicing isoforms, were first mapped to the corresponding gene. Next, a global GI network was built, by collecting all interactions reported

2 http://thebiogrid.org

5 http://www.ebi.ac.uk/intact

6 http://mint.bio.uniroma2.it/mint

in at least one dataset. We assigned to each edge a pair consisting of the average value of weights across the datasets that report that interaction and the percentage of datasets giving the interaction (dataset coverage). Average edge weights range from 0.131 to 1.

**Figure 1** depicts a Venn diagram of common gene interactions between PPI datasets. Interaction databases generally show low overlap, with only 25 interactions shared by all datasets and only 7,783 interactions in common between MINT, BioGRID, IntAct, and HPRD, which are the biggest ones. The final integrated network has 16,435 nodes, 175,841 edges and 17 connected components, with a high average diameter (9) and low clustering coefficient (0.289). The average degree is 21.398 and the degree distribution follows a power law (**Figure 2**).

3 http://dip.doe-mbi.ucla.edu/dip

4 http://www.hprd.org

**TABLE 2 | Features of expression datasets integrated in SPECTRA**.

### **2.2. Expression Datasets**

Gene expression data for various tissues and tumors were downloaded from ArrayExpress7 , GEO8 , ProteinAtlas9 , and TCGA (see text footnote 1). **Table 2** lists the gene expression datasets integrated within SPECTRA, the platform used to detect the expressions and the number of covered tissues and tumors.

**Figure 3** depicts a Venn diagram of common tissues and tumors across expression datasets. While tissue names are generally shared, tumor names are much differentiated, resulting in a poor overlap between datasets. In particular, TCGA contains data for very specific tumors and partially overlap only with E-MTAB-62 dataset, which is the richest one. Note that the numbers reported in **Figure 3** only refer to specific tumors and not to tumor classes. So, for instance, "breast carcinoma" and "breast adenocarcinoma" are considered distinct tumors, even though they belong to the same class of tumors, "breast cancer."

As regards the integration of expression data, we followed the work of Guo et al. (2013), where authors show that there is a positive correlation between normalized Affymetrix and RNA-Seq data. We performed RMA normalization (McCall et al., 2010) for datasets based on Affymetrix platforms (GDS181, GDS596, and GDS1096), using the corresponding R Bioconductor package. For GDS3113, ProteinAtlas, and TCGA, we first computed the log2 of the number of fragments and then we normalized values using the quantile normalization method of Bolstad et al. (2003). For GDS181, GDS596, and GDS1096 datasets, normalized values were computed from raw data. Then, probes, which were present in a particular microarray dataset, were mapped to the corresponding genes. Finally, the expression of a gene for a specific tissue was computed as the average expression value of probes mapping to that gene in the tissue. For GDS3113, ProteinAtlas, and TCGA, instead, we directly normalized gene expression values for different tissues. We did not normalize EMTAB-62 data, since its source values were already normalized with RMA.

Finally, we assigned to each pair gene–tissue a unique positive expression score, given by the average normalized expression value of the gene in that tissue, according to the different datasets. Expression scores in SPECTRA range from 3.566 to 17.366 for tissues and from 0.01 to 17.343 for tumors.

### **2.3. Data Schema**

SPECTRA database is structured as a MySQL relational database with six tables: *Genes*, *Tissues*, *Tumors*, *Interactions*, *Expr\_normal*, and *Expr\_tumor*.

The *Genes* table contains the list of all expressed and interacting genes. Each entry is identified by the gene symbol and contains associated data, including a description string, aliases, and cross references to Entrez Gene (if available).

The *Tissues* and *Tumors* tables have the same structure. Tissues and tumors are associated to different classes, depending on the organism part they refer to. Each entry is identified by a unique number and contains a description and the corresponding class.

7 http://www.ebi.ac.uk/arrayexpress

8 http://www.ncbi.nlm.nih.gov/geo

9 http://www.proteinatlas.org

SPECTRA contains 26 distinct classes of tissues and 32 distinct classes of tumors.

The *Interactions* table lists all the PPIs integrated in SPEC-TRA. Interactions are identified by a couple of gene symbols and the edge weight for each integrated dataset (when available) is stored, together with the average interaction weight across dataset reporting that interaction and the dataset coverage.

*Expr\_normal* and *Expr\_tumor* contain all the gene expressions in normal and cancer tissues. The unique identifier of *Expr\_normal* is a couple gene–tissue, while entries in *Expr\_tumor* are uniquely identified by the couple gene–tumor. In both tables, the normalized expression value for each integrated dataset (where available) and the average expression score are included as associated data.

### **2.4. An Algorithm for Differential Local Alignment of TS-PPI Networks**

TS-PPI networks are compared in SPECTRA for identifying patterns of differential gene expressions between multiple TS-PPI networks.

Our goal is to find conserved sub-regions in the TS-PPI networks, which maximize the difference of expression values of aligned genes. The problem is related to that of finding maximalscoring connected subgraphs, which is NP-hard, even in a common simpler setting where the aligning TS-PPI networks have the same set of nodes and edges (e.g., TS-PPI networks built starting from different expression data and the same interaction datasets) (Ideker et al., 2002).

In the case of two TS-PPI networks with the same set of nodes and edges (representing for instance case and control expression data), heuristic (Ideker et al., 2002; Sohler et al., 2004; Cabusora et al., 2005; Rajagopalan and Agarwal, 2005; Guo et al., 2007) and exact (Dittrich et al., 2008) solutions have been proposed. However, as far as we are concerned, no solutions are known for the multiple case. Here, we propose an approximate solution to the multiple differential alignment problem based on a modified version of the GASOLINE algorithm (Micale et al., 2014a). For simplicity, we consider TS-PPI networks with no multiple edges between two nodes.

### 2.4.1. The GASOLINE Algorithm

GASOLINE (Micale et al., 2014a) is a greedy and stochastic algorithm for multiple local alignment of protein–protein interaction networks. Flowchart in **Figure 4** provides a general description of GASOLINE.

Given *N* weighted PPI networks of different species, where edge weights are probabilities of interaction between proteins, local alignment aims at finding a set of connected subnetworks, one from each network, that are conserved in their sequence and interaction pattern. Such subnetworks could represent evolutionary conserved complexes or pathways across different organisms.

Such a problem is related to subgraph isomorphism, which is known to be NP-complete (Cook, 1971). GASOLINE proposes an approximate solution through a stochastic-greedy strategy consisting of two phases.

In the first step, called bootstrap phase, we look for orthologous proteins across the networks and build a set of seeds. The set of seeds initially consists of proteins, one from each network, and includes all the starting nodes of the suboptimal local network alignment.

The second step, called iterative phase, repeatedly adds (extension step) or removes (removal step) nodes in the network alignment, trying to maximize the final alignment score. Each extension step adds, in each network, a single node to the corresponding seed. During the extension step, the seeds grow up producing a set of subgraphs, one from each network. The extension process is regulated by a properly defined degree ratio measuring the average density of the aligned subgraphs with respect to their neighbors in the networks. The extension is performed until the degree ratio increases.

Each removal step replaces from the current alignment the set of proteins (one from each network) with minimum topology similarity score.

The bootstrap phase and each extension step are performed through Gibbs sampling (Geman and Geman, 1984). In both cases, the Gibbs sampling builds a chain, where each state represents a combination (i.e., alignment) of single proteins, one from each network. First, a random initial state is selected. Then, the sampling method iteratively performs a transition from a state to another, by replacing a randomly chosen protein of the current alignment with a protein of the same network, according to a properly defined transition probability distribution.

Due to its non-deterministic nature, different iterations of GASOLINE may produce different local alignments. The above steps are iterated to produce a set of local networks alignments, which are then ranked according to an Index of Structural Conservation (ISC) score. ISC score measures the percentage of conserved interactions in the final alignment. The higher is ISC, the better is the alignment.

GASOLINE implements preprocessing and post-processing steps. During preprocessing, the search space for potential seeds is reduced. This is obtained by marking only proteins having orthologs in all aligning networks and with a significant interaction degree in each network. All marked nodes in each network *Gi*(1 *≤ i ≤ N*) are added to a set called *Si* . These sets will be used in the initial phase and will be updated at each iteration. Finally, a post-processing filters the final set of local alignments returned by GASOLINE by removing highly overlapping complexes.

GASOLINE does not allow many-to-many mapping between aligned nodes. However, experimental results show that the algorithm can produce more reliable results than methods implementing many-to-many mapping. Moreover, GASOLINE is clearly faster than the state-of-art algorithms (Micale et al., 2014a).

### 2.4.2. The Adapted GASOLINE

We implemented a customized version of GASOLINE to compare two or more Tissue-Specific PPI (TS-PPI) networks for local differential alignment problem. GASOLINE algorithm was extended to deal with gene expressions as weights to the nodes.

Let*A*and *B* two genes and *Expr(A)* and *Expr(B)*their expression values, with *Expr(A) ≥ Expr(B)*. In order to evaluate the expression difference between *A* and *B*, we compute the *log fold change*, defined as follows:

$$\text{LogFold}(A, B) = \log\_2 \left( \frac{\text{Exp} r(A)}{\text{Exp} r(B)} \right) \tag{1}$$

Given *N* TS-PPI networks and a set of aligned genes *G* = {*G*1, *G*2, *. . .* ,*GN*}, one for each TS-PPI network, MaxLogFold is the maximum value of LogFold function among all pairs of genes in*G*:

$$\text{MaxLogFold}(G) = \max\{\text{LogFold}(G\_i, G\_j) \; \forall \; 1 \le i < j \le n\} \tag{2}$$

We applied the following changes to original GASOLINE algorithm:

*•* We included the LogFold function in the Gibbs sampling procedure of bootstrap and iterative phases, by multiplying it by the topology and homology scores in the computation of node similarities;

*•* The number of iterations of Gibbs sampling both in the bootstrap and in the extension phase is governed by a new parameter, is *α*, which is a probability threshold related to *N*, the number of networks, according to the following formula:

$$k = \max\left\{ k' : \left( \frac{N-1}{N} \right)^{k'} > \alpha \right\} \tag{3}$$

where *P* = ( *N−*1 *N* )*k ′* is the probability that a gene is never selected in *k ′* consecutive iterations of Gibbs sampling. The idea is to stop Gibbs sampling when an alignment does not change for *k* consecutive iterations. The lower is *α*, the higher is *k*, so the more precise and slower will be the sampling procedure:

### **3. Results**

SPECTRA is a framework for retrieving and analyzing protein–protein interaction data specific for a given set of normal or cancer tissues. The underlying graph model in SPECTRA is the Tissues-Specific PPI network (or TS-PPI network), in which the genes of corresponding interacting proteins are both expressed in one or more tissues. The architecture of SPECTRA is composed by (i) the *searching tool*, which allows to build TS-PPIs; (ii) the *comparison tool* to look for shared differential expressions patterns between genes of two or more TS-PPI networks. Results can be graphically visualized by using Cytoscape.js or downloaded as text files.

### **3.1. SPECTRA Search Tool: Building TS-PPI Networks in SPECTRA**

SPECTRA builds TS-PPI networks starting from a user-defined set of genes, tissues, expression data, and interaction data. **Figure 5** depicts the search interface of SPECTRA.

In the "Gene data" section (**Figure 5A**), the user can look for all genes expressed in a set of tissues or restrict the search to a specific list of genes. Genes can be provided with their official names or Aliases (e.g., Ensembl Gene, Entrez Gene, Affy).

In the "Expression data" section (**Figure 5B**), the user limits the search to a set of tissues/tumors and to a set of expression datasets or uploads a text file with custom expression data. Note that the two options are mutually exclusive, that is, all the settings concerning datasets and tissue/tumors will be ignored if the user

the interactions that are present at least in Havugimana and HPRD, involving

coverage, and expression scores.

provides a custom text file. Available tissues and tumors in SPEC-TRA are listed in a table and can be easily included in the input query list with a double click in each entry. When no data are provided, all the tissues and tumors in SPECTRA are considered. Tissues and tumors are also mutually exclusive, meaning that a TS-PPI network built-in SPECTRA cannot contains interactions defined on both normal and tumor tissues. However, two TS-PPI networks defined upon a specific set of tissues and tumors, respectively, can be always compared for differential analysis with the adapted GASOLINE. The user can also select one or more datasets from which the expression have to be reported. When the expression is in other datasets it will be also given. When no dataset is selected, all expression data in SPECTRA are considered.

Finally, a further filter on genes can be applied by indicating a threshold for the minimum normalized value of gene expressions to be considered.

The "Interaction data" section (**Figure 5C**) contains the parameters for filtering interaction data. As above, the user can select one or more datasets where protein interactions have to be reported. If no interaction dataset is selected, all PPIs in SPECTRA are considered. A threshold can be provided to select interaction weights above a given value and a minimum dataset coverage.

When all input parameters have been specified, the user clicks on the "Search" button. At the end of the process, all the TS-PPIs found are listed in a result table (**Figure 6**). For each TS-PPI, we show the interacting genes, the tissues where they are expressed, the expression values of genes in tissues, the average interaction weights and dataset coverages of corresponding proteins. Results are ordered by dataset coverage and average interaction weight. Expression values and interaction weights are depicted with colored progress bar, where colors range from cyan (low values) to red (high values).

By selecting a specific TS-PPI in the result table, additional data about the interaction and the interacting genes are shown (**Figures 6** and **7**). A list of datasets reporting the interaction and the corresponding interaction weight is reported on the right of the result table (**Figure 6**). Below the result table, two panels with details about the interacting genes are shown (**Figure 7**). For each gene, description and aliases are provided, together with the lists of tissues and tumors where the gene is expressed, according to the different expression datasets, ordered by expression score.

### **3.2. SPECTRA Comparison Part: Compare TS-PPI Subnetworks**

TS-PPI networks can be compared in SPECTRA for identifying patterns of differential gene expressions between multiple TS-PPI networks. The goal is to find conserved sub-regions in the TS-PPI networks, which maximize the difference of expression values of aligned genes.

**Figure 8** shows the "Compare" tabbed panel in SPECTRA. Before running the adapted GASOLINE, the user has to upload at least two TS-PPI networks. For each network, the number of nodes and edges are reported. Networks can also be renamed by double clicking on the corresponding cell. Note that uploaded TS-PPI networks with multi-edges between nodes will be always treated as simple networks, where multi-edges are replaced by a single edge with weight equals to the average weight of multiedges and label given by the concatenation of the multi-edge labels.

Once the networks have been uploaded, the user can click on "Run GASOLINE" button to set the input parameters for the adapted GASOLINE (**Figure 8**).

We briefly describe their meaning (default values are reported in brackets):

table we report, the tissue, the average expression scores of interacting genes, and the total interaction weight. Expression scores, weights, and dataset coverage are represented with a colored progress bar (from cyan to red). By

detailed data about the interaction are shown to the right. For each dataset, the corresponding interaction weight (when available) is reported (for example, 0.71 for IntAct database).

**FIGURE 7 | The panel with detailed information of a gene**. When an interaction is selected from the result table (**Figure 6**), two panels with additional data, one for each interacting gene, are shown. This example refers to the detailed panel for gene NCBP2, which appears when the row table of **Figure 6** is selected. In the detailed panel, the gene symbol, the