$$1 + 0.00528y^6 + 0.00157y^7 + 0.00042y^8 + 0.00008y^9 + 0.00002y^{10}.\tag{14}$$

From this result, we have the whole distribution of *N*: support is {0, 1, . . . , 10}, **P**(*N* = 0) = 0.33369, **P**(*N* = 1) = 0.31148, ..., **P**(*N* = 10) = 0.00002. We can also easily derive moments of *N* from this distribution: **E**[*N*] = 1.28, *σ*[*N*] = 1.29.

**Lemma 5** (derivatives of the pgf)**.** For any *k* 0, the order *k* derivative of the pgf *G* is given by:

$$G^{(k)}(y) = k! [z^k] \mathbf{u}(\mathbf{P} + y\mathbf{Q} + z\mathbf{Q})^{\ell - d} \mathbf{v} \tag{15}$$

where the [*zk*] operator denotes the extraction of the coefficient of *zk* in the expression.

*Proof.* The formal proof can be found in Nuel (2010) in a slightly less general case. Here we prove it only for the first two derivatives in the particular case where − *d* = 3. Starting from *G*(*y*) = **u**(**P** + *y***Q**)3**v** we get:

$$G'(y) = \mathbf{u} \left( \mathbf{Q} (\mathbf{P} + y\mathbf{Q})^2 + (\mathbf{P} + y\mathbf{Q})\mathbf{Q} (\mathbf{P} + y\mathbf{Q}) + (\mathbf{P} + y\mathbf{Q})^2 \mathbf{Q} \right) \mathbf{v} \tag{16}$$

and

$$\mathbf{G}'''(y) = 2\mathbf{u}\left(\mathbf{Q}^2(\mathbf{P} + y\mathbf{Q}) + \mathbf{Q}(\mathbf{P} + y\mathbf{Q})\mathbf{Q} + (\mathbf{P} + y\mathbf{Q})\mathbf{Q}^2\right)\mathbf{v} \tag{17}$$

which are easily connected to the terms coefficients of *z*1 and *z*2 in **u**(**P** + *y***Q** + *z***Q**)−*d***v**.

If we denote for all *k* 0 the *k*-th *factorial moment* of *N* by *Fk* = **E**[*N*!/(*N* − *k*)!], then, by the definition of the pgf, it is clear that *Fk* = *G*(*k*)(0), and thanks to Lemma 5 we get:

$$F\_k = k! [z^k] \mathbf{u}(\mathbf{T} + z\mathbf{Q})^{\ell - d} \mathbf{v}. \tag{18}$$

And if we now denote the moment generating function (mgf) of *N* by *M*(*t*) = **E**[*etN*] = *G*(*et* ), and the cumulant generating function (cgf) of *N* by Λ(*t*) = log **E**[*etN*] = log *M*(*t*) = log *G*(*et* ), we get directly the *k*-th moment of *N*: **E**[*Nk*] = *M*(*k*)(0); and the *k*-th cumulant of *N*: *κk* = Λ(*k*)(0).

**Corollary 6** (characteristics moments)**.** If we denote by *μ* = *κ*1 the *expectation* of *N*, by *σ* = √*κ*2 the *standard deviation* of *N*, by *γ*1 = *κ*3/*σ*3 the *skewness* of *N*, and by *γ*1 = *κ*4/*σ*4 the *excess kurtosis* of *N*, then we get: *μ* = *F*1, *σ*2 = *F*2 + *F*1 − *F*2 1 ,

$$\gamma\_1 = \frac{3F\_2 - 3F\_1^2 + F\_3 - 3F\_1 F\_2 + 2F\_1^3 + F\_1}{\sigma^3},\tag{19}$$

and

12 Will-be-set-by-IN-TECH

*Proof.* The first equality derives directly from the definition of *G*(*y*). For the second equality now, it is clear that **u**(**P** + **Q**)−*d* gives the marginal distribution of *Z*. We then connect this distribution to *N* by counting the number of times we use the transitions of **Q** with the dummy variable *y* so that **u**(**P** + *y***Q**)−*d* gives the joint distribution of (*Z*, *N*). Finally, we sum up the

For example, let us consider M = G(G|C)G and *X*1:12 generated by a M0 model with

0.6 0.4 0 0 0 0 0.2 0 0.4 0.4 0 0 0.6 0 0 0 0.4*y* 0 0.2 0 0.4 0 0 0.4*y* 0.2 0 0.4 0.4 0 0 0.2 0 0.4 0 0 0.4*y*

= 0.33369 + 0.31148*y* + 0.19357*y*2 + 0.09681*y*3 + 0.04140*y*4 + 0.01569*y*5

From this result, we have the whole distribution of *N*: support is {0, 1, . . . , 10}, **P**(*N* = 0) = 0.33369, **P**(*N* = 1) = 0.31148, ..., **P**(*N* = 10) = 0.00002. We can also easily derive moments

**Lemma 5** (derivatives of the pgf)**.** For any *k* 0, the order *k* derivative of the pgf *G* is given

*Proof.* The formal proof can be found in Nuel (2010) in a slightly less general case. Here we prove it only for the first two derivatives in the particular case where − *d* = 3. Starting from

**Q**(**P** + *y***Q**)2 + (**P** + *y***Q**)**Q**(**P** + *y***Q**)+(**P** + *y***Q**)2**Q**

**Q**2(**P** + *y***Q**) + **Q**(**P** + *y***Q**)**Q** + (**P** + *y***Q**)**Q**2

where the [*zk*] operator denotes the extraction of the coefficient of *zk* in the expression.

which are easily connected to the terms coefficients of *z*1 and *z*2 in **u**(**P** + *y***Q** + *z***Q**)−*d***v**.

definition of the pgf, it is clear that *Fk* = *G*(*k*)(0), and thanks to Lemma 5 we get:

If we denote for all *k* 0 the *k*-th *factorial moment* of *N* by *Fk* = **E**[*N*!/(*N* − *k*)!], then, by the

And if we now denote the moment generating function (mgf) of *N* by *M*(*t*) = **E**[*etN*] = *G*(*et*

and the cumulant generating function (cgf) of *N* by Λ(*t*) = log **E**[*etN*] = log *M*(*t*) = log *G*(*et*

we get directly the *k*-th moment of *N*: **E**[*Nk*] = *M*(*k*)(0); and the *k*-th cumulant of *N*: *κk* =

+0.00528*y*6 + 0.00157*y*7 + 0.00042*y*8 + 0.00008*y*9 + 0.00002*y*10. (14)

⎞

12

×

(*y*) = *k*![*zk*]**u**(**P** + *y***Q** + *z***Q**)−*d***v** (15)

*Fk* = *k*![*zk*]**u**(**T** + *z***Q**)−*d***v**. (18)

⎛

⎞

⎟⎟⎟⎟⎟⎠

�

�

**v** (16)

**v** (17)

),

),

(13)

⎜⎜⎜⎜⎜⎝

⎟⎟⎟⎟⎟⎠

parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40. Proposition 4 hence gives:

contributions of all states using the product with **v**.

100000 �

of *N* from this distribution: **E**[*N*] = 1.28, *σ*[*N*] = 1.29.

*G*(*k*)

×

⎛

⎜⎜⎜⎜⎜⎝

*G*(*y*) = �

*G*(*y*) = **u**(**P** + *y***Q**)3**v** we get:

*G*�

(*y*) = **u**

*G*��(*y*) = 2**u**

�

�

by:

and

Λ(*k*)(0).

$$\gamma\_2 = \frac{7F\_2 - 7F\_1^2 + 6F\_3 - 18F\_1F\_2 + 12F\_1^3 + F\_4 - 4F\_1F\_3 - 3F\_2^2 + 12F\_1^2F\_2 - 6F\_1^4 + F\_1}{\sigma^4}.\tag{20}$$

*Proof.* On just need to compute the derivatives Λ(1)(0), Λ(2)(0), Λ(3)(0), and Λ(4)(0).

If we consider again M = G(G|C)G and *X*1:12 generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40. Eq. (18) hence gives:

$$\sum\_{k<0} \frac{F\_k}{k!} = \begin{pmatrix} 1 \ 0 \ 0 \ 0 \ 0 \ 0 \end{pmatrix} \times \begin{pmatrix} 0.6 \ 0.4 & 0 \ 0 & 0 & 0\\ 0.2 & 0 \ 0.4 \ 0.4 & 0 & 0\\ 0.6 & 0 & 0 & 0.4 + 0.4y & 0\\ 0.2 & 0 & 0.4 & 0 & 0 & 0.4 + 0.4y\\ 0.2 & 0 & 0.4 \ 0.4 & 0 & 0\\ 0.2 & 0 & 0.4 & 0 & 0 & 0.4 + 0.4y \end{pmatrix}^{12} \times \begin{pmatrix} 1\\ 1\\ 1\\ 1\\ 1 \end{pmatrix} \tag{21}$$

$$= 1 + 1.28y + 1.01683y^2 + 0.61211y^3 + 0.29709y^4 + 0.11835y^5$$

$$+ 0.03845y^6 + 0.00992y^7 + 0.00193y^8 + 0.00025y^9 + 0.00002y^{10}. \tag{22}$$

From this result, we can get all factorial moments of *N*: **E**(1) = *F*0 = 1, **E**(*N*) = *F*1 = 1.28, **E**(*N*(*N* − 1)) = *F*2 = 2.033664, **E**(*N*(*N* − 1)(*N* − 2)) = *F*3 = 3.6726374, **E**(*N*(*N* − 1)(*N* − 2)(*N* − 3)) = *F*4 = 7.1302266, ..., **E**(*N*!/(*N* − 10)!) = *F*10 = 60.881161. Thanks to Corollary 6 we get the following characteristic moments: *μ* = 1.28, *σ* = 1.294320, *γ*1 = 1.163783, *γ*2 = 1.492661.

#### **3.5 Exact computations**

As we have seen above, Proposition 4 provides a way to obtain the whole distribution of *N* by computing *G*(*y*) = **u**(**P** + *y***Q**)�−*d***v** from which we can easily derive *S*(*n*) for any *n* 0:

$$S(n) = \begin{cases} +\log\_{10}\left(\sum\_{k=0}^{n} [y^k] G(y)\right) \text{ if } n \leqslant \mathbb{E}[N] \\ -\log\_{10}\left(\sum\_{k=n}^{+\infty} [y^k] G(y)\right) \text{ if } n > \mathbb{E}[N] \end{cases}$$

.

From the algorithmic point of view, there are basically two approaches to compute *S*(*n*) using Expression (12). The first one, called *power*, consists in computing (**P** + *y***Q**)�−*d* using the power method and a binary decomposition of � − *d*. Ex: if � − *d* = 1097 then � − *d* = 210 + 26 + 23 + 20. We then just have to recursively compute **D***k*(*y*)=(**P** + *y***Q**)2*k* using the relation **D***k*+1(*y*) = **D***k*(*y*) × **D***k*(*y*) for all *k* 0. Since in the computation of *S*(*n*) we are only interested in terms of degree *n* or less (or *n* or more), we can easily truncate 2 all polynomials at degree *n* thus dramatically reducing the computational costs of polynomial products. We end

2 In the case of over-representation, all contributions of degree *n* or more are summed into the term of degree *n*.

in Biological Sequences 15

Significance Score of Motifs in Biological Sequences 187

Table 1. Characteristic moments the number *N* of occurrences of motif M = G(G|C)G in a sequence *X*1: generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and

Since the random count *N* is basically defined by Eq. (2) as large sum of Bernouilli variables, the idea of approximating the distribution of *N* using Gaussian approximation sounds appealing. Indeed, Gaussian approximations are historically the first ones to have been suggested for this problem (Cowan, 1991; Kleffe & Borodovski, 1997; Pevzner et al., 1989; Prum et al., 1995). From the theoretical point of view, Central Limit Theorems (CLT) for weakly dependent variables ensure that *N* is asymptotically normal distributed. On Table 1, we can see the characteristic moments of *N* for motif M = G(G|C)G and various value of the sequence lengths . According to theory, we observe that the skewness and excess kurtosis both decease toward 0 when grows (a normal distribution has null skewness and excess kurtosis). But it is also clear that *N* is not normally distributed for small values of . As a consequence, the quality of a Gaussian approximation for *S*(*n*) is expected to be questionable

In order to overcome this issue, Nuel (2010) suggested to consider near Gaussian approximations instead of simple Gaussian approximations for this problem. The idea is simply to perform a higher order asymptotic development that exploits more than the two first moments of *N*. This technique is known as the Edgeworth's expansion. Blinnikov & Moessner (1998) gives a general (and rather complicated) formula for this expansion. For

probability distribution function (pdf) of a standard Gaussian, then for all *n* 0 we have

6 *H*3(*z*) *C*2(*z*) = *S*4

where *μ* = **E**[*N*], *σ* = **V**[*N*], *z* = (*n* − *μ*)/*σ*, *Sk* = *κk*/*σ*2*k*−2 for all *k* 1, and where *Hk*(*z*) are the Hermite polynomials defined recursively by *H*0(*z*) = 1 and *Hk*(*z*) = *zHk*−1(*z*) −

*C*0(*z*) + *σC*1(*z*) + *σ*2*C*2(*z*) + *σ*3*C*3(*z*)

144 *H*7(*z*) + *S*3

24 *H*4(*z*) + *S*2

3

3

√2*π* the

(24)

1296 *H*9(*z*) (26)

72 *H*6(*z*) (25)

*π*(C) = *π*(G) = 0.40. Computation performed using the power approach.

practical purpose, we present the result only up to order 3 expansions.

*σ*

**P**(*N* = *n*) � *ϕ*(*z*)

*C*0(*z*) = 1 *C*1(*z*) = *S*3

*C*3(*z*) = *S*5

**Proposition 7** (Edgeworth's expansion)**.** If we denote by *ϕ*(*z*) = exp(−*z*2/2)/

120 *H*5(*z*) + *S*3*S*4

**3.6 Near-Gaussian approximations**

at finite distance.

with

*H*�

the following approximation:

*k*−1(*z*) for all *k* 1.

 expectation std. dev. skewness e. kurtosis time (s) 12 1.280000 1.294320 1.163783 1.492661 0.01 120 15.104000 4.585724 0.361328 0.149974 0.02 1200 153.344000 14.648033 0.113920 0.014936 0.03 12000 1535.744000 46.367282 0.036014 0.001492 0.04<PHONE_NUMBER>9.744000 146.640798 0.011394 −0.000410 0.05

up with a *O*(log2 × *n*2 × *L*3) complexity in time where *L* is the order of the transition matrix **T** = **P** + **Q**. The corresponding memory complexity is *O*(log2 × *n* × *L*2). Since the length of the dataset appears in a logarithmic scale in these complexity, the power approach is obviously suitable for large datasets (ex: = 106 or = 109). Unfortunately, the cubic complexity with *L* (quadratic in memory) prevents the approach to deal with complex motifs with high *L*. One should also note that the quadratic complexity in *n* could really be a problem when dealing with frequent motifs and/or large datasets. In order to overcome this problem, Ribeca & Raineri (2008) suggested to use fast Fourier transforms (FFT) to perform all polynomial product hence replacing *n*2 by *n* log2 *n* in the time complexity. However appealing at first glance, this approach is not recommended in practice since the FFT products in floating-point arithmetics induce numerical instabilities that make totally unreliable the smallest coefficients of the polynomials. And unfortunately, these coefficients are precisely the one needed to study the tail distribution of *N*.

Another interesting approach called *full recursion*, consists in computing **v***i* = (**P** + *y***Q**)*i* **v** for all 0 *i* − *d* recursively using the relation **v***i*+1 = (**P** + *y***Q**)**v***i*. There are two main interests for this approach: 1) we have only products between polynomials of degree 1 and polynomials of degree *n* (by dropping terms of degree greater than *n* like in the power approach); 2) we can take full advantage of the sparse structure (only *L* × |A| non-zero terms in the worst case) of the transition matrix **T** = **P** + **Q**. The resulting complexity is *O*( × *L* × |A| × *n*) in time, and *O*(*L* × *n*) in memory. Since these complexities are linear with *L*, this approach is able to handle very complex motifs. The drawback is that the approach can be very slow when dealing with large and *n*. It exists a sophisticated version of this recursion called *partial recursion* (see Nuel & Dumas, 2010) which allows to replace × *n* by log × *n*2 in the time complexity. However, the quadratic complexity in *n* and numerical instabilities in floating-point arithmetic restrains its use to small *n* (ex: *n* 10).

For completeness, let us point out another approach to the problem. The idea is that we can derive from Expression (12) the following expression:

$$G(y, z) = \sum\_{n \ge 0} \sum\_{\ell \ge d} \mathbb{P}(N\_{\ell} = n) y^n z^{\ell} = \mathbf{u} z^d (\mathbf{I} - \mathbf{P}z + yz \mathbf{Q})^{-1} \mathbf{v} \tag{23}$$

where **I** is the identity matrix and *N* the number of matching position in *X*1:. It is then possible to obtain **P**(*N* = *n*) for any and *n* using (fast) Taylor expansions of *G*(*y*, *z*). For the mathematician, this approach is so "natural" that it is often referred as the "golden" approach to the problem of motif significance (Nicodème et al., 2002). However, this approach suffers several severe drawbacks that dramatically limits its practical interest: 1) the approach needs sophisticated computer algebra systems to be implemented (rather than simple floating point arithmetic for the previous approaches); 2) the explicit computation of (**I** − **P***z* + *yz***Q**)−1 could be very time (and memory) consuming; 3) even if the explicit computation of the inverse matrix is avoided (which is highly advisable), the coefficient extraction using state of the art techniques (like high-order lifting for example) is often slower than the much simpler alternative developed above (see Nuel & Dumas, 2010, for details).

Considering either the power or the recursion approaches we obtain easy to implement algorithms allowing to compute the exact value of *S*(*n*) in all cases except when dealing with high complexity motifs (large *L*) and/or frequent motifs (large *n*). But even if we stick to more tractable cases, exact computations could be slow. The question hence is: is it possible to compute fast and reliable approximations of *S*(*n*) ?

14 Will-be-set-by-IN-TECH

up with a *O*(log2 × *n*2 × *L*3) complexity in time where *L* is the order of the transition matrix **T** = **P** + **Q**. The corresponding memory complexity is *O*(log2 × *n* × *L*2). Since the length of the dataset appears in a logarithmic scale in these complexity, the power approach is obviously suitable for large datasets (ex: = 106 or = 109). Unfortunately, the cubic complexity with *L* (quadratic in memory) prevents the approach to deal with complex motifs with high *L*. One should also note that the quadratic complexity in *n* could really be a problem when dealing with frequent motifs and/or large datasets. In order to overcome this problem, Ribeca & Raineri (2008) suggested to use fast Fourier transforms (FFT) to perform all polynomial product hence replacing *n*2 by *n* log2 *n* in the time complexity. However appealing at first glance, this approach is not recommended in practice since the FFT products in floating-point arithmetics induce numerical instabilities that make totally unreliable the smallest coefficients of the polynomials. And unfortunately, these coefficients are precisely the one needed to study

Another interesting approach called *full recursion*, consists in computing **v***i* = (**P** + *y***Q**)*i*

all 0 *i* − *d* recursively using the relation **v***i*+1 = (**P** + *y***Q**)**v***i*. There are two main interests for this approach: 1) we have only products between polynomials of degree 1 and polynomials of degree *n* (by dropping terms of degree greater than *n* like in the power approach); 2) we can take full advantage of the sparse structure (only *L* × |A| non-zero terms in the worst case) of the transition matrix **T** = **P** + **Q**. The resulting complexity is *O*( × *L* × |A| × *n*) in time, and *O*(*L* × *n*) in memory. Since these complexities are linear with *L*, this approach is able to handle very complex motifs. The drawback is that the approach can be very slow when dealing with large and *n*. It exists a sophisticated version of this recursion called *partial recursion* (see Nuel & Dumas, 2010) which allows to replace × *n* by log × *n*2 in the time complexity. However, the quadratic complexity in *n* and numerical instabilities in floating-point arithmetic restrains

For completeness, let us point out another approach to the problem. The idea is that we can

where **I** is the identity matrix and *N* the number of matching position in *X*1:. It is then possible to obtain **P**(*N* = *n*) for any and *n* using (fast) Taylor expansions of *G*(*y*, *z*). For the mathematician, this approach is so "natural" that it is often referred as the "golden" approach to the problem of motif significance (Nicodème et al., 2002). However, this approach suffers several severe drawbacks that dramatically limits its practical interest: 1) the approach needs sophisticated computer algebra systems to be implemented (rather than simple floating point arithmetic for the previous approaches); 2) the explicit computation of (**I** − **P***z* + *yz***Q**)−1 could be very time (and memory) consuming; 3) even if the explicit computation of the inverse matrix is avoided (which is highly advisable), the coefficient extraction using state of the art techniques (like high-order lifting for example) is often slower than the much simpler

Considering either the power or the recursion approaches we obtain easy to implement algorithms allowing to compute the exact value of *S*(*n*) in all cases except when dealing with high complexity motifs (large *L*) and/or frequent motifs (large *n*). But even if we stick to more tractable cases, exact computations could be slow. The question hence is: is it possible

**P**(*N* = *n*)*ynz* = **u***zd*(**I** − **P***z* + *yz***Q**)−1**v** (23)

**v** for

the tail distribution of *N*.

its use to small *n* (ex: *n* 10).

derive from Expression (12) the following expression:

*n*0 ∑ *d*

alternative developed above (see Nuel & Dumas, 2010, for details).

to compute fast and reliable approximations of *S*(*n*) ?

*G*(*y*, *z*) = ∑

Table 1. Characteristic moments the number *N* of occurrences of motif M = G(G|C)G in a sequence *X*1: generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40. Computation performed using the power approach.

#### **3.6 Near-Gaussian approximations**

Since the random count *N* is basically defined by Eq. (2) as large sum of Bernouilli variables, the idea of approximating the distribution of *N* using Gaussian approximation sounds appealing. Indeed, Gaussian approximations are historically the first ones to have been suggested for this problem (Cowan, 1991; Kleffe & Borodovski, 1997; Pevzner et al., 1989; Prum et al., 1995). From the theoretical point of view, Central Limit Theorems (CLT) for weakly dependent variables ensure that *N* is asymptotically normal distributed. On Table 1, we can see the characteristic moments of *N* for motif M = G(G|C)G and various value of the sequence lengths . According to theory, we observe that the skewness and excess kurtosis both decease toward 0 when grows (a normal distribution has null skewness and excess kurtosis). But it is also clear that *N* is not normally distributed for small values of . As a consequence, the quality of a Gaussian approximation for *S*(*n*) is expected to be questionable at finite distance.

In order to overcome this issue, Nuel (2010) suggested to consider near Gaussian approximations instead of simple Gaussian approximations for this problem. The idea is simply to perform a higher order asymptotic development that exploits more than the two first moments of *N*. This technique is known as the Edgeworth's expansion. Blinnikov & Moessner (1998) gives a general (and rather complicated) formula for this expansion. For practical purpose, we present the result only up to order 3 expansions.

**Proposition 7** (Edgeworth's expansion)**.** If we denote by *ϕ*(*z*) = exp(−*z*2/2)/ √2*π* the probability distribution function (pdf) of a standard Gaussian, then for all *n* 0 we have the following approximation:

$$\mathbb{P}(N=n) \simeq \frac{\varrho(z)}{\sigma} \left( \mathbb{C}\_0(z) + \sigma \mathbb{C}\_1(z) + \sigma^2 \mathbb{C}\_2(z) + \sigma^3 \mathbb{C}\_3(z) \right) \tag{24}$$

with

$$\mathbb{C}\_{0}(z) = 1 \quad \mathbb{C}\_{1}(z) = \frac{\mathbb{S}\_{3}}{6} H\_{3}(z) \quad \mathbb{C}\_{2}(z) = \frac{\mathbb{S}\_{4}}{24} H\_{4}(z) + \frac{\mathbb{S}\_{3}^{2}}{72} H\_{6}(z) \tag{25}$$

$$\mathcal{C}\_{3}(z) = \frac{S\_{5}}{120}H\_{7}(z) + \frac{S\_{3}S\_{4}}{144}H\_{7}(z) + \frac{S\_{3}^{3}}{1296}H\_{9}(z) \tag{26}$$

where *μ* = **E**[*N*], *σ* = **V**[*N*], *z* = (*n* − *μ*)/*σ*, *Sk* = *κk*/*σ*2*k*−2 for all *k* 1, and where *Hk*(*z*) are the Hermite polynomials defined recursively by *H*0(*z*) = 1 and *Hk*(*z*) = *zHk*−1(*z*) − *H*� *k*−1(*z*) for all *k* 1.

in Biological Sequences 17

Significance Score of Motifs in Biological Sequences 189

−10

Fig. 7. Reliability of CB and BR approximations for M = G(G|C)G on a random sequence *X*1:� generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40, and with � = 1200. The error CB(*n*) − *S*(*n*) or BR(*n*) − *S*(*n*) is given on Figure (a); and the relative error (log-scale) − log10 |CB(*n*) − *S*(*n*)|/|*S*(*n*)| or − log10 |BR(*n*) − *S*(*n*)|/|*S*(*n*)| on Figure (b). The horizontal rule indicates the null error on Figure (a), and the threshold

From the computational point of view, the order *h* approximation requires the cumulants of *N* up to order *h* + 2. Using the power approach, the resulting complexity is hence *O*(log2 � × *h*2 × *L*3) in time and *O*(log2 � × (*h* + 2) × *L*2) in memory. Using the recursion, the complexity resulting complexity is *O*(� × *L* × |A| × *h*) in time, and *O*(*L* × *h*) in memory. In both cases,

Thanks to NG approximations, we hence have a fast and reliable way to compute an approximation of *S*(*n*) when *n* falls in the center of the distribution (ex: |*S*(*n*)| 3.0), but NG approximations unfortunately remain totally unreliable for tail distribution events (ex: |*S*(*n*)| > 3.0), which are moreover often precisely the event of interest. Fortunately we have a

We want here to study specifically the tail distribution of *N* with events on the form **P**(*N n*) with large *n* (or **P**(*N n*) with small *n*). For all *t* > 0 let us first notice that we can use the Markov inequality to write: **P**(*N n*) = **P**(*etN etn*) **E**[*etN*]/*etn* = exp(Λ(*t*) − *tn*). By taking the smallest of these bounds for *t* > 0 we hence get: log **P**(*N n*) Λ(*τ*) − *τn* with

surprisingly sharp for events located in the tail distribution. By dealing symmetrically with

(*τ*) = *n*. This upper bound, known as the Chernoff's Bound (CB), is often

log(10) (28)

*τn* − Λ(*τ*)

the computational time drops significantly from the exact computations.

**P**(*N n*) and *t* < 0 we hence obtain the following approximation for *S*(*n*):

where *δn* = −1 if *n* **E**[*N*], and *δn* = +1 if *n* > **E**[*N*].

CB(*n*) = *δn*

−log10(relative error)

0 100 200 300 400

(b) Relative error

Cramer Bahadur−Rao

n

0 100 200 300 400

(a) Error

corresponding to two correct digits on Figure (b).

n

Cramer Bahadur−Rao

−1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

solution to this problem.

**3.7 Bahadur-Rao**

*τ* defined by Λ�

error

Fig. 6. Reliability of NG approximations for M = G(G|C)G on a random sequence *X*1:� generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40, and with � = 1200. The error NG*h*(*n*) − *S*(*n*) is given on Figure (a); and the relative error (log-scale) − log10 |NG*h*(*n*) − *S*(*n*)|/|*S*(*n*)| on Figure (b). The horizontal rule indicates the null error on Figure (a), and the threshold corresponding to two correct digits on Figure (b).

For *h* ∈ {0, 1, 2, 3} we define the Near Gaussian (NG) approximation of order *h* of *S*(*n*) by:

$$\mathrm{NG}\_{h}(n) = \begin{cases} +\log\_{10}\left(\sum\_{k=0}^{n}\frac{1}{\sigma}\varphi\left(\frac{k-\mu}{\sigma}\right)\sum\_{j=0}^{h}\sigma^{j}\mathbb{C}\_{j}\left(\frac{k-\mu}{\sigma}\right)\right) & \text{if } n \leqslant \mathbb{E}[N] \\\ -\log\_{10}\left(\sum\_{k=n}^{+\infty}\frac{1}{\sigma}\varphi\left(\frac{k-\mu}{\sigma}\right)\sum\_{j=0}^{h}\sigma^{j}\mathbb{C}\_{j}\left(\frac{k-\mu}{\sigma}\right)\right) & \text{if } n > \mathbb{E}[N] \end{cases} . \tag{27}$$

We can see on Figure 6 the reliability of NG approximations. In solid black, the order 0 approximation corresponds to the classical Gaussian approximation. Unsurprisingly, this central limit approximation is accurate for the center of the distribution (*n* close to the expectation *μ* = 153.3), the reliability quickly deceases when |*n* − *μ*| increases.

Central limit theorems (CLT) for *N* have established long ago that *N* should be asymptotically Gaussian distributed. The problem however with CLT theorems is that the quality of the resulting approximation dramatically decreases at finite distance when considering tail distribution events. Here we try to overcome the issue by considering Near-Gaussian approximations that exploits higher moments of *N* to improve the quality of the approximations. In order to do this, a critical problem is first to obtain the first *k*-th moments of *N*. Of course we can access these moments by computing the full distribution of *N*, but if it is possible to do so, why bothering with approximations. We hence need an method to compute the moments of *N* whose complexity should be somehow significantly smaller than the complete exact computations. With higher order approximation, we can see a dramatic improvement of reliability of the results, with a noticeable increase of the region where at least two digits are correct (up to *n* ∈ [80; 240] for NG3).

16 Will-be-set-by-IN-TECH

02468

−log10(relative error)

Fig. 6. Reliability of NG approximations for M = G(G|C)G on a random sequence *X*1:� generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40, and

with � = 1200. The error NG*h*(*n*) − *S*(*n*) is given on Figure (a); and the relative error (log-scale) − log10 |NG*h*(*n*) − *S*(*n*)|/|*S*(*n*)| on Figure (b). The horizontal rule indicates the null error on Figure (a), and the threshold corresponding to two correct digits on Figure (b).

For *h* ∈ {0, 1, 2, 3} we define the Near Gaussian (NG) approximation of order *h* of *S*(*n*) by:

� *h* ∑ *j*=0 *σj Cj*

� *h* ∑ *j*=0 *σj Cj*

� *k* − *μ σ*

� *k* − *μ σ*

We can see on Figure 6 the reliability of NG approximations. In solid black, the order 0 approximation corresponds to the classical Gaussian approximation. Unsurprisingly, this central limit approximation is accurate for the center of the distribution (*n* close to the

Central limit theorems (CLT) for *N* have established long ago that *N* should be asymptotically Gaussian distributed. The problem however with CLT theorems is that the quality of the resulting approximation dramatically decreases at finite distance when considering tail distribution events. Here we try to overcome the issue by considering Near-Gaussian approximations that exploits higher moments of *N* to improve the quality of the approximations. In order to do this, a critical problem is first to obtain the first *k*-th moments of *N*. Of course we can access these moments by computing the full distribution of *N*, but if it is possible to do so, why bothering with approximations. We hence need an method to compute the moments of *N* whose complexity should be somehow significantly smaller than the complete exact computations. With higher order approximation, we can see a dramatic improvement of reliability of the results, with a noticeable increase of the region where at

50 100 150 200 250 300

order 0 order 1 order 2 order 3

n

(b) Relative error

� *k* − *μ σ*

� *k* − *μ σ*

�⎞

�⎞

⎠ if *n* **E**[*N*]

. (27)

⎠ if *n* > **E**[*N*]

50 100 150 200 250 300

order 0 order 1 order 2 order 3

n

(a) Error

−10

NG*h*(*n*) =

⎧ ⎪⎪⎪⎪⎪⎪⎨

+ log10

− log10

least two digits are correct (up to *n* ∈ [80; 240] for NG3).

⎛ ⎝ *n* ∑ *k*=0

⎛ ⎝ +∞ ∑ *k*=*n* 1 *σ ϕ*

1 *σ ϕ*

expectation *μ* = 153.3), the reliability quickly deceases when |*n* − *μ*| increases.

⎪⎪⎪⎪⎪⎪⎩

error

Fig. 7. Reliability of CB and BR approximations for M = G(G|C)G on a random sequence *X*1:� generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40, and with � = 1200. The error CB(*n*) − *S*(*n*) or BR(*n*) − *S*(*n*) is given on Figure (a); and the relative error (log-scale) − log10 |CB(*n*) − *S*(*n*)|/|*S*(*n*)| or − log10 |BR(*n*) − *S*(*n*)|/|*S*(*n*)| on Figure (b). The horizontal rule indicates the null error on Figure (a), and the threshold corresponding to two correct digits on Figure (b).

From the computational point of view, the order *h* approximation requires the cumulants of *N* up to order *h* + 2. Using the power approach, the resulting complexity is hence *O*(log2 � × *h*2 × *L*3) in time and *O*(log2 � × (*h* + 2) × *L*2) in memory. Using the recursion, the complexity resulting complexity is *O*(� × *L* × |A| × *h*) in time, and *O*(*L* × *h*) in memory. In both cases, the computational time drops significantly from the exact computations.

Thanks to NG approximations, we hence have a fast and reliable way to compute an approximation of *S*(*n*) when *n* falls in the center of the distribution (ex: |*S*(*n*)| 3.0), but NG approximations unfortunately remain totally unreliable for tail distribution events (ex: |*S*(*n*)| > 3.0), which are moreover often precisely the event of interest. Fortunately we have a solution to this problem.

#### **3.7 Bahadur-Rao**

We want here to study specifically the tail distribution of *N* with events on the form **P**(*N n*) with large *n* (or **P**(*N n*) with small *n*). For all *t* > 0 let us first notice that we can use the Markov inequality to write: **P**(*N n*) = **P**(*etN etn*) **E**[*etN*]/*etn* = exp(Λ(*t*) − *tn*). By taking the smallest of these bounds for *t* > 0 we hence get: log **P**(*N n*) Λ(*τ*) − *τn* with *τ* defined by Λ� (*τ*) = *n*. This upper bound, known as the Chernoff's Bound (CB), is often surprisingly sharp for events located in the tail distribution. By dealing symmetrically with **P**(*N n*) and *t* < 0 we hence obtain the following approximation for *S*(*n*):

$$\text{CB}(n) = \delta\_n \frac{\tau n - \Lambda(\tau)}{\log(10)} \tag{28}$$

where *δn* = −1 if *n* **E**[*N*], and *δn* = +1 if *n* > **E**[*N*].

in Biological Sequences 19

Significance Score of Motifs in Biological Sequences 191

−10

Fig. 8. Relative error in log-scale for various approximations of *S*(*n*) (*n* = 0, . . . , 200) in a sequence *X*1: generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and

In this chapter we deliberately left aside the Poisson-based approximations and considered only two of these approximations: the (Near-) Gaussian approximations with NG*h*(*n*), and the large deviations based approximations with CB(*n*) and BR(*n*). The reason why Poisson-based approximations are not considered here is basically practical, these approximations cannot be directly derived from the formalism of this manuscript and require the introduction of many tedious notions like clumps, overlapping words and so on. However, we compare here the performance of all these approximations (including compound Poisson approximations) in the case where *X*1: generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40 i.i.d. DNA sequence, and for two motifs: the frequent G(G|C)G, and the

We can see on Figure 8 the relative error (in log-scale) for all approximations. For Gaussian approximations, performances are only good in the very center of the distribution (for *n* very close to **E**(*n*)) for the frequent motif G(G|C)G, and performances are poor almost everywhere for the rare motif T(A|T)T. This observation to consistent with the well known claim that "Gaussian approximations a more suitable for frequent motif" (Lothaire, 2005). It has however to be pointed out that even in the most favorable case (with highly frequent motif), Gaussian approximations totally fail to capture the tail distribution of *N* and hence not suitable for the highly significant observations we usually encounter in biological sequences (Nuel, 2006b). If we consider now the near-Gaussian approximation, taking into account more moments of *N* dramatically improve the result for both motifs, but the failure to deal with extreme

Compound Poisson approximations are known to be extremely sensitive to the relative abundance of the motif of interest in the sequence, being more accurate for rare motifs (Lothaire, 2005; Roquain & Schbath, 2007). It is hence not a surprise to see that Poisson approximations are totally unreliable for the frequent motif G(G|C)G. For the rare motif T(A|T)T we naturally obtain much better results but like for Gaussian approximations, and even in this favorable case, reliability decreases in the tail distribution. Considering that Poisson

−log10(relative error)

0 20 40 60

S(n)

(b) A(A|T)A, rare motif

Gaussian Near−Gaussian Poisson Cramer Bahadur−Rao

−40 −20 0 20 40

S(n)

(a) G(G|C)G, frequent motif

02468

*π*(C) = *π*(G) = 0.40.

rare A(A|T)A.

distribution events remains.

Gaussian Near−Gaussian Poisson Cramer Bahadur−Rao

−log10(relative error)

From the computational point of view, the solution *τ* of Λ� (*τ*) = *n* can be easily determined numerically using (for example) using the Newton-Raphson sequence (Press et al., 1992). Starting for a first guess *t*0 (ex: *t*0 = 0), one performs *ti*+1 = *ti* + (*n* − Λ� (*ti*))/Λ��(*ti*) for *i* 0 until convergence to *τ*. The computation of Λ, Λ� , and Λ�� being possible thanks to Lemma 5 and the following formulas:

$$
\Lambda(t) = G(\varepsilon^t) \quad \Lambda'(t) = \frac{\varepsilon^t G'(\varepsilon^t)}{G(\varepsilon^t)} \quad \Lambda''(t) = \frac{\varepsilon^{2t} G''(\varepsilon^t)}{G(\varepsilon^t)} - \frac{\varepsilon^{2t} G'(t)^2}{G(\varepsilon^t)^2} + \frac{\varepsilon^t G'(\varepsilon^t)}{G(\varepsilon^t)} \tag{29}
$$

with *G*(*et* )=[*z*0]**u**(**P** + *et* **Q** + *z***Q**)−*d***v** = **u**(**P** + *et* **Q**)−*d***v**, *G*� (*et* )=[*z*1]**u**(**P** + *et* **Q** + *z***Q**)−*d***v**, and *G*��(*et* ) = 2[*z*2]**u**(**P** + *et* **Q** + *z***Q**)−*d***v**.

Moreover, this bound can be further refined using the Bahadur-Rao Theorem (Bahadur & Rao, 1960) and gives the following approximation for *S*(*n*):

$$\text{BR}(n) = \text{CB}(n) + \delta\_n \log 10 \left( (1 - e^{-|\tau|}) \sqrt{2\pi \Lambda''(\tau)} \right). \tag{30}$$

From the computational point of view, CB(*n*) and BR(*n*) can be computed either with the power approach with complexities *O*(log2 × *L*3) in time and *O*(log2 × *L*3) in memory; or with the recursion approach with complexities *O*( × *L* × |A|) in time and *O*(*L* × |A|) in memory.

On Figure 7 we can see the reliability of the approximations CB(*n*) and BR(*n*). Unsurprisingly, the farther from the center of the distribution, the better are both approximations. We also observe that BR(*n*) is a dramatic improvement over CB(*n*) since it obtains at least two correct digits of *S*(*n*) for all *n* but on [120, 200]. At the end of previous section, we have seen that the order 3 NG approximation achieves the same precision for region [80; 240], hence, by combining both NG3(*n*) (for the center of the distribution) and BR(*n*) (for the tail distributions), one can achieve at least two correct digits of *S*(*n*) on the whole bulk of the distribution for a modest computational cost.

#### **4. Discussion**

Obtaining the distribution of motif count in random sequences is a very challenging problem that has attracted considerable attention from mathematicians and computer scientists in the last fifty years. Recently however, a significant advance has been obtained by connecting the well-known theory of pattern matching and automata to the Markov chain embedding technique Lladser (2007); Nuel (2008a); Nuel & Prum (2007). Thanks to this finding, it is now possible to deal with simple (runs of 1 in binary sequences, single words, etc.) or complex motifs (PROSITE signature, gapped motifs, etc.) using the same general framework.

Using exact approaches, it is possible to obtain efficiently the first moments of any motif count *N*, and even the complete distribution of *N*. As a consequence, the computation of *S*(*n*) is now tractable for a wide range of motif problems including large datasets or complex motifs. However, the case of complex frequent motifs in large datasets remains an open problem (Nuel & Dumas, 2010).

As an alternative to exact computations, a wide range of approximations have been developed (see Lothaire, 2005; Nuel, 2006b; Reignier, 2000, for a review). We can basically classify these approximations in three categories: 1) Gaussian approximations (Cowan, 1991; Kleffe & Borodovski, 1997; Nuel, 2010; Pevzner et al., 1989; Prum et al., 1995); 2) Poisson approximations Erhardsson (2000); Geske et al. (1995); Godbole (1991); Reinert & Schbath (1999); Roquain & Schbath (2007); 3) large deviations approximations Denise et al. (2001); Nuel (2004).

18 Will-be-set-by-IN-TECH

numerically using (for example) using the Newton-Raphson sequence (Press et al., 1992).

Moreover, this bound can be further refined using the Bahadur-Rao Theorem (Bahadur & Rao,

 (1 − *e*

From the computational point of view, CB(*n*) and BR(*n*) can be computed either with the power approach with complexities *O*(log2 × *L*3) in time and *O*(log2 × *L*3) in memory; or with the recursion approach with complexities *O*( × *L* × |A|) in time and *O*(*L* × |A|) in

On Figure 7 we can see the reliability of the approximations CB(*n*) and BR(*n*). Unsurprisingly, the farther from the center of the distribution, the better are both approximations. We also observe that BR(*n*) is a dramatic improvement over CB(*n*) since it obtains at least two correct digits of *S*(*n*) for all *n* but on [120, 200]. At the end of previous section, we have seen that the order 3 NG approximation achieves the same precision for region [80; 240], hence, by combining both NG3(*n*) (for the center of the distribution) and BR(*n*) (for the tail distributions), one can achieve at least two correct digits of *S*(*n*) on the whole bulk of the

Obtaining the distribution of motif count in random sequences is a very challenging problem that has attracted considerable attention from mathematicians and computer scientists in the last fifty years. Recently however, a significant advance has been obtained by connecting the well-known theory of pattern matching and automata to the Markov chain embedding technique Lladser (2007); Nuel (2008a); Nuel & Prum (2007). Thanks to this finding, it is now possible to deal with simple (runs of 1 in binary sequences, single words, etc.) or complex

Using exact approaches, it is possible to obtain efficiently the first moments of any motif count *N*, and even the complete distribution of *N*. As a consequence, the computation of *S*(*n*) is now tractable for a wide range of motif problems including large datasets or complex motifs. However, the case of complex frequent motifs in large datasets remains an open problem

As an alternative to exact computations, a wide range of approximations have been developed (see Lothaire, 2005; Nuel, 2006b; Reignier, 2000, for a review). We can basically classify these approximations in three categories: 1) Gaussian approximations (Cowan, 1991; Kleffe & Borodovski, 1997; Nuel, 2010; Pevzner et al., 1989; Prum et al., 1995); 2) Poisson approximations Erhardsson (2000); Geske et al. (1995); Godbole (1991); Reinert & Schbath (1999); Roquain & Schbath (2007); 3) large deviations approximations Denise et al. (2001);

motifs (PROSITE signature, gapped motifs, etc.) using the same general framework.

*G*��(*et* ) *G*(*et*) − *e*2*t*

(*et*

2*π*Λ��(*τ*)

**Q**)−*d***v**, *G*�

−|*τ*| )

Starting for a first guess *t*0 (ex: *t*0 = 0), one performs *ti*+1 = *ti* + (*n* − Λ�

*G*� (*et* ) *G*(*et*) Λ��(*t*) = *e*2*t*

**Q** + *z***Q**)−*d***v**.

BR(*n*) = CB(*n*) + *δn* log 10

**Q** + *z***Q**)−*d***v** = **u**(**P** + *et*

(*τ*) = *n* can be easily determined

, and Λ�� being possible thanks to

*et G*� (*et* ) *G*(*et*) (29)

*G*� (*t*)2 *G*(*et*)2 +

)=[*z*1]**u**(**P** + *et*

(*ti*))/Λ��(*ti*) for

**Q** + *z***Q**)−*d***v**,

. (30)

From the computational point of view, the solution *τ* of Λ�

*i* 0 until convergence to *τ*. The computation of Λ, Λ�

(*t*) = *et*

1960) and gives the following approximation for *S*(*n*):

distribution for a modest computational cost.

Lemma 5 and the following formulas:

*t* ) Λ�

)=[*z*0]**u**(**P** + *et*

) = 2[*z*2]**u**(**P** + *et*

Λ(*t*) = *G*(*e*

with *G*(*et*

and *G*��(*et*

memory.

**4. Discussion**

(Nuel & Dumas, 2010).

Nuel (2004).

Fig. 8. Relative error in log-scale for various approximations of *S*(*n*) (*n* = 0, . . . , 200) in a sequence *X*1: generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40.

In this chapter we deliberately left aside the Poisson-based approximations and considered only two of these approximations: the (Near-) Gaussian approximations with NG*h*(*n*), and the large deviations based approximations with CB(*n*) and BR(*n*). The reason why Poisson-based approximations are not considered here is basically practical, these approximations cannot be directly derived from the formalism of this manuscript and require the introduction of many tedious notions like clumps, overlapping words and so on. However, we compare here the performance of all these approximations (including compound Poisson approximations) in the case where *X*1: generated by a M0 model with parameters *π*(A) = *π*(T) = 0.10 and *π*(C) = *π*(G) = 0.40 i.i.d. DNA sequence, and for two motifs: the frequent G(G|C)G, and the rare A(A|T)A.

We can see on Figure 8 the relative error (in log-scale) for all approximations. For Gaussian approximations, performances are only good in the very center of the distribution (for *n* very close to **E**(*n*)) for the frequent motif G(G|C)G, and performances are poor almost everywhere for the rare motif T(A|T)T. This observation to consistent with the well known claim that "Gaussian approximations a more suitable for frequent motif" (Lothaire, 2005). It has however to be pointed out that even in the most favorable case (with highly frequent motif), Gaussian approximations totally fail to capture the tail distribution of *N* and hence not suitable for the highly significant observations we usually encounter in biological sequences (Nuel, 2006b). If we consider now the near-Gaussian approximation, taking into account more moments of *N* dramatically improve the result for both motifs, but the failure to deal with extreme distribution events remains.

Compound Poisson approximations are known to be extremely sensitive to the relative abundance of the motif of interest in the sequence, being more accurate for rare motifs (Lothaire, 2005; Roquain & Schbath, 2007). It is hence not a surprise to see that Poisson approximations are totally unreliable for the frequent motif G(G|C)G. For the rare motif T(A|T)T we naturally obtain much better results but like for Gaussian approximations, and even in this favorable case, reliability decreases in the tail distribution. Considering that Poisson

in Biological Sequences 21

Significance Score of Motifs in Biological Sequences 193

Cornish-Bowden (1985). IUPAC-IUB symbols for nucleotide nomenclature, *Nucl. Acids Res.*

Cowan (1991). Expected frequencies of dna patterns using whittle's formula, *J. Appl. Prob.*

Crochemore, M. & Stefanov, V. (2003). Waiting time and complexity for matching patterns

Denise, A., RA©gnier, M. & Vandenbogaert, M. (2001). Assessing the statistical significance of overrepresented oligonucleotides, *Lecture Notes in Computer Science* 2149: 85–97. El Karoui, M., Biaudet, V., Schbath, S. & Gruss, A. (1999). Characteristics of chi distribution

Erhardsson, T. (2000). Compound Poisson approximation for counts of rare patterns in

Fatemi, M., Pao, M., Jeong, S., Gal-Yam, E., Egger, G., Weisenberger, D. & Jones, P. (2005).

Fu, J. C. (1996). Distribution theory of runs and patterns associated with a sequence of

Geske, M. X., Godbole, A. P., Schaffner, A. A., Skrolnick, A. M. & Wallstrom, G. L. (1995).

Godbole, A. P. (1991). Poissons approximations for runs and patterns of rare events, *Adv. Appl.*

Green, T. J., Gupta, A., Miklau, G., Onizuka, M. & Suciu, D. (2004). Processing xml

Hampson, S., Kibler, D. & Baldi, P. (2002). Distribution patterns of over-represented k-mers in

Hopcroft, J. E., Motwani, R. & Ullman, J. D. (2001). *Introduction the automata theory, languages,*

Karlin, S., Burge, C. & Campbell, A. (1992). Statistical analyses of counts and distributions of restriction sites in DNA sequences, *Nucl. Acids. Res.* 20(6): 1363–1370. Kleffe, J. & Borodovski, M. (1997). First and second moment of counts of words in random

Leonardo Marino-Ramírez, John L. Spouge, G. C. K. & Landsman, D. (2004). Statistical

Liddle, A. R. (2007). Information criteria for astrophysical model selection, *Monthly Notices of*

Lladser, M. E. (2007). Mininal markov chain embeddings of pattern problems, *Information*

Lothaire, M. (ed.) (2005). *Applied Combinatorics on Words*, Cambridge University Press,

Nicodème, P., Salvy, B. & Flajolet, P. (2002). Motif statistics, *Theoretical Com. Sci.* 287(2): 593–617. Nuel, G. (2004). Ld-spatt: Large deviations statistics for patterns on markov chains, *J. Comp.*

analysis of over-represented words in human promoter sequences, *Nuc. Acids Res.*

Markov chains and extreme sojourns in birth-death chains, *Ann. Appl. Probab.*

Footprinting of mammalian promoters: use of a cpg dna methyltransferase revealing nucleosome positions at a single molecule level, *Nucleic Acids Res* 33(20): 176. Frith, M. C., Spouge, J. L., Hansen, U. & Weng, Z. (2002). Statistical significance of clusters

of motifs represented by position specific scoring matrices in nucleotide sequences,

Compound poisson approximations for word patterns under markovian hypotheses,

streams with deterministic automata and stream indexes, *ACM Trans. Database Syst.*

with automata, *Info. Proc. Letters* 87(3): 119–125.

on different bacterial genomes, *Res. Microbiol.* 150: 579–587.

13: 3021–3030.

10(2): 573–591.

*Prob.* 23.

29: 752–788.

32(3): 949–958.

Cambridge.

*Biol.* 11(6): 1023–1033.

*Nucl. Acids. Res.* 30(14): 3214–3224.

*J. Appl. Probab.* 32: 877–892.

multi-state trials, *Statistica Sinica* 6(4): 957–974.

non-coding yeast DNA, *Bioinformatics* 18(4): 513–528.

texts generated by markov chains, *Bioinformatics* 8(5): 433–441.

*and computation, 2d ed.*, ACM Press, New York.

*the Royal Astronomical Society: Letters* 377: 74–78.

*Theory and Applications Workshop*, pp. 251–255.

28: 886–892.

approximations are not easily generalizable to motifs defined by regular expressions, that their computations could be complicated and time consuming, and that their reliability is highly questionable in some configurations, it seems advisable to avoid their use is most cases.

With large deviations based approximations, we unsurprisingly get a low reliability in the center of the distribution, but a high reliability in the tail distribution. With Bahadur-Rao precise approximations, the improvement over the classical Chernoff's bound is quite impressive, and the complementarity with Near-Gaussian approximations clearly shows that a combination of both approaches could be a very efficient way to obtain reliable approximations of *S*(*n*) for all *n*.

In this chapter we gave all the necessary ingredients to assess the significance score of motif in a biological sequence using state of the art results, including several unpublished ones: Lemma 5 which is an extension of the results of Nuel (2010), and the complete "Bahadur-Rao" Section which provides interesting improvements over previous large deviations work (Denise et al., 2001; Nuel, 2004).

Let us finally point out that for the sake of compactness, we have left aside some interesting questions and extensions like: approximate matching Hopcroft et al. (2001), renewal occurrences (Nuel, 2006b; Roquain & Schbath, 2007), joint distributions (Nuel, 2008b; Stefanov & Szpankowski, 2007), dataset with many sequences (Nuel et al., 2010), and sensitivity to parameter estimation (Nuel, 2006c). Even if some results are already available for these problems, many questions still have to be answered in the exciting and challenging field of the distribution of motifs in random sequences.

#### **5. References**

20 Will-be-set-by-IN-TECH

approximations are not easily generalizable to motifs defined by regular expressions, that their computations could be complicated and time consuming, and that their reliability is highly questionable in some configurations, it seems advisable to avoid their use is most cases. With large deviations based approximations, we unsurprisingly get a low reliability in the center of the distribution, but a high reliability in the tail distribution. With Bahadur-Rao precise approximations, the improvement over the classical Chernoff's bound is quite impressive, and the complementarity with Near-Gaussian approximations clearly shows that a combination of both approaches could be a very efficient way to obtain reliable

In this chapter we gave all the necessary ingredients to assess the significance score of motif in a biological sequence using state of the art results, including several unpublished ones: Lemma 5 which is an extension of the results of Nuel (2010), and the complete "Bahadur-Rao" Section which provides interesting improvements over previous large deviations work

Let us finally point out that for the sake of compactness, we have left aside some interesting questions and extensions like: approximate matching Hopcroft et al. (2001), renewal occurrences (Nuel, 2006b; Roquain & Schbath, 2007), joint distributions (Nuel, 2008b; Stefanov & Szpankowski, 2007), dataset with many sequences (Nuel et al., 2010), and sensitivity to parameter estimation (Nuel, 2006c). Even if some results are already available for these problems, many questions still have to be answered in the exciting and challenging field of

Allauzen, C. & Mohri, M. (2006). A unified construction of the glushkov, follow, and

Antzoulakos, D. L. (2001). Waiting times for patterns in a sequence of multistate trials, *J. Appl.*

Bahadur, R. R. & Rao, R. R. (1960). On deviations of the sample mean, *The Annals of Math.*

Beaudoing, E., Freier, S., Wyatt, J., Claverie, J.-M. & Gautheret, D. (2000). Patterns of variant polyadenylation signal usage in human genes, *Genome Res.* 10(7): 1001–1010. Blinnikov, S. & Moessner, R. (1998). Expansions for nearly Gaussian distributions, *Astron.*

Boeva, V., ClA©ment, J., RA©gnier, M. & Vandenbogaert, M. (2005). Assessing the significance

Brazma, A., Jonassen, I., Vilo, J. & Ukkonen, E. (1998). Predicting gene regulatory elements in

Bryne, J., Valen, E., Tang, M., Marstrand, T., Winther, O., da Piedade, I., Krogh, A.,

Chan, H. P., Zhang, N. R. & Chen, L. H. Y. (2010). Importance sampling of word patterns in

Chang, Y.-M. (2005). Distribution of waiting time until the rth occurrence of a compound

dna and protein sequences, *J. of Comput. Biol.* 17(12): 1697–1709.

silico on a genomic scale, *Genome Res.* 8(11): 1202–1215.

pattern, *Statistics and Probability Letters* 75(1): 29–38.

of sets of words, *Combinatorial Pattern Matching 05, Lecture Notes in Computer Science,*

Lenhard, B. & Sandelin, A. (2008). JASPAR, the open access database of transcription factor-binding profiles: new content and tools in the 2008 update., *Nucleic Acids Res.*

antimirov automata, *in* R. KrA?lovic & P. Urzyczyn (eds), *Mathematical Foundations of Computer Science 2006*, Vol. 4162 of *Lecture Notes in Computer Science*, Springer Berlin

approximations of *S*(*n*) for all *n*.

(Denise et al., 2001; Nuel, 2004).

**5. References**

the distribution of motifs in random sequences.

/ Heidelberg, pp. 110–121.

*Statistics.* 31(4): 1015–1027.

*vol. 3537*, Springer-Verlag.

36: 102–106.

*Astrophys. Suppl. Ser.* 130: 193–205.

*Prob.* 38: 508–518.

**10**

*1Portugal 2UK*

**A Systematic and Thorough Search for Domains**

*1IBMC – Instituto de Biologia Molecular e Celular and ICBAS – Instituto de Ciências*

The biological function of proteins is largely determined by their individual component domains, which are segments within the protein sequence that are self-contained and spatially arranged. These can be catalytic or structural, and define a number of different features of proteins such as their enzymatic activity, interactions with other proteins, sugars or lipids, and determine the cellular localization of the proteins that contain them. A number of intracellular three-dimensionally-arranged domains, such as Src-homology (SH) or Pleckstrin-homology (PH) domains, define the nature of protein interactions with other components of the cell, and enable them to interact with their substrates or binding partners. The specificity of interactions that is given by the domain is unique to its protein. Similarly, the extracellular part of most membrane-bound or secreted proteins of eukaryotic cells is also organized in semi-autonomously-arranged blocks that potentially confer multiple diverse functions to a particular protein. These domains have been classified and grouped into protein superfamilies depending on the similarity they have with domains of prototypical proteins, for example immunoglobulin, fibronectin or C-type lectin domains. Members of these groups are believed to be homologous and to have arisen by divergent evolution from a common ancestor. Many membrane-bound or extracellular proteins are comprised of several domains of the same type, but it is not uncommon to find mosaic

The scavenger receptor cysteine-rich (SRCR) superfamily comprises a group of proteins that contain one or multiple domains structurally similar to the membrane distal domain of the type I scavenger receptor expressed by human macrophages (Freeman et al., 1990). Proteins classified as belonging to this superfamily may contain other types of domains additionally to the dominant SRCR modules, such as EGF, CUB, LCCL, or other domains. In mammals, SRCR proteins are typically expressed in cells of the immune system (Resnick et al., 1994), although some members can be also expressed in non-immune cells and organs, including liver, kidney, placenta, stomach, brain and heart (Sarrias et al., 2004). Group A domaincontaining SRCR proteins are present in phyla from the most primitive metazoan to

proteins containing domains from different superfamilies.

**1. Introduction**

**of the Scavenger Receptor Cysteine-Rich**

**Group-B Family in the Human Genome**

*2The Weatherall Institute of Molecular Medicine, University of Oxford*

Alexandre M. Carmo1 and Vattipally B. Sreenu2

*Biomédicas de Abel Salazar, Universidade do Porto*

### **A Systematic and Thorough Search for Domains of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome**

Alexandre M. Carmo1 and Vattipally B. Sreenu2

*1IBMC – Instituto de Biologia Molecular e Celular and ICBAS – Instituto de Ciências Biomédicas de Abel Salazar, Universidade do Porto 2The Weatherall Institute of Molecular Medicine, University of Oxford 1Portugal 2UK*

#### **1. Introduction**

22 Will-be-set-by-IN-TECH

194 Bioinformatics – Trends and Methodologies

Nuel, G. (2006a). Effective p-value computations using Finite Markov Chain Imbedding

Nuel, G. (2006b). Numerical solutions for patterns statistics on markov chains, *Stat. App. in*

Nuel, G. (2006c). Pattern statistics on markov chains and sensitivity to parameter estimation,

Nuel, G. (2008a). Pattern Markov chains: optimal Markov chain embedding through

Nuel, G. (2008b). Waiting time distribution for pattern occurrence in a constrained sequence:

Nuel, G. (2010). On the first k moments of the random count of a pattern in a multi-states sequence generated by a markov source, *Journal of Applied Probability* 47: 1–19. Nuel, G. & Dumas, J.-G. (2010). Sparse approaches for the exact distribution of patterns in long

Nuel, G. & Prum, B. (2007). Analyse statistique des séquences biologiques: modélisation

Nuel, G., Regad, L., Martin, J. & Camproux, A.-C. (2010). Exact distribution of a pattern in

Pevzner, P., Borodovski, M. & Mironov, A. (1989). Linguistic of nucleotide sequences:

Prum, B., Rodolphe, F. & de Turckheim, E. (1995). Finding words with unexpected frequencies

Reignier, M. (2000). A unified approach to word occurrences probabilities, *Discrete Applied*

Reinert, G. & Schbath, S. (1999). Compound poisson and poisson process approximations for occurrences of multiple words in markov chains, *J. of Comp. Biol.* 5: 223–254. Ribeca, P. & Raineri, E. (2008). Faster exact Markovian probability functions for motif occurrences: a DFA-only approach, *Bioinformatics* 24(24): 2839–2848. Roberts, R., Vincze, T., Posfai, J. & Macelis, D. (2010). REBASE – a database for dna restriction and modification: enzymes, genes and genomes, *Nucl. Acids Res.* 38: 234–236. Roquain, E. & Schbath, S. (2007). Improved compound poisson approximation for the number

Sigrist, C., Cerutti, L., de Castro, E., Langendijk-Genevaux, P., Bulliard, V., Bairoch, A. & Hulo,

Stefanov, V. T. & Szpankowski, W. (2007). Waiting Time Distributions for Pattern Occurrence

van Helden, J., André, B. & Collado-Vides, J. (1998). Extracting regulatory sites from

frequencies of occurrence of words, *J. Biomol. Struct. Dyn.* 6: 1013–1026. Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (1992). *Numerical Recipes in C*,

deterministic finite automata, *J. of Applied Prob.* 45(1): 226–243.

markovienne, alignements et motifs, *Hermes editions, Paris*.

*Biology* 1(1): 5.

*Science* 10: 3.

arXiv:1006.3246v1.

*Genet. and Mol. Biol.* 5(1): 26.

*Algorithms for Molecular Biology* 1(1): 17.

data, *Algorithms for Molecular Biology* 5: 15.

in dna sequences, *J. R. Statist. Soc. B* 11: 190–192.

Cambridge University Press.

*Mathematics* 104(1): 259–280.

*Probab.* 39(1): 128–140.

9(1): 305–320.

annotation, *Nucleic Acids Res.* 38.

frequencies, *J. Mol. Biol.* 281(5): 827–842.

(FMCI): application to local score and to pattern statistics, *Algorithms for Molecular*

an embedding markov chain approach, *Discrete Mathematics and Theoretical Computer*

multi-states sequences generated by a markov source, *submitted to J. Applied. Prob.* .

a set of random sequences generated by a markov source: applications to biological

The significance of deviation from mean statistical characteristics and prediction of

of occurrences of any rare word family in a stationary markov chain, *Adv. in Appl.*

N. (2010). PROSITE, a protein domain database for functional characterization and

in a Constrained Sequence, *Discrete Mathematics and Theoretical Computer Science*

the upstream region of yeast genes by computational analysis of oligonucleotide

The biological function of proteins is largely determined by their individual component domains, which are segments within the protein sequence that are self-contained and spatially arranged. These can be catalytic or structural, and define a number of different features of proteins such as their enzymatic activity, interactions with other proteins, sugars or lipids, and determine the cellular localization of the proteins that contain them. A number of intracellular three-dimensionally-arranged domains, such as Src-homology (SH) or Pleckstrin-homology (PH) domains, define the nature of protein interactions with other components of the cell, and enable them to interact with their substrates or binding partners. The specificity of interactions that is given by the domain is unique to its protein. Similarly, the extracellular part of most membrane-bound or secreted proteins of eukaryotic cells is also organized in semi-autonomously-arranged blocks that potentially confer multiple diverse functions to a particular protein. These domains have been classified and grouped into protein superfamilies depending on the similarity they have with domains of prototypical proteins, for example immunoglobulin, fibronectin or C-type lectin domains. Members of these groups are believed to be homologous and to have arisen by divergent evolution from a common ancestor. Many membrane-bound or extracellular proteins are comprised of several domains of the same type, but it is not uncommon to find mosaic proteins containing domains from different superfamilies.

The scavenger receptor cysteine-rich (SRCR) superfamily comprises a group of proteins that contain one or multiple domains structurally similar to the membrane distal domain of the type I scavenger receptor expressed by human macrophages (Freeman et al., 1990). Proteins classified as belonging to this superfamily may contain other types of domains additionally to the dominant SRCR modules, such as EGF, CUB, LCCL, or other domains. In mammals, SRCR proteins are typically expressed in cells of the immune system (Resnick et al., 1994), although some members can be also expressed in non-immune cells and organs, including liver, kidney, placenta, stomach, brain and heart (Sarrias et al., 2004). Group A domaincontaining SRCR proteins are present in phyla from the most primitive metazoan to

A Systematic and Thorough Search for Domains

may be the only unifying feature of the family.

**2.2 Structure and organization of SRCR domains**

2009; Padilla et al., 2002).

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 197

(Graversen et al., 2002; Kristiansen et al., 2001). The remaining three members of the group B SRCR family are secreted glycoproteins of different sizes and structural complexity. DMBT1, which was identified on the basis of its deletion in a medulloblastoma cell line, is the largest member of the family, comprising 14 SRCR domains separated by SRCR-interspacing domains (Mollenhauer et al., 1997). Apart from being secreted, DMBT1 is also found in association with the plasma membrane of macrophages, although it is not clear whether there is a specific receptor or the poorly characterized DMBT1 gene may encode a transmembrane sequence. Once in the membrane, DMBT1 is a ligand for Surfactant protein D (SP-D), a C-type lectin that binds to exposed carbohydrates (Holmskov et al., 1999). The SRCR soluble proteins S4D-SRCRB and SSc5D have four and five group B domains, respectively, and little is known of their functional or binding properties (Gonçalves et al.,

However, it has been recently suggested that Sp (Sarrias et al., 2005), DMBT1 (Bikker et al., 2002), CD163 (Fabriek et al., 2009), CD5 (Vera et al., 2009) and CD6 (Sarrias et al., 2007) are capable of detecting microbe-associated molecular patterns, and could bind and clear bacteria or fungi, reaffirming a scavenger-like role for this group of molecules. These developments notwithstanding, SRCR superfamily proteins may prove to have very diverse functions, to the extent that the structural properties of the highly conserved SRCR domains

Typically, the 100-110 amino acid-long SRCR domains possess a characteristic pattern of cysteine residues that establish intra-domain disulfide bridges and contribute to the overall architecture of the compact domain. The number of cysteine residues and their distribution, together with the organization of the genomic sequence encoding each domain, divide the SRCR family into two groups, A and B. Group A domains are encoded by split exons, and typically have six cysteine residues establishing three disulfide bonds. Group B domains, on the other hand, are encoded by a single exon and have eight cysteine residues, whose

So far, eight human SRCR group B proteins have been described, Sp, CD5, CD6, S4D, SSc5D, CD163, M160 and DMBT1 that contain three to fourteen SRCR domains. Their encoding genes are dispersed throughout the genome, however a few highly similar pairs such as CD5-CD6 and CD163-M160 are located on the same chromosome. The identity between individual domains of different SRCR group B proteins varies from 20 to 80%, and phylogenetic analysis suggests that they have evolved by sequential intragene duplication, although there are examples that suggest they may have evolved in some cases by interprotein domain shuffling. Only four SRCR domains have been characterized by X-ray crystallography, and of these three are group A SRCR domains, those of hepsin, a cell surface serine protease involved in cell growth and maintenance of cellular morphology (Somoza et al., 2003), M2bp, a tumor associated antigen and matrix protein (Hohenester et al., 1999), and MARCO, a trimeric SRCR group A protein expressed by macrophages and dendritic cells that recognizes polyanionic particles and pathogens (Ojala et al., 2007). The crystal structure of the membrane proximal domain of CD5 (Rodamilans et al., 2007), together with an NMR solution structure of domain 1 of CD5 (Garza-Garcia et al., 2008) constitute the only sources of structural information of SRCR group B domains. Comparing the structures, it is however apparent that the 3D assembly of the different domains in the

distribution is remarkably conserved in nearly all known domains (Fig. 1).

two groups is overall conserved, all displaying a very similar fold.

vertebrates, whereas group B domain containing SRCR proteins are only found in vertebrates. Intriguingly, although SRCR proteins can include other domains, no proteins have been reported to contain group A and B domains simultaneously.

In mammalian species, SRCR group B orthologs are usually very well conserved and regarding some of the proteins, a high level of conservation is extended to birds and fish. However, in some cases a human SRCR protein apparently has no corresponding ortholog in some mammals, and conversely, there are examples of SRCR group B proteins that are well characterized in a few mammalian species, that have not been described in humans. By analyzing the human genome, we can now identify all the remaining, still undescribed genes encoding SRCR group B domains, which will allow us to perform phylogenetic analysis of the complete set of group B domains. By comprehensive and systematic whole genome analysis we have found two new putative transcriptional units containing clusters of potential SRCR domains, and additionally a further putative gene that contains a single domain. After our thorough search, we are now confident that all proteins containing group B SRCR domains in the human genome have been identified.

#### **2. The scavenger receptor cysteine-rich group B family**

#### **2.1 Biological function of SRCR group B proteins**

The cell surface antigens CD5 and CD6, which function in T lymphocytes, are probably the most well characterized of the family, each containing three extracellular SRCR domains (Aruffo et al., 1991; Jones et al., 1986). CD5 and CD6 co-associate with each other at the surface of T cells (Castro et al., 2003; Gimferrer et al., 2003), and are involved in the regulation of T cell receptor-mediated activation. The extensive characterization of the interaction of CD6 with its ligand CD166, expressed by antigen presenting cells (Aruffo et al., 1997), and the identification of different binding partners for CD5 (Biancone et al., 1996; Calvo et al., 1999; Pospisil et al., 2000; Van de Velde et al., 1991), had initially suggested that SRCR group B domains participate in intercellular contacts *via* protein-protein interactions. Also, the three SRCR domain-containing soluble protein Sp (Gebe et al., 1997) has been reported to bind to cells of myeloid and lymphoid origin. Also known as AIM (apoptosis inhibitor expressed by macrophages), API6 (apoptosis inhibitor 6) or CD5L (CD5-like molecule), Sp is best known for promoting macrophage survival. Therefore, this sub-group of small SRCR-containing proteins may be described as having a role in cellular communication, differentiation and activation. However, for most of the remaining members of the family no such clear function has been established. In particular the lack of cellular ligands for most of these proteins raises the possibility that a totally different function for SRCR domains may exist, if indeed SRCR domain proteins share any common function.

In addition to CD5, CD6 and Sp, the group B SRCR family presently contains five other proteins, of which two, CD163 and M160, are membrane bound and expressed by macrophages. CD163 (Law et al., 1993) and M160 (CD163L1) (Gronlund et al., 2000), which were both identified in human monocytes, are considered a subgroup of the SRCR group B molecules. No definitive function has been established for these molecules, although CD163 has been described as binding to, and internalizing, tumor necrosis factor-like weak inducer of apoptosis (TWEAK), thus having a potential role in atherosclerosis (Moreno et al., 2009). Additionally, CD163 has a detoxifying role in iron metabolism, where by binding to hemoglobin-complexed haptoglobin it is able to remove hemoglobin from the plasma

vertebrates, whereas group B domain containing SRCR proteins are only found in vertebrates. Intriguingly, although SRCR proteins can include other domains, no proteins

In mammalian species, SRCR group B orthologs are usually very well conserved and regarding some of the proteins, a high level of conservation is extended to birds and fish. However, in some cases a human SRCR protein apparently has no corresponding ortholog in some mammals, and conversely, there are examples of SRCR group B proteins that are well characterized in a few mammalian species, that have not been described in humans. By analyzing the human genome, we can now identify all the remaining, still undescribed genes encoding SRCR group B domains, which will allow us to perform phylogenetic analysis of the complete set of group B domains. By comprehensive and systematic whole genome analysis we have found two new putative transcriptional units containing clusters of potential SRCR domains, and additionally a further putative gene that contains a single domain. After our thorough search, we are now confident that all proteins containing group

The cell surface antigens CD5 and CD6, which function in T lymphocytes, are probably the most well characterized of the family, each containing three extracellular SRCR domains (Aruffo et al., 1991; Jones et al., 1986). CD5 and CD6 co-associate with each other at the surface of T cells (Castro et al., 2003; Gimferrer et al., 2003), and are involved in the regulation of T cell receptor-mediated activation. The extensive characterization of the interaction of CD6 with its ligand CD166, expressed by antigen presenting cells (Aruffo et al., 1997), and the identification of different binding partners for CD5 (Biancone et al., 1996; Calvo et al., 1999; Pospisil et al., 2000; Van de Velde et al., 1991), had initially suggested that SRCR group B domains participate in intercellular contacts *via* protein-protein interactions. Also, the three SRCR domain-containing soluble protein Sp (Gebe et al., 1997) has been reported to bind to cells of myeloid and lymphoid origin. Also known as AIM (apoptosis inhibitor expressed by macrophages), API6 (apoptosis inhibitor 6) or CD5L (CD5-like molecule), Sp is best known for promoting macrophage survival. Therefore, this sub-group of small SRCR-containing proteins may be described as having a role in cellular communication, differentiation and activation. However, for most of the remaining members of the family no such clear function has been established. In particular the lack of cellular ligands for most of these proteins raises the possibility that a totally different function for SRCR domains may exist, if indeed SRCR domain proteins share any common

In addition to CD5, CD6 and Sp, the group B SRCR family presently contains five other proteins, of which two, CD163 and M160, are membrane bound and expressed by macrophages. CD163 (Law et al., 1993) and M160 (CD163L1) (Gronlund et al., 2000), which were both identified in human monocytes, are considered a subgroup of the SRCR group B molecules. No definitive function has been established for these molecules, although CD163 has been described as binding to, and internalizing, tumor necrosis factor-like weak inducer of apoptosis (TWEAK), thus having a potential role in atherosclerosis (Moreno et al., 2009). Additionally, CD163 has a detoxifying role in iron metabolism, where by binding to hemoglobin-complexed haptoglobin it is able to remove hemoglobin from the plasma

have been reported to contain group A and B domains simultaneously.

B SRCR domains in the human genome have been identified.

**2.1 Biological function of SRCR group B proteins**

function.

**2. The scavenger receptor cysteine-rich group B family**

(Graversen et al., 2002; Kristiansen et al., 2001). The remaining three members of the group B SRCR family are secreted glycoproteins of different sizes and structural complexity. DMBT1, which was identified on the basis of its deletion in a medulloblastoma cell line, is the largest member of the family, comprising 14 SRCR domains separated by SRCR-interspacing domains (Mollenhauer et al., 1997). Apart from being secreted, DMBT1 is also found in association with the plasma membrane of macrophages, although it is not clear whether there is a specific receptor or the poorly characterized DMBT1 gene may encode a transmembrane sequence. Once in the membrane, DMBT1 is a ligand for Surfactant protein D (SP-D), a C-type lectin that binds to exposed carbohydrates (Holmskov et al., 1999). The SRCR soluble proteins S4D-SRCRB and SSc5D have four and five group B domains, respectively, and little is known of their functional or binding properties (Gonçalves et al., 2009; Padilla et al., 2002).

However, it has been recently suggested that Sp (Sarrias et al., 2005), DMBT1 (Bikker et al., 2002), CD163 (Fabriek et al., 2009), CD5 (Vera et al., 2009) and CD6 (Sarrias et al., 2007) are capable of detecting microbe-associated molecular patterns, and could bind and clear bacteria or fungi, reaffirming a scavenger-like role for this group of molecules. These developments notwithstanding, SRCR superfamily proteins may prove to have very diverse functions, to the extent that the structural properties of the highly conserved SRCR domains may be the only unifying feature of the family.

#### **2.2 Structure and organization of SRCR domains**

Typically, the 100-110 amino acid-long SRCR domains possess a characteristic pattern of cysteine residues that establish intra-domain disulfide bridges and contribute to the overall architecture of the compact domain. The number of cysteine residues and their distribution, together with the organization of the genomic sequence encoding each domain, divide the SRCR family into two groups, A and B. Group A domains are encoded by split exons, and typically have six cysteine residues establishing three disulfide bonds. Group B domains, on the other hand, are encoded by a single exon and have eight cysteine residues, whose distribution is remarkably conserved in nearly all known domains (Fig. 1).

So far, eight human SRCR group B proteins have been described, Sp, CD5, CD6, S4D, SSc5D, CD163, M160 and DMBT1 that contain three to fourteen SRCR domains. Their encoding genes are dispersed throughout the genome, however a few highly similar pairs such as CD5-CD6 and CD163-M160 are located on the same chromosome. The identity between individual domains of different SRCR group B proteins varies from 20 to 80%, and phylogenetic analysis suggests that they have evolved by sequential intragene duplication, although there are examples that suggest they may have evolved in some cases by interprotein domain shuffling. Only four SRCR domains have been characterized by X-ray crystallography, and of these three are group A SRCR domains, those of hepsin, a cell surface serine protease involved in cell growth and maintenance of cellular morphology (Somoza et al., 2003), M2bp, a tumor associated antigen and matrix protein (Hohenester et al., 1999), and MARCO, a trimeric SRCR group A protein expressed by macrophages and dendritic cells that recognizes polyanionic particles and pathogens (Ojala et al., 2007). The crystal structure of the membrane proximal domain of CD5 (Rodamilans et al., 2007), together with an NMR solution structure of domain 1 of CD5 (Garza-Garcia et al., 2008) constitute the only sources of structural information of SRCR group B domains. Comparing the structures, it is however apparent that the 3D assembly of the different domains in the two groups is overall conserved, all displaying a very similar fold.

A Systematic and Thorough Search for Domains

**2.3 Homology between SRCR domains**

non-placental species (Table 1).

been reported in man.

SCART 1 and 2 molecules present in the mouse.

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 199

The level of amino acid identity among human SRCR group B domains from different molecules varies from 20% to 80%, but within the same molecule this level can be higher and even be identical in some domains (e.g. domains 3 and 7, and 10 and 11 of DMBT1). Similarly, some molecules are remarkably conserved between species, especially among mammals, although it appears that some level of conservation can be extended to birds, fish and amphibians in a few specific cases. There are good indications for there being orthologs of CD6 in the genomes of *T. guttata* and *D. rerio*, and some other examples. Nevertheless, the structure of SRCR group B-containing molecules is best preserved in mammalian species. The strong homology of SSc5D domains dates back to the divergence of egg- and non-egglaying mammals, while CD163 has clearly conserved orthologs in all mammals, including

Table 1. Homology between human and other mammalian CD163 domains. Numbers represent percentage of identity between each domain, compared to the human sequence. The significantly conserved homology of some SRCR orthologs is suggestive of profound functional constraints acting on these proteins. On the other hand, it appears that not all human SRCR B group molecules have described orthologs in all mammalian species, and conversely, that there are some SRCR proteins described in different animals that have not

Noticeably, bovine WC1 (Wijngaard et al., 1992) does not have a human counterpart, nor do the mouse SCART molecules (Kisielow et al., 2008). Similarly, the human macrophage specific receptor M160 is not found in all mammalian species, while the closely related molecule CD163, also specific to the monocytic/macrophage lineage, is clearly present in all genomes that we have examined. We have compared the similarity between individual domains of known and characterized members of the SRCR B group, and the corresponding domains in the bovine proteins (Table 2). While proteins such as S4D, SSc5D and CD163 show high levels of identity between human and bovine sequences, others like CD5 and Sp are more distantly related. M160 does not have a straightforward ortholog in cattle, so the bovine sequence used was of the related molecule M160-like, that is related in turn to the

Fig. 1. Sequence alignment of domains from group B SRCR superfamily members.

SRCR domains are typically sequences of 100-110 amino acids in length compacted into a heart-shaped fold, where a six/seven-stranded -sheet cradles an -helix. Strands 1, 3 and 4, together with 7, form a curved sheet that wraps around the core 1 helix. From 4 onwards, the structures start to diverge. It is the sequence of amino acids between the beginning of the domain and the 4 strand that is best conserved between group A and group B domains, and that roughly corresponds in the group A proteins to the first of two exons that encode a full SRCR A domain, and in group B proteins to the first 50 amino acids of the domain.

#### **2.3 Homology between SRCR domains**

198 Bioinformatics – Trends and Methodologies

Fig. 1. Sequence alignment of domains from group B SRCR superfamily members.

of the domain.

SRCR domains are typically sequences of 100-110 amino acids in length compacted into a heart-shaped fold, where a six/seven-stranded -sheet cradles an -helix. Strands 1, 3 and 4, together with 7, form a curved sheet that wraps around the core 1 helix. From 4 onwards, the structures start to diverge. It is the sequence of amino acids between the beginning of the domain and the 4 strand that is best conserved between group A and group B domains, and that roughly corresponds in the group A proteins to the first of two exons that encode a full SRCR A domain, and in group B proteins to the first 50 amino acids The level of amino acid identity among human SRCR group B domains from different molecules varies from 20% to 80%, but within the same molecule this level can be higher and even be identical in some domains (e.g. domains 3 and 7, and 10 and 11 of DMBT1). Similarly, some molecules are remarkably conserved between species, especially among mammals, although it appears that some level of conservation can be extended to birds, fish and amphibians in a few specific cases. There are good indications for there being orthologs of CD6 in the genomes of *T. guttata* and *D. rerio*, and some other examples. Nevertheless, the structure of SRCR group B-containing molecules is best preserved in mammalian species. The strong homology of SSc5D domains dates back to the divergence of egg- and non-egglaying mammals, while CD163 has clearly conserved orthologs in all mammals, including non-placental species (Table 1).

Table 1. Homology between human and other mammalian CD163 domains. Numbers represent percentage of identity between each domain, compared to the human sequence.

The significantly conserved homology of some SRCR orthologs is suggestive of profound functional constraints acting on these proteins. On the other hand, it appears that not all human SRCR B group molecules have described orthologs in all mammalian species, and conversely, that there are some SRCR proteins described in different animals that have not been reported in man.

Noticeably, bovine WC1 (Wijngaard et al., 1992) does not have a human counterpart, nor do the mouse SCART molecules (Kisielow et al., 2008). Similarly, the human macrophage specific receptor M160 is not found in all mammalian species, while the closely related molecule CD163, also specific to the monocytic/macrophage lineage, is clearly present in all genomes that we have examined. We have compared the similarity between individual domains of known and characterized members of the SRCR B group, and the corresponding domains in the bovine proteins (Table 2). While proteins such as S4D, SSc5D and CD163 show high levels of identity between human and bovine sequences, others like CD5 and Sp are more distantly related. M160 does not have a straightforward ortholog in cattle, so the bovine sequence used was of the related molecule M160-like, that is related in turn to the SCART 1 and 2 molecules present in the mouse.

A Systematic and Thorough Search for Domains

achieved. Clearly, either hypothesis does not exclude the other.

comprises five SRCR group B domains (Gonçalves et al., 2009).

novel proteins (Table 3).

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 201

SCART (CD163c), of which there are two genes found in the mouse, and the bovine gene M160-like (Herzig et al., 2010). These sets of genes are related to WC1 genes expressed in cattle, sheep and swine. To obtain a better idea of the relationship between these families of genes, we aligned the full sequences of known human and bovine proteins containing SRCR group B domains using ClustalW and drew the corresponding phylogram (Fig. 2). As can be seen, there are no direct links between human M160, bovine M160L and bovine WC1, raising the possibility that either some genes were lost during mammalian evolution, or that the complete characterization and annotation of the genomes has still not been fully

**3. A systematic and thorough search for SRCR domains in the genome**

Our hope is that the evolution and function of SRCR domains would emerge when all members of this protein family have been identified. The advent of the human genome sequence has allowed us to screen, using bioinformatics-based approaches, for new SRCR proteins still not described or characterized. We decided to focus on group B molecules, given that proteins of this type are more conserved, restricted in number, and their specialized function, in this case immune-related, seems better defined. We performed searches for new members of the SRCR-SF in the completed human genome sequence by interrogating the genome using TBLASTN 2.2.20+ (Altschul et al., 1997; http://www.ncbi.nlm.nih.gov/BLAST). Initially, we screened for new sequences exhibiting similarity with any or all of the SRCR domains comprising the then known SRCR superfamily proteins (Gonçalves et al., 2009). We expected that, for a given TBLASTN run, *bona fide* new SRCR domains would have smaller E values than the best matches of the search sequence with Group A SRCR domains. According to this criterion, the search identified the sequences encoding domains within already known and characterized proteins *i.e.* CD5, CD6, Sp, S4D, CD163, M160 and DMBT1. Additionally, we identified a cluster of five new SRCR domains, which we further investigated and that later resulted in the cloning and characterization of SSc5D, a molecule secreted by macrophages and that

A *caveat* in our methodology was that not all group B domains were identified using this strategy. The most divergent domains, namely those of CD5, were not retrieved in all searches, and in particular CD5-d1 was rarely identified as having a clear homology with any other group B domain. Sequence alignment of all group B domains (Fig. 1) highlights the striking differences of CD5 sequences, and also to some extent of the CD6 domains, when compared with other sequences that are remarkably similar to each other. In order to perform a more rigorous search for all putative SRCR group B domains, we conducted a comprehensive systematic search using PSI-BLAST (Altschul et al., 1997) to find distant homologs. All known SRCR domains were used as queries to search iteratively against human non-redundant database with the target sequence length set to 250. BLOSUM62 amino acid substitution matrix with gap open penalty 11, and extension penalty 1, was used. Sequence masking was disabled and the PSI-BLAST threshold was set to 0.005. While searching, each PSI-BLAST query was iterated including new hits from the previous search until it converged, i.e. no new hits were found in subsequent searches. After each search iteration, results were checked for new SRCR proteins. This meticulous and robust method picked up all known SRCR domain-containing proteins along with

Table 2. Similarity between human and bovine corresponding SRCR domains. Percentage identity between each domain is indicated.

Fig. 2. Relationships between human and bovine SRCR molecules.

There are three groups of genes in the CD163 family: CD163 itself (CD163a), present in all mammals; M160 (CD163b), so far only found in the genomes of primates and in horses; and

Table 2. Similarity between human and bovine corresponding SRCR domains. Percentage

There are three groups of genes in the CD163 family: CD163 itself (CD163a), present in all mammals; M160 (CD163b), so far only found in the genomes of primates and in horses; and

identity between each domain is indicated.

Fig. 2. Relationships between human and bovine SRCR molecules.

SCART (CD163c), of which there are two genes found in the mouse, and the bovine gene M160-like (Herzig et al., 2010). These sets of genes are related to WC1 genes expressed in cattle, sheep and swine. To obtain a better idea of the relationship between these families of genes, we aligned the full sequences of known human and bovine proteins containing SRCR group B domains using ClustalW and drew the corresponding phylogram (Fig. 2). As can be seen, there are no direct links between human M160, bovine M160L and bovine WC1, raising the possibility that either some genes were lost during mammalian evolution, or that the complete characterization and annotation of the genomes has still not been fully achieved. Clearly, either hypothesis does not exclude the other.

#### **3. A systematic and thorough search for SRCR domains in the genome**

Our hope is that the evolution and function of SRCR domains would emerge when all members of this protein family have been identified. The advent of the human genome sequence has allowed us to screen, using bioinformatics-based approaches, for new SRCR proteins still not described or characterized. We decided to focus on group B molecules, given that proteins of this type are more conserved, restricted in number, and their specialized function, in this case immune-related, seems better defined. We performed searches for new members of the SRCR-SF in the completed human genome sequence by interrogating the genome using TBLASTN 2.2.20+ (Altschul et al., 1997; http://www.ncbi.nlm.nih.gov/BLAST). Initially, we screened for new sequences exhibiting similarity with any or all of the SRCR domains comprising the then known SRCR superfamily proteins (Gonçalves et al., 2009). We expected that, for a given TBLASTN run, *bona fide* new SRCR domains would have smaller E values than the best matches of the search sequence with Group A SRCR domains. According to this criterion, the search identified the sequences encoding domains within already known and characterized proteins *i.e.* CD5, CD6, Sp, S4D, CD163, M160 and DMBT1. Additionally, we identified a cluster of five new SRCR domains, which we further investigated and that later resulted in the cloning and characterization of SSc5D, a molecule secreted by macrophages and that comprises five SRCR group B domains (Gonçalves et al., 2009).

A *caveat* in our methodology was that not all group B domains were identified using this strategy. The most divergent domains, namely those of CD5, were not retrieved in all searches, and in particular CD5-d1 was rarely identified as having a clear homology with any other group B domain. Sequence alignment of all group B domains (Fig. 1) highlights the striking differences of CD5 sequences, and also to some extent of the CD6 domains, when compared with other sequences that are remarkably similar to each other. In order to perform a more rigorous search for all putative SRCR group B domains, we conducted a comprehensive systematic search using PSI-BLAST (Altschul et al., 1997) to find distant homologs. All known SRCR domains were used as queries to search iteratively against human non-redundant database with the target sequence length set to 250. BLOSUM62 amino acid substitution matrix with gap open penalty 11, and extension penalty 1, was used. Sequence masking was disabled and the PSI-BLAST threshold was set to 0.005. While searching, each PSI-BLAST query was iterated including new hits from the previous search until it converged, i.e. no new hits were found in subsequent searches. After each search iteration, results were checked for new SRCR proteins. This meticulous and robust method picked up all known SRCR domain-containing proteins along with novel proteins (Table 3).

A Systematic and Thorough Search for Domains

Fig. 3. Multiple sequence alignment of all SRCR domains.

family (Fig. 4).

The overall height of the stack indicates the sequence conservation at that position, while the height of symbols within the stack indicates the relative frequency of each amino acid at that position. It is apparent from the WebLogo that, although sequences vary substantially between SRCR domains, all cysteine residues (colored in red) are conserved across the

Fig. 4. Multiple sequence alignment of SRCR domains in WebLogo format.

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 203

Table 3. List of SRCR-containing proteins. \* - DMBT1 has been described as containing 14 SRCR domains; # - annotated as a pseudogene; in red denotes new SRCR domains from uncharacterized proteins.

From this genome-wide search we obtained a total of 76 SRCR group B domains distributed in 11 genes, each putatively encoding a varying number of SRCR domains. The eleven genes are spread across the genome on seven different chromosomes: chromosome 10 contains three SRCR group B-encoding genes, chromosomes 11 and 12 contain two each, and chromosomes 1, 7, 14 and 19 each contain a copy of just one SRCR-encoding gene. Among these genes and in addition to known domains from characterized genes, our search has uncovered 23 new putative group B domains, three of which represent previously unreported domains localized within the DMBT1 gene. The DMBT1 gene thus putatively encodes a maximum of 17 SRCR domains. Some controversy has existed on the number of SRCR domains within the DMBT1 molecule. Like other SRCR-containing proteins (Castro et al., 2007; Padilla et al., 2002), DMBT1 can be expressed as different isoforms arising by alternative splicing, which include or exclude individual SRCR domains (Mollenhauer et al., 1999). It is possible that the new DMBT1 domains have not been previously reported because they are not expressed in the tissues or cells investigated, however it is also plausible that the exons coding for these domains have been silenced during evolution and are now non-functional.

The remaining new domains belong to 3 new putative genes, one 8 domain-encoding gene (8D), one gene, annotated as a DMBT1-like pseudo-gene, that encodes 11 fragments of SRCR domains of variable lengths (D11), and a gene encoding a putative Hedgehog interacting protein-like 1 molecule (HHIP-like 1), which contains a single SRCR domain.

In order to analyze the sequence conservation and diversity of SRCR domains, we aligned all individual 76 SRCR group B domains using ClustalW2 (Thompson et al., 1994) with the default substitution matrix (Gonnet series) and gap opening and extension penalties of 10 and 0.2 respectively. Due to the sequence diversity in SRCR domains, several insertions and deletions were found in the multiple sequence alignment (Fig. 3).

To locate sequence patterns as well as conserved amino acids, the multiple sequence alignment was used to create a WebLogo (Crooks et al., 2004; http://weblogo.berkeley.edu).

Protein Number of SRCR Domains Chromosome

Sp 3 3 S4D 4 7 DMBT1 17\* 10 CD5 3 11 CD6 3 11 CD163 9 12 M160 12 12 SSc5D 5 19 8D 8 10 D11# 11 10 HHIPL1 1 14 Table 3. List of SRCR-containing proteins. \* - DMBT1 has been described as containing 14 SRCR domains; # - annotated as a pseudogene; in red denotes new SRCR domains from

From this genome-wide search we obtained a total of 76 SRCR group B domains distributed in 11 genes, each putatively encoding a varying number of SRCR domains. The eleven genes are spread across the genome on seven different chromosomes: chromosome 10 contains three SRCR group B-encoding genes, chromosomes 11 and 12 contain two each, and chromosomes 1, 7, 14 and 19 each contain a copy of just one SRCR-encoding gene. Among these genes and in addition to known domains from characterized genes, our search has uncovered 23 new putative group B domains, three of which represent previously unreported domains localized within the DMBT1 gene. The DMBT1 gene thus putatively encodes a maximum of 17 SRCR domains. Some controversy has existed on the number of SRCR domains within the DMBT1 molecule. Like other SRCR-containing proteins (Castro et al., 2007; Padilla et al., 2002), DMBT1 can be expressed as different isoforms arising by alternative splicing, which include or exclude individual SRCR domains (Mollenhauer et al., 1999). It is possible that the new DMBT1 domains have not been previously reported because they are not expressed in the tissues or cells investigated, however it is also plausible that the exons coding for these domains have been silenced during evolution and

The remaining new domains belong to 3 new putative genes, one 8 domain-encoding gene (8D), one gene, annotated as a DMBT1-like pseudo-gene, that encodes 11 fragments of SRCR domains of variable lengths (D11), and a gene encoding a putative Hedgehog interacting

In order to analyze the sequence conservation and diversity of SRCR domains, we aligned all individual 76 SRCR group B domains using ClustalW2 (Thompson et al., 1994) with the default substitution matrix (Gonnet series) and gap opening and extension penalties of 10 and 0.2 respectively. Due to the sequence diversity in SRCR domains, several insertions and

To locate sequence patterns as well as conserved amino acids, the multiple sequence alignment was used to create a WebLogo (Crooks et al., 2004; http://weblogo.berkeley.edu).

protein-like 1 molecule (HHIP-like 1), which contains a single SRCR domain.

deletions were found in the multiple sequence alignment (Fig. 3).

uncharacterized proteins.

are now non-functional.

Fig. 3. Multiple sequence alignment of all SRCR domains.

The overall height of the stack indicates the sequence conservation at that position, while the height of symbols within the stack indicates the relative frequency of each amino acid at that position. It is apparent from the WebLogo that, although sequences vary substantially between SRCR domains, all cysteine residues (colored in red) are conserved across the family (Fig. 4).

Fig. 4. Multiple sequence alignment of SRCR domains in WebLogo format.

A Systematic and Thorough Search for Domains

Carmo, unpublished).

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 205

Fig. 5. Maximum likelihood phylogenetic estimation of human SRCR domains. Internal branch reliability was assessed using bootstrapping method (100 bootstrap replicates). Branches observed in more than 75 of 100 bootstrapped re-sampling are shown in red.

The second cluster of SRCR domains we have uncovered is located on chromosome 10 and includes 8 such domains, thus we provisionally termed it 8D. Using the predicted exonderived protein sequence, we BLAST-searched other mammalian genomes and the proteins that we retrieved which were most similar to human 8D were mouse SCART1, bovine M160L, and mouse SCART2, whose ClustalW alignment scores were 70, 66 and 57, respectively (Fig. 6). Human 8D and mouse SCART1 have 64% identity for the entire sequence, while some individual SRCR domains share identities of close to, and even above 80%. We thus believe that 8D is the human ortholog of mouse SCART1. It remains to be seen whether human 8D can be expressed and produce a mature and functional protein, although we have detected several 8D transcripts of different sizes (C Gonçalves and A

In order to estimate evolutionary relationships among SRCR domains, we performed a phylogenetic analysis utilizing the maximum likelihood phylogenetic reconstruction method (proml) available in the phylip phylogenetic package (Felsenstein, 1989). Gaps were trimmed from the SRCR multiple sequence alignment prior to tree building and the Jones-Taylor-Thornton probability model employed with constant rate of change among sites. The reliability of internal branches was subsequently evaluated using 100 bootstrap samplings. SRCR domains exhibit very complex evolutionary relationships (Fig. 5). In the reconstructed phylogenetic tree, intra-protein domain clustering as well as inter-protein domain clusters were observed. Intra-domain clustering, as in the case of DMBT1, strongly suggests the evolution of these domains via sequential intragenic duplication. At the same time it is difficult to understand the inter-protein domain similarities. CD5, CD6, Spα, M160 and CD163 exhibit more diverse relationships. Among them, SRCR domains show greater inter than intra protein similarities. Given their low sequence similarities, it is uncertain whether the domains evolved through gene duplication and accumulated mutations have reduced sequence similarity, or if it is through a convergent evolution mechanism subsequent to domain shuffling. The similarities of domain pairs M160\_d4-CD163\_d1, M160\_d7-CD163\_d4 M160\_d8-CD163\_d5, M160\_d9-CD163\_d6, M160\_d10-CD163\_d7, M160\_d11-CD163\_d8, M160\_d12–CD163\_d9, CD5\_d1-CD6\_d3 and CD5\_d2-S4D\_d4 are strongly suggestive of inter-protein domain shuffling.

#### **4. Concluding remarks - the completion of the SRCR group B family**

In contrast to the complexity and variety of large protein families such as the G protein coupled receptor (GPCR) superfamily, which has nearly 800 genes in the human genome, corresponding to roughly 4% of the full protein-encoding genome, group B of the scavenger receptor cysteine-rich superfamily appears to be much more limited. So far it includes only 8 members in the entire human genome, although there are additionally 3 proteins described in other mammalians; SCART1 and SCART2 initially found in mice, and the 11 SRCR domain-containing protein WC1 expressed in cattle, sheep and swine. Also, the function of mammalian SRCR proteins seems to be restricted to the immune system, although the exact nature or biological role of the family is still to be fully determined.

In this study we set out to identify the remaining members of the SRCR group B family in order to obtain a clear understanding of the biological significance of this important group of proteins and to clarify some as yet unresolved questions regarding their evolution in mammalian species. We searched the human genome for the presence of SRCR-encoding genes using as probes the amino acid sequence of all reported human SRCR domains. Interestingly, one of the new members we have identified, HHIP-like 1, contains a single SRCR domain, which is unknown in the family. Moreover, the amino acid sequence corresponding to the SRCR domain constitutes only a small fraction of the total of the putative protein (13%). This is in contrast with most other members, whose amino acid content corresponding to SRCR domains relative to the whole of the protein is significantly higher, varying between 32% (SSc5D) and 92% (Sp). HHIP-like 1 is related to Hedgehog interacting protein, a regulatory component of the Hedgehog signaling pathway (Chuang and McMahon, 1999). However, unlike HHIP-like 1, HHIP does not contain an SRCR domain.

In order to estimate evolutionary relationships among SRCR domains, we performed a phylogenetic analysis utilizing the maximum likelihood phylogenetic reconstruction method (proml) available in the phylip phylogenetic package (Felsenstein, 1989). Gaps were trimmed from the SRCR multiple sequence alignment prior to tree building and the Jones-Taylor-Thornton probability model employed with constant rate of change among sites. The reliability of internal branches was subsequently evaluated using 100 bootstrap samplings. SRCR domains exhibit very complex evolutionary relationships (Fig. 5). In the reconstructed phylogenetic tree, intra-protein domain clustering as well as inter-protein domain clusters were observed. Intra-domain clustering, as in the case of DMBT1, strongly suggests the evolution of these domains via sequential intragenic duplication. At the same time it is difficult to understand the inter-protein domain similarities. CD5, CD6, Spα, M160 and CD163 exhibit more diverse relationships. Among them, SRCR domains show greater inter than intra protein similarities. Given their low sequence similarities, it is uncertain whether the domains evolved through gene duplication and accumulated mutations have reduced sequence similarity, or if it is through a convergent evolution mechanism subsequent to domain shuffling. The similarities of domain pairs M160\_d4-CD163\_d1, M160\_d7-CD163\_d4 M160\_d8-CD163\_d5, M160\_d9-CD163\_d6, M160\_d10-CD163\_d7, M160\_d11-CD163\_d8, M160\_d12–CD163\_d9, CD5\_d1-CD6\_d3 and CD5\_d2-S4D\_d4 are strongly suggestive of

**4. Concluding remarks - the completion of the SRCR group B family**

nature or biological role of the family is still to be fully determined.

In contrast to the complexity and variety of large protein families such as the G protein coupled receptor (GPCR) superfamily, which has nearly 800 genes in the human genome, corresponding to roughly 4% of the full protein-encoding genome, group B of the scavenger receptor cysteine-rich superfamily appears to be much more limited. So far it includes only 8 members in the entire human genome, although there are additionally 3 proteins described in other mammalians; SCART1 and SCART2 initially found in mice, and the 11 SRCR domain-containing protein WC1 expressed in cattle, sheep and swine. Also, the function of mammalian SRCR proteins seems to be restricted to the immune system, although the exact

In this study we set out to identify the remaining members of the SRCR group B family in order to obtain a clear understanding of the biological significance of this important group of proteins and to clarify some as yet unresolved questions regarding their evolution in mammalian species. We searched the human genome for the presence of SRCR-encoding genes using as probes the amino acid sequence of all reported human SRCR domains. Interestingly, one of the new members we have identified, HHIP-like 1, contains a single SRCR domain, which is unknown in the family. Moreover, the amino acid sequence corresponding to the SRCR domain constitutes only a small fraction of the total of the putative protein (13%). This is in contrast with most other members, whose amino acid content corresponding to SRCR domains relative to the whole of the protein is significantly higher, varying between 32% (SSc5D) and 92% (Sp). HHIP-like 1 is related to Hedgehog interacting protein, a regulatory component of the Hedgehog signaling pathway (Chuang and McMahon, 1999). However, unlike HHIP-like 1, HHIP does not contain an SRCR

inter-protein domain shuffling.

domain.

Fig. 5. Maximum likelihood phylogenetic estimation of human SRCR domains. Internal branch reliability was assessed using bootstrapping method (100 bootstrap replicates). Branches observed in more than 75 of 100 bootstrapped re-sampling are shown in red.

The second cluster of SRCR domains we have uncovered is located on chromosome 10 and includes 8 such domains, thus we provisionally termed it 8D. Using the predicted exonderived protein sequence, we BLAST-searched other mammalian genomes and the proteins that we retrieved which were most similar to human 8D were mouse SCART1, bovine M160L, and mouse SCART2, whose ClustalW alignment scores were 70, 66 and 57, respectively (Fig. 6). Human 8D and mouse SCART1 have 64% identity for the entire sequence, while some individual SRCR domains share identities of close to, and even above 80%. We thus believe that 8D is the human ortholog of mouse SCART1. It remains to be seen whether human 8D can be expressed and produce a mature and functional protein, although we have detected several 8D transcripts of different sizes (C Gonçalves and A Carmo, unpublished).

A Systematic and Thorough Search for Domains

Vattipally Sreenu is funded by the University of Oxford.

lymphocytes. J Exp Med *184*, 811-819.

surface ligand. Eur J Immunol *29*, 2119-2129.

synapse. Journal of Immunology *178*, 4351-4361.

generator. Genome Res *14*, 1188-1190.

search programs. Nucleic Acids Res *25*, 3389-3402.

surface and secreted proteins. J Exp Med *174*, 949-952.

**5. Acknowledgements**

Today *18*, 498-504.

32109-32115.

Biology *73*, 183-190.

*113*, 887-892.

Science *246*, 941-942.

**6. References**

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 207

We thank Dr. Simon Lee for reviewing this manuscript. The work in Alexandre Carmo's laboratory is funded by FEDER through the Programa Operacional Factores de Competitividade – COMPETE, and by FCT – Fundação para a Ciência e a Tecnologia.

Altschul, S.F., Madden, T.L., Schäffer, A.A., Zhang, J., Zhang, Z., Miller, W., and Lipman,

Aruffo, A., Bowen, M.A., Patel, D.D., Haynes, B.F., Starling, G.C., Gebe, J.A., and Bajorath, J.

Aruffo, A., Melnick, M.B., Linsley, P.S., and Seed, B. (1991). The lymphocyte glycoprotein

Biancone, L., Bowen, M.A., Lim, A., Aruffo, A., Andres, G., and Stamenkovic, I. (1996).

Bikker, F.J., Ligtenberg, A.J., Nazmi, K., Veerman, E.C., van't Hof, W., Bolscher, J.G.,

Calvo, J., Places, L., Padilla, O., Vilà, J.M., Vives, J., Bowen, M.A., and Lozano, F. (1999).

Castro, M.A.A., Nunes, R.J., Oliveira, M.I., Tavares, P.A., Simões, C., Parnes, J.R., Moreira,

Castro, M.A.A., Oliveira, M.I., Nunes, R.J., Fabre, S., Barbosa, R., Peixoto, A., Brown, M.H.,

Chuang, P.T., and McMahon, A.P. (1999). Vertebrate Hedgehog signalling modulated by

Crooks, G.E., Hon, G., Chandonia, J.M., and Brenner, S.E. (2004). WebLogo: a sequence logo

Fabriek, B.O., van Bruggen, R., Deng, D.M., Ligtenberg, A.J., Nazmi, K., Schornagel, K.,

Felsenstein, J. (1989). Mathematics vs. Evolution: Mathematical Evolutionary Theory.

induction of a Hedgehog-binding protein. Nature *397*, 617-621.

D.J. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database

(1997). CD6-ligand interactions: a paradigm for SRCR domain function? Immunol

CD6 contains a repeated domain structure characteristic of a new family of cell

Identification of a novel inducible cell-surface ligand of CD5 on activated

Poustka, A., Nieuw Amerongen, A.V., and Mollenhauer, J. (2002). Identification of the bacteria-binding peptide domain on salivary agglutinin (gp-340/DMBT1), a member of the scavenger receptor cysteine-rich superfamily. J Biol Chem *277*,

Interaction of recombinant and natural soluble CD5 forms with an alternative cell

A., and Carmo, A.M. (2003). OX52 is the rat homologue of CD6: evidence for an effector function in the regulation of CD5 phosphorylation. Journal of Leukocyte

Parnes, J.R., Bismuth, G., Moreira, A.*, et al.* (2007). Extracellular Isoforms of CD6 generated by alternative splicing regulate targeting of CD6 to the immunological

Vloet, R.P., Dijkstra, C.D., and van den Berg, T.K. (2009). The macrophage scavenger receptor CD163 functions as an innate immune sensor for bacteria. Blood

Fig. 6. Sequence alignment of mouse SCART1, human 8D and bovine M160-like protein.

The last new set of domains we have identified is located in a gene also on chromosome 10, but has been annotated as a non-coding pseudogene. Analysis of its putative sequence derived from the exon-like sequences in fact reveal that some stretches of several of the SRCR domains are missing, adding to a number of frameshifts and premature stop codons. Curiously, 11 SRCR-like domains can be identified, exactly the same number as the typical bovine WC1 protein. Comparison between the two sequences has failed however to definitely determine whether these two genes have the same evolutionary origin, as individually identifiable SRCR or SRCR-type domains seem to have already drifted apart significantly.

With the recognition of the three new genes, albeit none of them proven to be functional as yet together with the detection of three new putative SRCR-encoding sequences present in the DMBT1 gene, we are confident that we have completed the identification of the full set of scavenger receptor cysteine-rich group B domains in the human genome.

#### **5. Acknowledgements**

We thank Dr. Simon Lee for reviewing this manuscript. The work in Alexandre Carmo's laboratory is funded by FEDER through the Programa Operacional Factores de Competitividade – COMPETE, and by FCT – Fundação para a Ciência e a Tecnologia. Vattipally Sreenu is funded by the University of Oxford.

#### **6. References**

206 Bioinformatics – Trends and Methodologies

Fig. 6. Sequence alignment of mouse SCART1, human 8D and bovine M160-like protein.

significantly.

The last new set of domains we have identified is located in a gene also on chromosome 10, but has been annotated as a non-coding pseudogene. Analysis of its putative sequence derived from the exon-like sequences in fact reveal that some stretches of several of the SRCR domains are missing, adding to a number of frameshifts and premature stop codons. Curiously, 11 SRCR-like domains can be identified, exactly the same number as the typical bovine WC1 protein. Comparison between the two sequences has failed however to definitely determine whether these two genes have the same evolutionary origin, as individually identifiable SRCR or SRCR-type domains seem to have already drifted apart

With the recognition of the three new genes, albeit none of them proven to be functional as yet together with the detection of three new putative SRCR-encoding sequences present in the DMBT1 gene, we are confident that we have completed the identification of the full set

of scavenger receptor cysteine-rich group B domains in the human genome.

A Systematic and Thorough Search for Domains

Nat Genet *17*, 32-39.

Biol Chem *282*, 16654-16666.

Lymphoma *36*, 353-365.

USA *104*, 11724-11729.

37.

domain. J Biol Chem *282*, 12669-12677.

recognition receptor. J Biol Chem *280*, 35391-35398.

atherosclerosis. Atherosclerosis *207*, 103-110.

6240.

of the Scavenger Receptor Cysteine-Rich Group-B Family in the Human Genome 209

Mollenhauer, J., Holmskov, U., Wiemann, S., Krebs, I., Herbertz, S., Madsen, J., Kioschis, P.,

Mollenhauer, J., Wiemann, S., Scheurlen, W., Korn, B., Hayashi, Y., Wilgenbus, K.K., von

Moreno, J.A., Muñoz-García, B., Martín-Ventura, J.L., Madrigal-Matute, J., Orbe, J., Páramo,

Ojala, J.R., Pikkarainen, T., Tuuttila, A., Sandalova, T., and Tryggvason, K. (2007). Crystal

Padilla, O., Pujana, M.A., López-de la Iglesia, A., Gimferrer, I., Arman, M., Vilà, J.M., Places,

Resnick, D., Pearson, A., and Krieger, M. (1994). The SRCR superfamily: a family

Rodamilans, B., Muñoz, I.G., Bragado-Nilsson, E., Sarrias, M.R., Padilla, O., Blanco, F.J.,

Sarrias, M.R., Roselló, S., Sánchez-Barbero, F., Sierra, J.M., Vila, J., Yélamos, J., Vives, J.,

Sarrias, M.R., Farnós, M., Mota, R., Sánchez-Barbero, F., Ibáñez, A., Gimferrer, I., Vera, J.,

Sarrias, M.R., Grønlund, J., Padilla, O., Madsen, J., Holmskov, U., and Lozano, F. (2004). The

Somoza, J.R., Ho, J.D., Luong, C., Ghate, M., Sprengeler, P.A., Mortara, K., Shrader, W.D.,

scavenger receptor cysteine-rich (SRCR) domain. Structure *11*, 1123-1131.

reminiscent of the Ig superfamily. Trends Biochem Sci *19*, 5-8.

mapping to human chromosome 7q11.23. Immunogenetics *54*, 621-634. Pospisil, R., Silverman, G.J., Marti, G.E., Aruffo, A., Bowen, M.A., and Mage, R.G. (2000).

Coy, J.F., and Poustka, A. (1999). The genomic structure of the DMBT1 gene: evidence for a region with susceptibility to genomic instability. Oncogene *18*, 6233-

Deimling, A., and Poustka, A. (1997). DMBT1, a new member of the SRCR superfamily, on chromosome 10q25.3-26.1 is deleted in malignant brain tumours.

J.A., Ortega, L., Egido, J., and Blanco-Colio, L.M. (2009). The CD163-expressing macrophages recognize and internalize TWEAK: potential consequences in

structure of the cysteine-rich domain of scavenger receptor MARCO reveals the presence of a basic and an acidic cluster that both contribute to ligand recognition. J

L., Vives, J., Estivill, X., and Lozano, F. (2002). Cloning of S4D-SRCRB, a new soluble member of the group B scavenger receptor cysteine-rich family (SRCR-SF)

CD5 is a potential selecting ligand for B-cell surface immunoglobulin: a possible role in maintenance and selective expansion of normal and malignant B cells. Leuk

Lozano, F., and Montoya, G. (2007). Crystal structure of the third extracellular domain of CD5 reveals the fold of a group B scavenger cysteine-rich receptor

Casals, C., and Lozano, F. (2005). A role for human Sp alpha as a pattern

Fenutría, R., Casals, C., Yélamos, J.*, et al.* (2007). CD6 binds to pathogen-associated molecular patterns and protects from LPS-induced septic shock. Proc Natl Acad Sci

Scavenger Receptor Cysteine-Rich (SRCR) domain: an ancient and highly conserved protein module of the innate immune system. Crit Rev Immunol *24*, 1-

Sperandio, D., Chan, H., McGrath, M.E.*, et al.* (2003). The structure of the extracellular region of human hepsin reveals a serine protease domain and a novel

Freeman, M., Ashkenas, J., Rees, D.J., Kingsley, D.M., Copeland, N.G., Jenkins, N.A., and

Garza-Garcia, A., Esposito, D., Rieping, W., Harris, R., Briggs, C., Brown, M.H., and Driscoll,

Gebe, J.A., Kiener, P.A., Ring, H.Z., Li, X., Francke, U., and Aruffo, A. (1997). Molecular

Gimferrer, I., Farnós, M., Calvo, M., Mittelbrunn, M., Enrich, C., Sánchez-Madrid, F., Vives,

Gonçalves, C.M., Castro, M.A.A., M., Henriques, T., Oliveira, M.I., Pinheiro, H.C., H.,

Graversen, J.H., Madsen, M., and Moestrup, S.K. (2002). CD163: a signal receptor scavenging

Gronlund, J., Vitved, L., Lausen, M., Skjodt, K., and Holmskov, U. (2000). Cloning of a novel

Herzig, C.T., Waters, R.W., Baldwin, C.L., and Telfer, J.C. (2010). Evolution of the CD163

Hohenester, E., Sasaki, T., and Timpl, R. (1999). Crystal structure of a scavenger receptor

Holmskov, U., Mollenhauer, J., Madsen, J., Vitved, L., Gronlund, J., Tornoe, I., Kliem, A.,

Kisielow, J., Kopf, M., and Karjalainen, K. (2008). SCART scavenger receptors identify a novel subset of adult gammadelta T cells. J Immunol *181*, 1710-1716. Kristiansen, M., Graversen, J.H., Jacobsen, C., Sonne, O., Hoffman, H.J., Law, S.K., and

Law, S.K., Micklem, K.J., Shaw, J.M., Zhang, X.P., Dong, Y., Willis, A.C., and Mason, D.Y.

membrane of lymphoid T cells. J Biol Chem *278*, 8564-8571.

rich superfamily. Molecular Immunology *46*, 2585-2596.

by human macrophages. J Immunol *165*, 6406-6415.

lymphocyte glycoprotein T1/Leu-1. Nature *323*, 346-349.

scavenger receptor superfamily. Eur J Immunol *23*, 2320-2325.

receptors. Proc Natl Acad Sci USA *87*, 8810-8814.

proteins. J Biol Chem *272*, 6151-6158.

*378*, 129-144.

314.

232.

Evol Biol *10*, 181.

*409*, 198-201.

Krieger, M. (1990). An ancient, highly conserved family of cysteine-rich protein domains revealed by cloning type I and type II murine macrophage scavenger

P.C. (2008). Three-dimensional solution structure and conformational plasticity of the N-terminal scavenger receptor cysteine-rich domain of human CD5. J Mol Biol

cloning, mapping to human chromosome 1 q21-q23, and cell binding characteristics of Spalpha, a new member of the scavenger receptor cysteine-rich (SRCR) family of

J., and Lozano, F. (2003). The accessory molecules CD5 and CD6 associate on the

Oliveira, C., Sreenu, V.B., Evans, E.J., Davis, S.J., Moreira, A.*, et al.* (2009). Molecular cloning and analysis of SSc5D, a new member of the scavenger receptor cysteine-

haptoglobin-hemoglobin complexes from plasma. Int J Biochem Cell Biol *34*, 309-

scavenger receptor cysteine-rich type I transmembrane molecule (M160) expressed

family and its relationship to the bovine gamma delta T cell co-receptor WC1. BMC

cysteine-rich domain sheds light on an ancient superfamily. Nat Struct Biol *6*, 228-

Reid, K.B., Poustka, A., and Skjodt, K. (1999). Cloning of gp-340, a putative opsonin receptor for lung surfactant protein D. Proc Natl Acad Sci USA *96*, 10794-10799. Jones, N.H., Clabby, M.L., Dialynas, D.P., Huang, H.J., Herzenberg, L.A., and Strominger,

J.L. (1986). Isolation of complementary DNA clones encoding the human

Moestrup, S.K. (2001). Identification of the haemoglobin scavenger receptor. Nature

(1993). A new macrophage differentiation antigen which is a member of the

**11**

*U.S.A.*

**Assessing Multiple Sequence**

*1Department of Computer Science and Engineering*

*2School of Biological Sciences and 3Center for Plant Science Innovation University of Nebraska-Lincoln,*

**Alignments Using Visual Tools**

Catherine L. Anderson1, Cory L. Strope2 and Etsuko N. Moriyama2,3

Bioinformatics and molecular evolutionary analyses most often start with comparing DNA or amino acid sequences by aligning them. Pairwise alignment, for example, is used to measure the similarities between a query sequence and each of those in a database in BLAST similarity search, the most used bioinformatics tool (Altschul *et al.*, 1990; Camacho *et al.*, 2009). Evolutionary history among sequences can be reflected better when more than two sequences are aligned, in a multiple sequence alignment (MSA). When building an MSA, we assume that the sequences compared are derived from a common ancestral sequence. Then the process of MSA building is to infer homologous positions between the input sequences and place gaps in the sequences in order to align these homologous positions. These gaps represent evolutionary events of their own. Gaps (also called indels) are caused by either insertions or deletions of characters (nucleotides or amino acids) on a particular lineage of sequences during the evolution. Building an MSA is, therefore, to reconstruct the evolutionary history of the sequences involved. While it is easy to understand that the quality of MSAs affects the quality of phylogenetic tree reconstruction, the effect of MSA quality reaches far beyond this. Some examples of bioinformatics methods that utilize information extracted from MSAs include: profile building in similarity search (*e.g*., PSI-BLAST: Altschul *et al.*, 1997), motif/profile recognition (*e.g*., PROSITE: Hulo *et al.*, 2008), profile hidden Markov models for protein families/domains (*e.g.*, Pfam: Finn *et al.*, 2010), and protein secondary-structure prediction (for review, see Pirovano & Heringa, 2010). There are numerous bioinformatics and molecular evolutionary analyses that are affected by