Samir Touzani

Daniel Busby

Multivariate

Multivariate wavelet kernel regression method

Keywords: Wavelet, Reproducing kernel, Nonparametric regression, ANOVA, Landweber iterations

come

Introduction

In this work we consider the classical multivariate nonparametric regression problem with the aim of recovering an unknown function f from noisy data y i = f (x i ) + ǫ i , i = 1, . . . , n where x i = (x i 1 , . . . , x i d ) T is d dimensional vector of inputs, the ǫ i are i.i.d. and N (0, σ 2 ) random errors. The scope of this work is to take advantage of the multiresolution analysis provided by wavelet decomposition in multivariate function estimation. Indeed, in recent years there has been an important development in the application of wavelet methods in statistics, especially in signal processing, in image and function representation methods, with many successes in the efficient analysis and compression of noisy data. The multiresolution analysis provides a good time frequency localization, which makes wavelet methods particularly effective to estimate functions with sharp spikes, and discontinuities. Thus wavelets are used in various nonparametric regression methods. However most of these methods are implemented only for one (signal) or two (image) dimensional problems. The reason for this is that these algorithms assumed that data is dyadic and with equally spaced points. Several algorithms have been proposed to overcome the setting of non-dyadic and non-equispaced design. Among them, [START_REF] Antoniadis | Random design wavelet curve smoothing[END_REF] transforms the random design into equispaced data via a binning method. [START_REF] Kovac | Extending the scope of wavelet regression methods by coefficient-dependent thresholding[END_REF] apply the linear transformation to Email addresses: samirtouzani.phd@gmail.com (Samir Touzani), daniel.busby@ifpen.fr (Daniel Busby) the data to map it to a dyadic and equispaced set of points. [START_REF] Kerkyacharian | Regression in random design and warped wavelets[END_REF] project the data on an unusual non-orthonormal basis, called warped wavelet basis. [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF] suggested a regularization method relying on wavelet kernel reproducing Hilbert spaces, which does not require a pre-processing of data. The method also achieves optimal convergence rates in the Besov spaces when the estimation error is calculated at the design points only no matter how irregular the design is. Given that, it seems that this method is well adapted to be generalized for the multivariate regression using wavelets. Inspired by the component selection and smoothing operator (COSSO) [START_REF] Lin | Component selection and smoothing in smoothing spline analysis of variance models[END_REF], which is based on ANOVA (ANalysis Of VAriance) decomposition, and the wavelet kernel penalized estimation for non-equispaced design regression proposed by [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF], we introduce a new approach in the estimation of ANOVA components. Given a wavelet type expansion of f we consider a class of wavelet estimators for the nonparametric regression problem using a penalized least-squares approach. The penalties are chosen in order to control the smoothness of the resulting estimator. For this we use the same penalty as the one used for COSSO, in other words the semi-norm penalty. So we take for penalty, a weighted sum of wavelet details norms. In this paper, we first shortly review some definitions on wavelet. Then we present a new nonparametric regression method, named Wavelet Kernel ANOVA (WK-ANOVA). A new iterative projected shrinkage algorithm based on Landweber iterations is introduced in section 4. Finally, numerical tests are presented and discussed.

Multiresolution analysis and wavelet kernel

Consider the univariate regression problem:

y i = f (x i ) + ǫ i , i = 1, . . . , n (1)
where (x i ) i=1,...,n is the irregular design, the ǫ i are i.i.d. and N (0, σ 2 ) random errors and f an unknown regression function to be estimated.

Multiresolution analysis

We define a multi resolution analysis as a sequence of closed subspaces V j , j ∈ Z in L 2 (R) and which possesses the following properties:

1. ∩ j∈Z V j = 0, 2. ∪ j∈Z V j = L 2 (R), 3. ∀f ∈ L 2 (R), ∀j ∈ Z, f ∈ V j if and only if f (2x) ∈ V j+1 ; 4. ∀f ∈ L 2 (R), ∀k ∈ Z, f ∈ V 0 if and only if f (x -k) ∈ V 0 , 5. there exist a scaling function φ ∈ V 0 whose integer-translates x → φ(x -k) k∈(Z) span the space V 0 .
If we define P j f as the projection of a function f onto the space V j , this is expressed by

P j f = P j-1 f + w j-1
where the function w j-1 represents the residual between the two approximations on V j and on V j-1 . This function can be written in terms of dilated and translated wavelets:

w j-1 = k∈Z f, ψ j-1,k ψ j-1,k
where {ψ j,k (x) = 2 j/2 ψ(2 j x -k); k ∈ Z} is a set of functions that are orthogonal to each function of V j and span the space W j which is a detail space. Hence, an important property of multiresolution analysis can be defined as:

V j = V j-1 ⊕ W j-1
By periodizing an orthonormal basis for L 2 (R) we construct an orthonormal wavelet basis for L 2 ([0, 1]) generated by dilatation and translation of compactly supported scaling function φ per , and a compactly supported wavelet ψ per , where:

φ per j,k (x) = l∈Z φ j,k (x + l), ψ per j,k (x) = l∈Z ψ j,k (x + l) and V per j = span{φ per j,k , k ∈ Z}, W per j = span{ψ per j,k , k ∈ Z} The resulting orthogonal basis provides an orthogonal decomposition L 2 ([0, 1]) = V per 0 ⊕ W per 0 ⊕ W per 1 ⊕ ...
where V per 0 (spanned by φ 0,0 = ψ -1,0 ) consists of constant function and W per j is a 2 j dimensional space. To enhance the interpretability we omit in what follows the index per. For any integer j 0 ≥ 0 any function f ∈ L 2 ([0, 1]) is expressed in the form:

f (t) = 2 j 0 -1 k=0 α j 0 ,k φ j 0 ,k (t) + j≥j 0 2 j 0 -1 k=0 β j,k ψ j,k (t), t ∈ [0, 1]
where the scaling coefficients are α j 0 ,k = f, φ j 0 ,k (k = 0, 1, ..., 2 j 0 -1) and the wavelet coefficients are β j,k = f, φ j,k (j ≥ j 0 , k = 0, 1, ..., 2 j 0 -1). For more details on the mathematical aspects of wavelets and their applications in statistical settings we refer to [START_REF] Daubechies | Ten lectures on wavelets[END_REF], [START_REF] Vidakovic | Statistical modeling by wavelets[END_REF], [START_REF] Ogden | Essential wavelets for statistical applications and data analysis[END_REF] and Antoniadis et al. (2001).

Wavelet kernel

Let

G -1 = {-1} × {0}, G 0 = {0} × {0, 1} and for each integer J ≥ 1 let G J = {J} × {k ∈ {0, ..., 2 J }; k/2 /
∈ Z}, i.e. G J is the index set of wavelets at resolution level J. The whole set of indexes pairs (j, k) that describes all wavelets will be denoted by G = j≥-1 G j . Therefore, any function f ∈ L 2 ([0, 1)) admits the infinite wavelet expansion:

f = g∈G f g ψ g
where ψ g is the wavelet basis function indexed by g ∈ G, f g is the corresponding expansion coefficient and ψ -1,0 = φ 0,0 . We now define a class of wavelet-based Hilbert spaces. For any function:

Γ : G → [0, ∞)
define the Hilbert space:

H Γ = {f ∈ L 2 ([0, 1]) : g∈G Γ(g)|f g | 2 < ∞} with scalar product: f, h Γ = g∈G f g h g Γ(g)
and let be • Γ the associated norm. As G J is a finite subset of G, we have V J ⊂ H Γ for every J ≥ 1 defined as:

V J = V 0 ⊕ J-1 j=0 W j Moreover, for any f ∈ H Γ , lim J→∞ f -P J (f ) Γ = 0 ( 2
)
where P J (f ) is the projection of a function f into the space V J . The space H Γ is a RKHS and the corresponding reproducing kernels are given by:

K Γ (x, y) = g∈G ψ g (x) Γ(g) ψ g (y), x, y ∈ [0, 1)
By definition of the index set G, the kernel K can also be written as a sum of the reproducing kernels:

K Γ j (x, y) = 2 j -1 k=0 ψ j,k (x) Γ(j, k) ψ j,k (y)
This implies that the RKHS H Γ , can be decomposed into a direct sum of wavelet RKHS's (spanned by a set of wavelets of scale j) as:

H Γ = V 0 ⊕ j≥0 W Γ j (3)
where each space W Γ j is the RKHS associated to the kernel K Γ j . This representation involves an infinite decomposition of the detail space, in practice we truncate (3) up to a maximum resolution J, in other words, the RKHS H J,Γ = V 0 ⊕ J j=0 W Γ j defines a multiresolution analysis of H Γ and the associated kernel is:

K Γ J (x, y) = g∈∪ 0≤j≤J G j ψ g (x) Γ(g) ψ g (y), x, y ∈ [0, 1) Furthermore, from (2) lim J→∞ K Γ -K Γ J ∞ = 0
We assume that Γ is only a function of j and equals to 2 2js on G j and s > 1/2, then H Γ equals to the Sobolev space B s 2,2 ([0, 1]) of index s. For more mathematical details we refer to [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF].

The wavelet kernel ANOVA

Consider now a multivariate regression problem:

y i = f (x i ) + ǫ i , i = 1, • • • , n
where x i = (x i 1 , . . . , x i d ) T is d dimensional vector of inputs, the ǫ i are i.i.d. and N (0, σ 2 ) random errors and f an unknown multivariate regression function to be estimated.

Definition

The idea behind the well known smoothing spline ANOVA model [START_REF] Wahba | Spline models for observational data[END_REF] is to construct a RKHS F = {f ∈ L 2 ([0, 1] d )} corresponding to the ANOVA decomposition:

f (X (1) , • • • , X (d) ) = f 0 + d l=1 f l (X (l) ) + l<m f lm (X (l) , X (m) ) + . . . + f 1,2,...,d (X (1) , • • • , X (d) ) (4)
where f 0 is a constant, f j 's are univariate functions representing the main effects, f jl 's are bivariate functions representing the two way interactions, and so on. Then the model space F is the tensor product space of H l Γ :

F = d l=1 H l Γ = {1} ⊕ H l Γ ⊗ H m Γ = (V l 0 ⊗ V m 0 ) J-1 j=0 Γ -2 j (W l j ⊗ W m j ) (7)
It's easy to see that V 0 is also the subspace of L 2 ([0, 1]) spanned by the constant function on [0, 1], one has V 0 = V l 0 ⊗ V m 0 = {1}. Thus, the function space F, which is a wavelet-based RKHS, can also be written as:

F = {1} ⊕ q γ=1 F γ (8)
where the F γ is an orthogonal subspaces of F and correspond to the subspaces Hl Γ , [ Hl Γ ⊗ Hm Γ ], etc. . . In the additive model q = d where d is the number of input parameters and in the model with two way interactions q = d(d + 1)/2. We assume that a second order ANOVA expansion gives a satisfactory approximation of f . We denote by P Γ γ f the orthogonal projection of f onto Γ -1 j W j and • the norm in the RKHS Γ -1 j W j . Under the framework of smoothing spline ANOVA one way to estimate f is to find f ∈ F that minimizes:

1 n n i=1 {y i -f (x i )} 2 + λ 2 q γ=1 J-1 j=0 2 j -1 k=0 θ -1 γ,j,k P Γ γ,j,k f 2 (9)
where θ γ,j,k ≥ 0. If θ γ,j,k = 0, then the minimizer is taken to satisfy P Γ γ,j,k f 2 = 0, using the convention 0/0 = 0. The parameter λ controls the trade-off between the first term in the above expression which discourages the lack of fit of f and the second one which penalizes the roughness of f . In analogy with COSSO [START_REF] Lin | Component selection and smoothing in smoothing spline analysis of variance models[END_REF] and wavelet kernel penalized estimation [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF] we propose the WK-ANOVA procedure, another way to estimate f , given by f ∈ F that minimize:

1

n n i=1 {y i -f (x i )} 2 + λ 2 R q (f ) (10) with R q (f ) = q γ=1 J-1 j=0 2 j -1 k=0 P Γ γ,j,k f is a sum of wavelet-based RKHS norms, instead of the squared RKHS norm employed in (9). We note that R q (f ) is not a norm in F but a pseudo-norm in the following sense: R q (f ) ≥ 0, R q (cf ) = |c|R q (f ), R q (f + h) ≤ R q (f ) + R q (h) ∀ f, h ∈ F, and, R q (f ) > 0 for any non constant f ∈ F. Furthermore B ≤ R q (f ) 2 ≤ qB (11) with B = q γ=1 J-1 j=0 2 j -1 k=0 P Γ γ,j,k f 2 (12)
Note that there is only one smoothing parameter λ which should be properly chosen, instead of multiple smoothing parameters θ's in (9).

The existence of the WK-ANOVA estimate, which is due to the convexity of ( 10), is guaranteed by adapting Theorem 1 of [START_REF] Lin | Component selection and smoothing in smoothing spline analysis of variance models[END_REF].

Theorem 3.1. Let F be the wavelet-based RKHS of functions over [0, 1] d . Assume that F can be decomposed as ( 8). There exists a minimizer of ( 10). Define • n as the Euclidian norm in R n . Under our previous assumption, s > 1/2 and Γ j = 2 2sj . The following theorem is equivalent to theorem 2 of [START_REF] Lin | Component selection and smoothing in smoothing spline analysis of variance models[END_REF] and shows that the WK-ANOVA estimator in the additive model has a rate of convergence n -s/(2s+1) , where s is the order of smoothness of the components.

Theorem 3.2. Consider the regression model y

i = f 0 (x i ) + ε i , i = 1, ..., n, where x i 's are given deterministic points in [0, 1] d , and the ε i 's are i.i.d. N (0, σ 2 ) noise variables. Assume f 0 lies in F = {1} ⊕ d l=1 H l Γ , with H γ Γ = {1} ⊕ Hγ Γ being the Sobolev space B s 2,2 ([0, 1]) of index s.
Consider the WK-ANOVA estimator f at the design points as defined by (10).

Then (i) if f 0 is not a constant, and λ -1 n = O p (n s/(2s-1) )R (2s-1)/(4s+2) q (f 0 ), we have f - f 0 n = O p (λ n )R 1/2 q (f 0 ); (ii) if f 0 is a constant, we have f -f 0 n = O p (max{n -s/(2s-1) λ -2/(2s-1) n , n -1/2 }).
The following Lemma shows that the solution of ( 10) is in finite dimensional space and the WK-ANOVA estimate can be computed directly from (10) by linear programming techniques.

Lemma 3.3. Let f = b + q γ=1 fγ be a minimizer of ( 10), with f γ ∈ F γ . Then fγ ∈ span{K γ (x i , •), i = 1, ..., n}, where K γ is the reproducing kernel of the space F γ Using the suggestion of [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF] for solving penalized problems with l 1 penalty, we can give an equivalent formulation of (10) for computational consideration. Consider the problem of finding θ = {θ γ,j,k , γ = 1, ..., q; j = 0, ..., J -1; k = 1, ..., 2 j -1} and f ∈ F to minimize:

1 n n i=1 {y i -f (x i )} 2 + λ 0 q γ=1 J-1 j=0 2 j -1 k=0 θ -1 γ,j,k P Γ γ,j,k f 2 + ν q γ=1 J-1 j=0 2 j -1 k=0 θ γ,j,k (13)
subject to θ γ,j,k ≥ 0, and where λ 0 is a fixed positive constant and ν is a smoothing parameter.

We fix λ 0 at some value. Then 13). (ii) On the other hand, if a pair ( θ, f ) minimizes ( 13), then f minimizes (10).

Lemma 3.4. Set ν = λ 4 /(4λ 0 ). (i) if f minimizes (10), set θγ,j,k = λ 1/2 0 ν -1/2 P Γ γ,j,k f , then the pair ( θ, f ) minimizes (
As already introduced by Amato et al. ( 2006) we can penalize the norm of coefficients by blocks, which allows reducing the number of θ's that need to be estimated and can provide a better regularization. Hence, as defined before:

K Γ jm (x, y) = k∈T jm ψ j,k (x) Γ j ψ j,k (y)
where m = 1, ..., M j with M j denotes the number of blocks at resolution j, and T jm the blocks of length L jm at resolution j. In the same way consider the decomposition:

H l Γ = V 0 ⊕ J-1 j=0 m Γ -1 j W l j,m
then replace (13) by:

1 n n i=1 {y i -f (x i )} 2 + λ 0 q γ=1 J-1 j=0 M j m=1 θ -1 γ,j,m P Γ γ,j,m f 2 + ν q γ=1 J-1 j=0 M j m=1 θ γ,j,m (14)
We can note that the form of ( 14) is similar to the smoothing spline ANOVA (9) with multiple smoothing parameters and an additional penalty on the θ's. There is only one smoothing parameter ν in ( 14) and θ's are part of the estimate, rather than three smoothing parameters.

For the WK-ANOVA procedure the sparsity on the detail components is controlled by the additional penalty on θ's in ( 14) makes possible to have some θ's to be zero, thus producing a sparse kernel estimate in sense of [START_REF] Gunn | Structural modeling with sparse kernels[END_REF].

Algorithm

We will use an iterative optimization algorithm which is equivalent to the one used in [START_REF] Lin | Component selection and smoothing in smoothing spline analysis of variance models[END_REF] and [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF]. On each step of iteration, for any fixed θ we minimize ( 14) with respect of f , and then for this choice of f we minimize ( 14) with respect of θ. Note that for any fixed θ (14) is equivalent to the smoothing spline ANOVA procedure. Therefore from [START_REF] Wahba | Spline models for observational data[END_REF] the solution f of ( 14) has the following form

f (x) = b + n i=1 c i q γ=0 J-1 j=0 M j m=1 θ γ,j,m K Γ γ,j,m (x i , x) (15)
Where

c = (c 1 , ..., c n ) T , b ∈ R, K Γ γ,j,m is the reproducing kernel of Γ -1 j W l γ,j,m if γ ≤ d and is the reproducing kernel of Γ -2 j W l γ,j,m ⊗ W p γ,j,m else.
In what follows, we denote by

K Γ γ,j,m the n × n matrix {K Γ γ,j,m (x i , x t )}, i = 1, ..., n, t = 1, ..., n, by K Γ θ the matrix q γ=0 J-1 j=0 M j m=1 θ γ,j,m K Γ γ,j,m (x i , x
) and 1 n be the column vector consisting of n ones. Then we can write f = K Γ θ c + b1 n , it follows that ( 14) can be expressed as

1 n Y - q γ=0 J-1 j=0 M j m=1 θ γ,j,m K Γ γ,j,m c -b1 n 2 n + λ 0 c T K Γ θ c + ν q γ=0 J-1 j=0 M j m=1 θ γ,j,m (16)
where θ γ,j,m ≥ 0, γ = 1, ..., q, j = 0, ..., J -1, m = 1, ..., M j . If θ's are fixed, then ( 16) can be written as

min c,b Y -K Γ θ c -b1 n 2 n + nλ 0 c T K Γ θ c (17)
which is similar to the smooting spline problem (a quadratic minimization problem) and the solution satisfy (for more details see [START_REF] Wahba | Spline models for observational data[END_REF]):

(K Γ θ + nλ 0 I)c + b1 n = Y 1 T n c = 0
where I is the identity matrix. Let's fix b and c at their values from (17), denote d γ,j,m = K Γ γ,j,m c, and let D be the n × ( γ j (2 j -1)) matrix with the (γ, j, m)th column being d γ,j,m . The θ that minimizes ( 16) is the same as the solution to

min θ z -Dθ 2 n + nν q γ=0 J-1 j=0 M j m=1 θ γ,j,m subject to θ γ,j,m ≥ 0 ( 18
)
where

z = Y -(1/2)nλ 0 c -b1 n .
By starting from a simpler estimate such as the one obtained by penalized least squares with quadratic penalties on the coefficients, a one step update procedure is sufficient to improve the WK-ANOVA estimator. Then we propose a one step update procedure:

1. Initialization: Fix θ γ,j,m = 1, γ = 1, ..., q, j = 0, ..., J -1, m = 1, ..., M j . 2. Tune λ 0 using v-fold-cross-validation.

3. Solve for c and b with (17).

4. For each fixed ν, solve (18) with the c and b obtained in step 3. Tune ν using v-foldcross-validation. The θ's corresponding to the best ν are the final solution at this step. 5. With the new θ tune λ 0 using v-fold-cross-validation. 6. With the new θ and λ 0 , solve for c and b with (17)

A discussion of a one step procedure and fully iterated procedure can be found in [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF]. The performance of the WK-ANOVA estimator depends on the smoothing parameter ν and the chosen resolution J. The choice of these parameters obviously involves an arbitrary decision. In our work we will fix J = log 2 n, but by varying the resolution level we can explore features of the data arising on different scales. We will use v fold cross validation to tune ν. It seems reasonable to take v equal to 5. We also choose to use compactly supported wavelets, it follows that the numerical algorithm for the kernel computation is based on Daubechies cascade procedures [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. Specifically, the cascade algorithm computes the values of wavelets at dyadic points. In order to evaluate the kernel matrices K Γ γ,j,m the values of the wavelets have been computed on a fine dyadic grid and stored in a table. Values of wavelets at arbitrary points, necessary for evaluation of K Γ γ,j,m , were then computed by considering the value at the closest point on the tabulated grid. The table construction of wavelet kernel matrices requires O(n 2 S) elementary operations where S denotes the length of the wavelet filter. However, the table is constructed once and stored in memory. In addition, as the dimension of the problem grows, the number of matrices K Γ γ,j,m also grows as well, and because of the v-fold-cross-validation these matrices must be re-computed several times. All this increases significantly the computational time, and therefore it is necessary to compute the matrices once and stored them in memory. The formulation in ( 18) is a high-dimensional nonnegative garrote (NNG) optimization problem introduced by [START_REF] Breiman | Better subset regression using the nonnegative garrote[END_REF] for which there exists a variety of algorithms to find the solution. However, we introduce in the next section a new iterative shrinkage algorithm. Even if the number of iterations may be high to compete with modern high-dimensional optimization algorithms, the proposed iterative algorithm is conceptually simple and easy to implement.

Iterative projected shrinkage algorithm

Recently, iterative shrinkage/theresholding (IST) algorithm tailored to solve the LASSO regression problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], has been proposed independently by several research groups (among them [START_REF] Figueiredo | An em algorithm for wavelet-based image restoration[END_REF] and [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]). This algorithm combines the Landweber iteration [START_REF] Landweber | An iterative formula for fredholm integral equations of the first kind[END_REF] and the soft thresholding method. The main advantages of this algorithm are its conceptual simplicity, its easiness of implementation and only involves matrix-vector multiplication. In this section, we propose a modified version of IST algorithm that solves the NNG regression problem.

Definition

Consider the (18) regression problem:

min θ z -Dθ 2 n + nν q γ=0 J-1 j=0 M j m=1
θ γ,j,m subject to θ γ,j,m ≥ 0

The functional ( 18) is convex since matrix D T D is symmetric and positive semidefinite and since the constraints θ γ,j,m > 0 define also a convex feasible set. For the convex optimization problem, the Karush-Kuhn-Tucker (KKT) conditions are necessary and sufficient for the optimal solution θ * , where θ * = argmin θ z -Dθ 2 + nν q γ=0 J-1 j=0 M j m=1 θ γ,j,m subject to θ γ,j,m ≥ 0. This KKT conditions are defined as:

{-d T γ,j,m (z -Dθ * ) + ν}θ * γ,j,m = 0 -d T γ,j,m (z -Dθ * ) + ν ≥ 0 θ * γ,j,m ≥ 0 which is equivalent to -d T γ,j,m (z -Dθ * ) + ν = 0, if θ * γ,j,m = 0 (19) -d T γ,j,m (z -Dθ * ) + ν > 0, if θ * γ,j,m = 0 ( 20
)
where d γ,j,m denotes the (γ, j, m)th column of D. Therefore, from ( 19) and ( 20) we can derive the fixed-point equation:

θ * = P Ω + (δ Sof t ν (θ * + D T (z -Dθ * ))) (21)
where P Ω + is the nearest point projection operator onto the nonnegative orthant (closed convex set) Ω + = {x : x ≥ 0} and δ Sof t λ is the soft-thresholding function defined as

δ Sof t ν (x) =    0 if |x| ≤ ν x -ν if x > ν x + ν if x < -ν
Thus, we propose an iterative algorithm that we name the iterative projected shrinkage algorithm (IPS) and which is defined by

θ [p+1] = P Ω + (δ Sof t ν (θ [p] + D T (z -Dθ [p] ))) (22)
The following theorem concerns the convergences of IPS algorithm:

Theorem 4.1. IPS algorithm defined by ( 22) converge to the solution of ( 18), whenever such solution exists, for any starting vector θ [0] .

The proof of this theorem can be found in the Appendix. We have assumed that λ max (D T D) ≤ 1 (where λ max is the maximum eigenvalue). Otherwise we solve the equivalent minimization problem

min θ z α - D α θ 2 + nν α q γ=0 J-1 j=0 M j m=1 θ γ,j,m subject to θ γ,j,m ≥ 0
where the positive constant α ensures that λ max (D T D) ≤ 1.

Stopping conditions IPS algorithm is an iterative procedure which produces a sequence of solutions θ [0] , θ [0] , . . . , θ [p] converging to the optimal solution θ * . There is a need to stop the algorithm when the solution θ [p] is sufficiently close to the optimal solution θ * . Several stopping conditions have been proposed in the literature (for example [START_REF] Defrise | A note on stopping rules for iterative regularization methods and filtered svd[END_REF]). We choose to use a stopping condition based on the KKT conditions, which are easy to evaluate. This ǫ-KKT conditions are defined as

d T γ,j,m (Y -Dθ * ) = ν -ǫ, if θ * γ,j,m = 0 d T γ,j,m (Y -Dθ * ) ≤ ν -ǫ, if θ * γ,j,m = 0
where ǫ > 0 is a constant which defines the precision of the solution.

Simulations

In this section we will study the empirical performance of WK-ANOVA, in terms of prediction accuracy. The measure of the prediction accuracy is given by Q 2 which is defined as

Q 2 = 1 - ntest i=1 (y i -f (x i )) 2 ntest i=1 (y i -ȳ) 2 , with n test = 500 (23)
where y i denotes the ith test observation of the test set, ȳ is their empirical mean and f (x i ) is the predicted value. We compare the obtained results with those obtained by COSSO and Gaussian Process (GP) [START_REF] Busby | Hierarchical adaptive experimental design for gaussian process emulators[END_REF]. We also compare the methods for different experimental design sizes, uniformly distributed on [0, 1] d and built by Latin Hypercube Design procedure [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] with maximin criterion (Santner et al., 2003) (maximinLHD). Moreover, different signal to noise ratio were applied SN R ≡ 1 : 3 (high noise) SN R ≡ 1 : 7 (medium noise) and SN R ≡ ∞ (without noise), with SN R = [V ar(f (X))]/σ 2 . For each setting of test examples 2 and 3, we perform 50 times the test. We fixed M j = 2 j -1 for γ = 1, . . . , d and M j = 1 for γ > d, in other words we penalize by the translation parameter k for the main effects and by resolution j for the interaction. This assumption permits us reducing significantly the computational time. The wavelets used in our tests were Daubechies wavelets with 3 vanishing moments [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. The WK-ANOVA was implemented using R. We run the simulations on a computer operated by 32 bits-Windows OS, this latter imposes limits on the total memory allocation. Knowing that the storage of the matrices K Γ γ,j,m is memory consuming, we limit the dimension to our examples to 8 and the sample size to estimate Q 2 to 500. To fit COSSO models we have used a modified version of the original Matlab code provided by Yi Lin and Hao Helen Zhang. We made our modification of this algorithm to make it faster. The GP code was implemented using R with a generalized power exponential Family [START_REF] Busby | Hierarchical adaptive experimental design for gaussian process emulators[END_REF].

Example 1

Let's consider an additive model with X ∈ [0, 1] 6 , with the following function

f (X) = g 1 (X (1) ) + g 2 (X (2) ) + g 3 (X (3) ) + g 4 (X (4) )
where g 1 (t) = t; g 2 (t) = (2t -1) 2 ; g 3 (t) = sin(2πt) 2 -sin(2πt) ;

g 4 (t) = 0.1 sin(2πt) 0.2 cos(2πt) + 0.3 sin 2 (2πt) + 0.4 cos 3 (2πt) + 0.5 sin 3 (2πt) Therefore X (5) , X (6) are uninformative. We use an experimental design of size n = 200, built by maxminLHD, and SN R ≡ ∞. Figure 1 gives the plot of data observation with the true ANOVA component f l and theirs WK-ANOVA estimates against inputs X (l) , l = 1, . . . , 6.

The Q 2 of this WK-ANOVA estimate is equal to 0.96 which is a good performance. However, we can note that the estimation of the linear function component f 1 does suffer from using a wavelet method. Part of the reason is the boundary effects caused by using periodic wavelets.

Example 2

In this first test case, consider an additive model with X = [0, 1] 8 , with the following function

f (X) = g 1 (X (1) ) + g 2 (X (2) ) + g 3 (X (3) ) + g 4 (X (4) ) + ǫ
where g 1 (t) = 0.1 sin(2πt) + 0.2 cos(2πt) + 0.3 sin(2πt) + 0.4 cos 3 (2πt) + 0.5 sin 3 (2πt)

g 2 (t) = (2t -1) 2 g 3 (t) = | sin(3πt)| + 0.5| sin(5πt)| 2 -sin(4πt) g 4 (t) = | sin(2πt)| 2 -sin(2πt)
Therefore X (5) , . . . , X (8) are uninformative. Note that all informative input of f have nonlinear response, in addition g 3 and g 4 have discontinuities on the derivative. The true ANOVA components f l , l = 1, . . . , n with their WK-ANOVA, COSSO and GP estimates are given in figure 2. These estimates were built with an experimental design of size n = 200 and with noise ratio SN R = 3. The WK-ANOVA has more fidelity to the reality than COSSO and GP especially for the components f 3 and f 4 . Indeed, WK-ANOVA captures more the discontinuities of the components f 3 and f 4 . This good fit is due to the properties of wavelets analysis. In other words, our algorithm based on wavelets is well suited to this type of functions (with discontinuities on the derevatives). We run the simulation 50 times for different sizes of experimental design (n = 50, 100, 200) and different signal to noise ratio SN R ≡ 1 : 3, SN R ≡ 1 : 7 and SN R ≡ ∞. The results are summarized in figure 3 each panel is a boxplot of the 50 estimations of Q 2 . As expected, the accuracy of WK-ANOVA estimates increases when the sample size raise. We can see that WK-ANOVA procedure outperforms COSSO and GP in all the studied settings. Moreover, even though there are much more parameters to estimate with WK-ANOVA comparing to COSSO, this procedure does not seem to suffer from small sample size effect. For this example, WK-ANOVA has shown better denoising properties and better predictivity. In addition, for n = 100 and n = 200 WK-ANOVA is the most robust.

Example 3

Consider the g-Sobol function, which is strongly nonlinear and is described by a nonmonotonic relationship. This function is a well-known test case in the studies of global sensitivity analysis. Figure 4 illustrates the g-Sobol function against the two most influential x ( (

0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 X (1) f1(X (1) ) f1(X (1) ) WK-ANOVA 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 X (2) f2(X (2) ) f2(X (2) ) WK-ANOVA 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 X (3) f3(X (3) ) f3(X (3) ) WK-ANOVA 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 X (4) f4(X (4) ) f4(X (4) ) WK-ANOVA 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 X (5) f5(X (5) ) f5(X (5) ) WK-ANOVA 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 X (6) f6(X (6) ) f6(X (6) ) WK-ANOVA
) )

f1( (x ( (1) ) ) ) WK-ANOVA COSSO GP 0.0 0.2 0.4 0.6 0.8 1.0 -0.5 0.0 0.5 1.0 x ( (2) ) f2( (x ( (2) )
) )

f2( (x ( (2) ) ) )

WK-ANOVA COSSO GP 0.0 0.2 0.4 0.6 0.8 1.0 -0.5 0.0 0.5

1.0 x ( (3) ) f3( (x ( (3) )
) )

f3( (x ( (3) ) ) )

WK parameters X (1) and X (2) . The g-Sobol function [START_REF] Saltelli | Sensitivity analysis[END_REF]) is defined for 8 inputs factors as

g Sobol (X (1) , . . . , X (8) ) = 8 k=1 g k (X (k) ) with g k (X (k) ) = 4X (k) -2 + a k 1 + a k
where {a 1 , . . . , a 8 } = {0, 1, 4.5, 9, 99, 99, 99, 99}. The contribution of each input X (k) to the variability of the model is represented by the weighting coefficient a k . The lower this coefficient a k , the more significant the variable X (k) . For example: k) is non significant.

           a k = 0 → x (k) is very important, a k = 1 → x (k) is relatively important, a k = 4.5 → x (k) is poorly important, a k = 9 → x (k) is non important, a k = 99 → x (
We run the simulation 50 times for different sizes of experimental design (n = 50, 100, 200) and different signal to noise ratio SN R ≡ 1 : 3, SN R ≡ 1 : 7 and SN R ≡ ∞. The results are summarized in figure 5 each panel is a boxplot of the 50 estimations of Q 2 . We can see that for all the tested experimental design sizes and noise ratios WK-ANOVA outperforms COSSO and GP. Moreover, the accuracy of the prediction is very good ( Q2 = 0.98, where Q2 is the average of the Q 2 's) even with n = 50 for the setting without noise, when a design with n = 200 design is necessary to perform a response surface with Q2 = 0.94 for GP and with Q2 = 0.90 for COSSO. Clearly, for this example WK-ANOVA has the best results in term of predictivity, denoising property and robustness.

Conclusion

In this article, we introduced a new regularized nonparametric regression method that we named WK-ANOVA. Differently than other wavelet methods, WK-ANOVA does require neither an equispaced experimental design points, nor a dyadic size of data. For the tested analytical examples which contains some discontinuities, WK-ANOVA outperforms COSSO and GP. However, in example 1 the wavelet methods have undesired boundary in the case of the estimation of nonperiodic function, which results from the use of periodic wavelets. One future investigation topic is to use boundary adapted wavelets.

Proofs

Proof of Theorem 3.1. Consider the following decomposition of the wavelet-based RKHS F F = {1} ⊕ H q where H q = q γ=1 F γ . Denote by A(f ) the functional to be minimized in (10). A(f ) is convex and continuous. By inequality (11) we have R q (f ) ≥ f F for any f ∈ H q . Let K Hq be the reproducing kernel on H q and •, • Hq be the inner product of H q . Denote by e n = max n i=1 (K Hq ) (1/2) (x i , x i ). By the definition of the reproducing kernel and the properties of Γ, we have for any f ∈ H q and i = 1, ..., n

|f (x i )| = | f, K Hq(x i ,•) Hq | ≤ e n f F ≤ e n R q (f )
Let consider D a closed, convex and bounded set defined as:

D = {f ∈ F; f = b + f 1 , with b ∈ {1}, f 1 ∈ H q , R q (f ) ≤ v, |b| ≤ v 1/2 + (e n + 1)v} where v = max i {y 2 i + |y i | + 1}.
Therefore by the theorem 4 of [START_REF] Tapia | Nomparametric Probability Density Estimation[END_REF]. There exist a minimizer f of (10) in D and A(f ) ≤ A(0) ≤ v. On the other hand, for any f

∈ F with R q (f ) > v clearly A(f ) ≥ R q (f ) > v; for any f ∈ F, f = b + f 1 with b ∈ {1}, f 1 ∈ H q , R q (f ) ≤ v and |b| > v 1/2 + (e n + 1)v, we therefore have |b + f 1 (x i ) -y i | > (v 1/2 + (e n + 1)v) -e n v -v = v 1/2
Hence A(f ) > v, and for any f / ∈ D, we have A(f ) > A(f ), wich proves that f is the minimizer of (10) in F.

Proof of Theorem 3.2. The condition on the unknown regression function f 0 are only active for its wavelets coefficients and do not include the V 0 scaling coefficients of f 0 . For any f ∈ F, write f

(x) = b + f 1 (x (1) ) + ... + f d (x (d) ) = b + g(x), such that n i=1 f l (x (l) i ) = 0, l = 1, ..., d and where b ∈ {1} and g ∈ d l=1 H l Γ . Similarly, write f 0 (x) = b 0 + g 0 (x), such that g 0 ∈ d l=1 H l Γ . By construction n i=1 {g 0 (x i ) -g(x i )} = 0, we can write A(f ) as : (b -b 0 ) 2 + 2 n (b -b 0 ) n i=1 ε i + 1 n ) n i=1 (g 0 (x i ) + ε i -g(x i )) 2 + λ 2 n R q (g)
Therefore, the minimizing b is b = b 0 + 1/n n i=1 ε i , which shows that b converges towards b 0 at rate n -1/2 . On the other hand, ĝ must minimize over d l=1 H l Γ the functional

1 n n i=1 {g 0 (x i ) + ε i -g(x i )} 2 + λ n R q (g) Let G = {g ∈ F : g(x) = f 1 (x (1) ) + ... + f d (x (d) ), with n i=1 f l (x (l)
i ) = 0, l = 1, ..., d}. then g 0 ∈ G and ĝ ∈ G. The conclusion of theorem 2 follows from the following Lemma.

Lemma Appendix .1. (Theorem 10.2 of [START_REF] Van De Geer | Empirical Processes in M-Estimation[END_REF], lemma 5.1 of [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF] and lemma 3 of [START_REF] Lin | Component selection and smoothing in smoothing spline analysis of variance models[END_REF] ) Let H ∞ (δ, G) be the δ-entropy of G for the supremum norm. Then H ∞ (δ, {g ∈ G : R q (g) ≤ 1} ≤ Ad (s+1)/s δ -1/s , for all δ > 0, n ≥ 1, and some A > 0 and 0 < 1/s < 2.

Proof of Lemma Appendix .1. Define G l as the set of univariate function of x (l) .

G l = {f l ∈ H l Γ : R q (f l ) ≤ 1, n i=1 f γ (x i ) (l) = 0}
It follows from Lemma 5.1 of [START_REF] Amato | Wavelet kernel penalized estimation for nonequispaced design regression[END_REF] that

H ∞ (δ, G l ) ≤ Aδ -1/s
for all δ > 0, and n ≥ 1, some A > 0 and 0 < 1/s < 2. By definition of G we see that in terms on the supreme norm, if each G l , l = 1, .., d can be covered by N balls of radius δ, then the set {g ∈ G : R q (g) ≤ 1} can be covered byN d balls with radius dδ, and we get :

H ∞ (dδ, {g ∈ G : R q (g) ≤ 1}) ≤ Adδ -1/s

Proof of Lemma 3.3. For any f ∈ F, write f = b + q γ=1 J-1 j=0

2 j -1 k=0 f γ,j,k with b ∈ {1} and f γ,j,k ∈ W Γ γ,j,k . Let the projection of f γ,j,k onto span{K Γ γ,j,k (x i , •), i = 1, ..., n} ⊂ W Γ γ,j,k be denoted by α γ,j,k and the orthonormal complement by β γ,j,k . Then f γ,j,k = α γ,j,k + β γ,j,k and (10) can be written as

1 n n i=1 {y i -b - q γ=1 J-1 j=0 2 j -1 k=0 K Γ γ,j,k (x i , •), α γ,j,k } 2 + λ 2 q γ=1 J-1 j=0 2 j -1 k=0 ( α γ,j,k 2 + β γ,j,k 2 ) 1/2
Therefore any minimizing f must be such that β γ,j,k = 0, and the result follows immediately.

Proof of Lemma 3.4. Denote the functional in ( 14) by B(θ, f ). For any γ = 1, ..., q; j = 0, ..., J -1; k = 1, ..., 2 j -1, we have

λ 0 θ -1 γ,j,k P Γ γ,j,k f 2 W Γ γ,j,k + νθ γ,j,k ≥ 2λ 1/2 0 ν 1/2 P Γ γ,j,k f Γ γ,j,k = λ 2 P Γ γ,j,k f W Γ γ,j,k
.

for any θ γ,j,k ≥ 0 and f ∈ F, and the equality holds if and only if θ γ,j,k = λ 1/2 0 ν -1/2 P Γ γ,j,k f W Γ γ,j,k . Therefore B(θ, f ) ≥ A(f ), where A(f ) denote the functional of (10) for any θ γ,j,k ≥ 0, γ = 1, ..., q; j = 0, ..., J -1; k = 1, ..., 2 j -1 and f ∈ F, with the equality holds only if θ γ,j,k = λ 1/2 0 ν -1/2 P Γ γ,j,k f W Γ γ,j,k

. The conclusion then follows.

Proof of Theorem 4.1. The orthogonal projection P Ω x of x onto Ω + is characterized by the following useful inequality: for all a ∈ Ω + and all x we have a -P Ω + x, P Ω + x -x ≥ 0 (.1)

From the inequality (.1) we can say that for any Ω + , x and z, we have P Ω + z -P Ω + x, P Ω + x -x ≥ 0 and P Ω + z -P Ω + x, z -P Ω + z ≥ 0 (.2)

Adding, we obtain P Ω + z -P Ω + x, z -x ≥ P Ω + z -P Ω + x 2 (.3)

From the Cauchy inequality we conclude that The convergence of the IPS algorithm follows from the following Lemma.

Lemma Appendix .2. (Lemma 3.4 of [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]) S ν is nonexpansive, i.e., for all x and z ∈ R,

δ Sof t ν (x) -δ Sof t ν (z) ≤ x -z (.6)
Since (.9) is convex the Karush-Kuhn-Tucker theorem suggests that a necessary and sufficient condition for θ * to be the solution of model ( 18) is that there is λ ≥ 0 such that, for any γ = 1, ..., q, j = 0, ... γ,j,m = 0 and that the condition (.9) is always satisfied. Now we consider the second possibility. If θ * γ,j,m is a fixed of the map T , with T x = P Ω + (S ν (x)), that is, T θ * γ,j,m = θ * γ,j,m . By (.10), we have

d T γ,j,m (Y -Dθ * ) -λ = 0
The conclusion follows immediately.

Figure 1 :

 1 Figure1: Plot of the six true functional components, f l , l = 1, . . . , 4 along with the data observations and theirs estimates given by WK-ANOVA for a realization from example 1

Figure 2 :Figure 3 :

 23 Figure 2: Plot of f l , l = 1, . . . 4 along with theirs estimates given by WK-ANOVA, COSSO and GP for a realization from example 2

Figure 4 :

 4 Figure4: Plot of g-Sobol function versus inputs X (1) and X (2) with other inputs fixed at 0.5

Figure 5 :

 5 Figure 5: Q2 results from example 3

P

  Ω + x -P Ω + z ≤ x -z (.4) Let E be nonempty set of all θ ∈ Ω + at which the functional (18) attains its minimum value over Ω + and θ * a member of E. Then θ * = P Ω + (δ Sof t ν (θ * )) andθ *θ [p+1] = P Ω + (δ Sof t ν (θ * )) -P Ω + (δ Sof t ν (θ [p] )) ≤ δ Sof t ν (θ * ) -δ Sof t ν (θ [p] ) (.5)

  , J -1, m = 1, ..., M j {-d T γ,j,m (Y -Dθ * ) + ν}θ * γ,m + d T γ,j,m (Y -Dθ [p] ) -ν (.10)It is not difficult to show that the KKT conditions are satisfied when θ [p+1]

d l=1 Hl Γ ⊕ l<m [ Hl Γ ⊗ Hm Γ ]...(5)where H l Γ = {1} ⊕ Hl Γ and Hl Γ is the RKHS associated to the first-order component functions f l of ANOVA expansion. The tensor products [ Hl Γ ⊗ Hm Γ ] is associated to the second-order component function f lm . We denote by W l j (a detail space at scale j) the RKHS associated to wavelet kernel K j = 2 j -1 k=0 ψ j,k (x)ψ j,k (y) and the variate X (l) , thereby for a fixed maximum resolution J the function space H l Γ can be written as:H l Γ = V 0 ⊕ J-1 j=0 Γ -1 j W l j(6)and the tensor product [H l Γ ⊗ H m Γ ] as:

Acknowledgements

The authors are grateful to Professor Anestis Antoniadis for many useful suggestions and helpful discussions.


Dr Sven Cichon
email: sven.cichon@uni-bonn.de

Lutz Priebe

Franziska A Degenhardt

Stefan Herms

Britta Haenisch

Manuel Mattheisen

Vanessa Nieratschker

Moritz Weingarten

Stephanie H Witt

René Breuer

T Paul

M Alblas

S Moebus

M Lathrop

M Leboyer

S Schreiber

M Grigoroiu-Serbanescu

W Maier

P Propping

M Rietschel

M M Nöthen

T W Mühleisen

BMP3, C4ORF22, CCDC102B, CHRNA9, CNTN4, ERBB2IP, GCNT2, LINGO2, MAD1L1, MAK, MBTPS1, MSRA Arhgap26

CDC16, CDH12, CDH13, CHRNA7, CLK4, COL23A1, COLEC12, CSMD1, CTSL1, CYP2E1, DAPK1, DEPDC6, DGKB, DOC2A, FAT1, GPC5, GPR115, GPR158, HCP5, IMMP2L, KCNT2, KDM4C, KIF25, MAD1L1, MAMDC2, MICA Cd36

ALG10, AMBRA1, APBA1, ATXN7L1, BTBD3, C14ORF156, C14ORF166B, C3ORF20, CD209, CD36, CDH12, CDH13, CFH, CHRNA7, CLK4, COL23A1, COLEC12, CRYZ, CSMD1, CTH, CTSL1, CUL2, CYFIP1, CYP2E1, DAPK1, DEPDC6, DGKB, DOCK4, DOCK8, DSCC1, E2F1, EGFR, F2, FAT1, FOLH1, GDAP1, GLDC, GPC5, GPR115, GPR158, GPRC5C, GRIP2, HYDIN, IMMP2L, JPH1, KCNE1, KCNE2, KCTD13, KDM4C, KIF25, LILRA3, LILRB5, LYN, MAD1L1, MAMDC2, MLLT4, MPHOSPH6, MTHFD1L, MTNR1A, NEGR1, NIPA1, NLRP12, NXN, NXPH1, PARK2, PDE8A, PDE8B, PDPN, PLEKHG1, PPP1R9B, PRKCG, PSMF1 Aldoa

ACACB, ALKBH1, AMBRA1, APBA1, ATP10A, BTBD3, C14ORF166B, CD36, CDC16, CDH12, CDH13, CHRNA7, CLK4, COL23A1, COLEC12, CSMD1, CTH, CTSL1, CYP2E1, DAPK1, DEPDC6, DGKB, DOC2A, F2, FAT1, GLDC, GPC5, GPR115, GPR158, HCP5, IMMP2L, KCNT2, KDM4C, KIF25, MAD1L1, MAMDC2, MAZ Abca4

Genome-wide survey implicates the influence of copy number variants (CNVs) in the development of early-onset bipolar disorder

Keywords: bipolar disorder, copy number variant, genome-wide burden, early age-at-onset, association ARHGAP26, CHRNA9, CNTN4, EEF1D, ERBB2IP, KLRG1, LINGO2, MAD1L1, MSRA, NAALADL2 ABCA4, ACACB, ALKBH1, AMBRA1, APBA1, ATP10A, BTBD3, C14ORF166B,

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction:

Bipolar disorder (BD) is a common, severe mood disorder that is characterized by recurring episodes of extreme exaltation (mania) and depression. Mood symptoms are often accompanied by disturbances in thinking and behavior. BD has a lifetime risk in the general population of 0.5-1.5%. [START_REF] Craddock | Genetics of bipolar disorder[END_REF] Twin, family, and adoption studies have provided strong evidence for a genetic predisposition to BD. [START_REF] Berrettiniw | Progress and pitfalls: bipolar molecular linkage studies[END_REF][START_REF] Kendler | A pilot Swedish twin study of affective illness, including hospital-and population-ascertained subsamples[END_REF] Heritability has been estimated to be between 60% and 85%. [START_REF] Burmeister | Psychiatric genetics: progress amid controversy[END_REF] Changes in the copy number of submicroscopic chromosomal segments, known as copy number variants (CNVs), are a major component of the difference between human genomes. [START_REF] Sebat | Major changes in our DNA lead to major changes in our thinking[END_REF] Based on their frequency, rare (<1%) and common (>1%) CNVs can be distinguished. Several recent studies have shown a strong influence of rare CNVs on the development of neuropsychiatric phenotypes such as autism [START_REF] Sebat | Strong association of de novo copy number mutations with autism[END_REF][START_REF] Weiss | Association between microdeletion and microduplication at 16p11.2 and autism[END_REF] and schizophrenia. [START_REF] Stefansson | Large recurrent microdeletions associated with schizophrenia[END_REF][9][10] To date, only a limited number of studies have been published for BD. [START_REF] Lachman | Increase in GSK3beta gene copy number variation in bipolar disorder[END_REF][START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF][START_REF] Yang | Genomic landscape of a three-generation pedigree segregating affective disorder[END_REF][START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF][START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] Lachman et al. investigated a mixed cohort of Caucasian patients and controls from the Czech Republic and the United States, and found that microdeletions and microduplications, affecting the gene glycogen synthase kinase 3 beta (GSK3B) were significantly increased in patients. [START_REF] Lachman | Increase in GSK3beta gene copy number variation in bipolar disorder[END_REF] Zhang et al. investigated singleton microdeletions (i.e. those occurring only once in the total dataset of patients and controls) of >100 kb in their European American sample of 1 001 BD patients and 1 034 controls, and found that they were overrepresented in patients. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] Recently, Yang et al. published a study of a three-generation Older Amish pedigree with segregating affective disorder. [START_REF] Yang | Genomic landscape of a three-generation pedigree segregating affective disorder[END_REF] They reported that a set of four CNVs on chromosomes 6q27, 9q21, 12p13, and 15q11 were enriched in affected family members, and that these altered the expression of neuronal genes.

Grozeva et al. [START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF] screened a sample of 1 868 patients with BD and 2 938 controls for large (>100 kb) and rare (found in <1% of the population) CNVs. No specific CNV was associated with BD, and the authors found no increased genome-wide burden of CNVs in patients compared to controls. The Wellcome Trust Case Control Consortium (WTCCC [START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] ) investigated common CNVs (found in >5% of the population) of >0.5 kb in a sample of 2007 patients with BD and 3 000 controls, and found no association between CNVs and the disease.

In the present study, we screened the genome-wide SNP data of 882 patients with BD and 872 population-based controls for predicted common and rare CNVs using more than 540 000 autosomal and X-chromosomal markers. All study participants were of German descent. We tested both the overall group of BD patients and a subgroup with an age-at-onset of ≤21 years since several studies have suggested that early-onset BD patients may represent a clinically and genetically more homogeneous subtype of BD. Clinical studies have demonstrated that early-onset BD is a more severe form of the disorder that is characterized by frequent psychotic features, more mixed episodes, greater psychiatric co-morbidity, and poorer response to prophylactic lithium treatment. [START_REF] Weissman | Psychiatric disorders in the relatives of probands with affective disorders. The Yale University-National Institute of Mental Health Collaborative Study[END_REF][START_REF] Schurhoff | Early and late onset bipolar disorders: two different forms of manic-depressive illness?[END_REF][START_REF] Strober | A family study of bipolar I disorder in adolescence. Early onset of symptoms linked to increased familial loading and lithium resistance[END_REF] Familial aggregation is more pronounced in relatives of early-onset BD patients than in relatives of later-onset BD patients. [START_REF] Weissman | Psychiatric disorders in the relatives of probands with affective disorders. The Yale University-National Institute of Mental Health Collaborative Study[END_REF][START_REF] Schurhoff | Early and late onset bipolar disorders: two different forms of manic-depressive illness?[END_REF][START_REF] Strober | A family study of bipolar I disorder in adolescence. Early onset of symptoms linked to increased familial loading and lithium resistance[END_REF][START_REF] Faraone | The genetics of pediatric-onset bipolar disorder[END_REF][START_REF] Grigoroiu-Serbanescu | Different familial transmission patterns in bipolar I disorder with onset before and after age 25[END_REF] Finally, the findings of several studies have suggested the existence of an intra-familial correlation for age-at-onset among bipolar siblings, [START_REF] Baron | Age at onset in bipolar-related major affective illness: clinical and genetic implications[END_REF][START_REF] Leboyer | Age at onset and gender resemblance in bipolar siblings[END_REF] and a segregation analysis has shown that BD is transmitted differentially in early-and later-onset BD families. [START_REF] Grigoroiu-Serbanescu | Different familial transmission patterns in bipolar I disorder with onset before and after age 25[END_REF] To control for the number of technical artifacts, we developed a stringent protocol for quality control (QC) and filtering. On the basis of these data, we conducted statistical tests for the genomewide burden of CNVs and for all specific common and rare CNVs that were found to be associated with BD.

Materials and Methods:

Sample description

Unrelated patients with a clinical history of bipolar disorder (post QC: type I, n=767; type II, n=102; not other specified, n=13) were recruited at two centers: the Central Institute of Mental Health, Mannheim, and the Department of Psychiatry and Psychotherapy of the University of Bonn.

The study was approved by the Institutional Review Boards, and all patients provided written informed consent prior to inclusion. DSM-IV life-time diagnoses of BD were assigned using a consensus best-estimate procedure that was based on all available information including the findings of a structured SCID-I interview. [START_REF] Spitzer | The Structured Clinical Interview for DSM-III-R (SCID) I: History, Rationale, and Description[END_REF] The same set of instruments was used by both centers. [START_REF] Fangerau | Computer-assisted phenotype characterization for genetic research in psychiatry[END_REF] Age-at-onset (AO) was defined as the age at which the first DSM-IV-criteria episode of either depression or mania had occurred. Post QC, the mean age of patients at the time of recruitment was 44.03 years with a standard deviation (SD) of 13.41, the mean age-at-onset was 27.90 years (SD=11.28; median=24; mode=19). The AO distribution of the total sample deviated to the right of the Gaussian distribution (Kolmogorov-Smirnov Z=5.10; P<0.001; positive skewness=1.10). The male/female ratio was 0.47.

Determining the cut-off point for early age-at-onset

Despite extensive debate over the past decade, no consensus has yet been reached concerning the cut-off point for the definition of early and late AO in BD. Authors who have applied the same expectation-maximization algorithm to different samples have described divergent cut-offs for the definition of the early AO. Some have reported that a three-AO-group distribution best fitted their data [START_REF] Bellivier | Age at onset in bipolar I affective disorder: further evidence for three subgroups[END_REF][START_REF] Kennedy | Incidence and distribution of first-episode mania by age: results from a 35-year study[END_REF] and others a two-AO-group distribution. [START_REF] Hamshere | Age-at-onset in bipolar-I disorder: mixture analysis of 1369 cases identifies three distinct clinical sub-groups[END_REF] This demonstrates that the results of an admixture analysis are sample-dependent.

We therefore performed a commingling analysis in the present sample before selecting the cut-off point for the definition of early AO. We used the SEGREG-subroutine of the software S.A.G.E.

(version 6.1, http://darwin.cwru.edu/sage/). [START_REF]Statistical Analysis for Genetic Epidemiology[END_REF] Commingling analysis reveals the distribution mixture of a trait through segregation analysis while allowing for ascertainment correction. Class D models are used as the regressive models in commingling analysis, which assume that the trait under investigation in the study probands is not conditional upon the trait in antecedent family members.

The model that fitted the data best was selected on the basis of the smallest value of the Akaike Information Criterion.

Although the two-AO-group and the three-AO-group models had fitted our data equally well in a preliminary analysis of a larger sample, [START_REF] Grigoroiu-Serbanescu | Two or three age-of-onset groups in bipolar I disorder? Findings of commingling analysis in Romanian and German bipolar I patients[END_REF] the best model in the present sample was a two-AO-group distribution. Prior to QC, the mean AO in the early onset group was 20.67 years (SD=8.40) and the mean AO in the late onset group was 33.20 years (SD=10.30). To select patients with a clear early AO, we only selected patients with an AO that was lower or equal to the rounded mean of the early onset group, i.e. age 21 (AO≤21y). In this AO≤21y subgroup, the mean AO was 17.44 years (SD=2.59).

Prior to QC, our overall sample consisted of 957 patients and 880 controls. Post QC, this was reduced to 882 patients and 872 controls. A total of 291 patients remained in the AO≤21y subgroup and for these patients the mean age-at-recruitment was 38.71 years (SD=12.88), the mean AO was 17.54 years (SD=2.52), and the male/female ratio was 0.43. In the later-onset subgroup with an AO>21 years (AO>21y, n=591), the mean age-at-recruitment was 46.84 years (SD=12.87), and the mean AO was 33.20 years (SD=10.30), and the male/female ratio was 0.49. Controls were drawn from two population-based epidemiological studies: (a) Population-based Recruitment of Patients and Controls for the Analysis (PopGen, n=497) [START_REF] Krawczak | PopGen: population-based recruitment of patients and controls for the analysis of complex genotypephenotype relationships[END_REF] from Schleswig-Holstein (Northern Germany); and (b) the Heinz Nixdorf Recall study (Risk Factors, Evaluation of Coronary Calcification, and Lifestyle; HNR, n=383) [START_REF] Schmermund | Assessment of clinically silent atherosclerotic disease and established and novel risk factors for predicting myocardial infarction and cardiac death in healthy middle-aged subjects: rationale and design of the Heinz Nixdorf RECALL Study. Risk Factors, Evaluation of Coronary Calcium and Lifestyle[END_REF] from Essen, Bochum, and Mülheim a. d. Ruhr (Ruhr area). Post QC, the mean age at recruitment of the control group was 47.98 (SD=11.42), and the male/female ratio was 0.51. All patients and controls reported that they were of German descent.

DNA extraction, genotyping, and quality control

Venous blood samples were collected from all patients and controls. Lymphocyte DNA was extracted either by salting out with saturated sodium chloride solution [START_REF] Miller | A simple salting out procedure for extracting DNA from human nucleated cells[END_REF] or by a Chemagic Magnetic Separation Module I (Chemagen, Baesweiler, Germany) according to the manufacturer's recommendations.

Individuals were genotyped using Illumina's HumanHap550v3 (HH550) or Human610-Quadv1 (H610Q) BeadArrays. The genotype data had been generated as part of a genome-wide association study of BD (Cichon, Mühleisen et al., unpublished data). The H610Q chip contains approximately 60 000 more probes and SNPs than the HH550 array. The majority of the excess content represents non-polymorphic CNV probes, which leads to an excess of CNV calls for individuals genotyped on the H610Q chip (data not shown). To avoid such a bias, we only analyzed SNPs that are present on both chips, i.e. a total of 541 524 SNPs (post QC).

Prior to computational CNV prediction, stringent QC criteria were applied to the genotype data at both the marker and the individual level. SNPs with a call rate of <98% were excluded. Individuals were excluded for the following reasons: (a) DNA call rate <97% (20 patients);(b) differences between X-chromosomally inferred and phenotypic sex (six patients); (c) DNA sample doublets identified by identity-by-state estimates (defined as IBS=2.0, two controls); (d) relatedness of individuals (1.6<=IBS<2.0, no individual excluded); and (e) population outlier according to multidimensional scaling with HapMap phase 2 (one patient and three controls were excluded prior to the present CNV study).

CNV detection

CNVs

were predicted using the program QuantiSNP (version 1.1, http://www.well.ox.ac.uk/QuantiSNP). [START_REF] Colella | QuantiSNP: an Objective Bayes Hidden-Markov Model to detect and accurately map copy number variation using SNP genotyping data[END_REF] The algorithm implemented in QuantiSNP uses an Objective Bayes Hidden-Markov Model to estimate the copy number. To evaluate the presence of a CNV, QuantiSNP uses the normalized intensity data (i.e. log R ratio) and allele frequency data (i.e.

B allele frequency) of each SNP. Both values were calculated by Illumina's BeadStudio Genotyping module (version 3.3.7, http://www.illumina.com/pages.ilmn?ID=169). Individuals were excluded if their SD from the log R ratio or their B allele frequency exceeded certain thresholds: log R ratio >0.36 or B allele frequency >0.12 (22 patients, six controls). We employed the normalization procedure for local GC content implemented in QuantiSNP to improve the accuracy of detection.

The Log Bayes Factor (LBF) was computed for each CNV. This factor indicates the confidence of each predicted CNV, with higher values indicating higher statistical reliability.

To minimize the number of false-positive CNV calls, we only considered CNVs with a LBF ≥30, that spanned a minimum of 30 consecutive SNPs, and which directly affect at least one RefSeq gene. [START_REF] Pruitt | NCBI Reference Sequence (RefSeq): a curated nonredundant sequence database of genomes, transcripts and proteins[END_REF] After applying these filters, we also excluded individuals with more than seven CNVs, since they were extreme outliers in terms of the number of CNV events (27 patients). Following quality control, 1 044 CNV calls from a total of 882 BD patients and 872 controls were statistically analyzed. A total of 291 of the patients (33%) had an AO≤21 years, and 591 patients had an AO>21 years.

To confirm our QuantiSNP-based CNV results, we additionally screened our dataset with PennCNV. [START_REF] Wang | PennCNV: an integrated hidden Markov model designed for high-resolution copy number variation detection in whole-genome SNP genotyping data[END_REF] Both PennCNV and QuantiSNP apply a Hidden-Markov Model to estimate the copy number of an individual. They also take the log R ratio and B allele frequency of each SNP into account, and correct for GC content. The CNV data of both algorithms are available upon request.

Statistical analysis of CNV burden and specific CNVs

All association tests for genome-wide CNV burden and association of specific CNVs were performed using PLINK (version 1.06, http://pngu.mgh.harvard.edu/purcell/plink/). [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analysis[END_REF] We conducted the burden tests for CNV frequency (PROP, RATE) as well as for CNV length (TOTKB, AVGKB),

i.e. the total number of CNVs in patients vs. controls (RATE); proportion of individuals with one or more CNVs in patients vs. controls (PROP); total length spanned by CNVs per individual in the patient group vs. the control group (TOTKB), and average size of CNVs per individual in the patient group vs. the control group (AVGKB). All P-values were generated using 50 000 permutations.

We defined three major comparison groups for the genome-wide burden analyses:

1. All patients vs. all controls 2. AO≤21y patients vs. controls

AO>21y patients vs. controls

We tested each of these three comparison groups for six different categories of CNVs:

1. All CNVs We performed a total of four different burden tests (RATE, PROP, TOTKB, and AVGKB) in 18 test groups, i.e. three major comparison groups for the six different categories of CNVs as described above, resulting in a total of 72 tests for association between CNV burden and BD. To account for all 72 tests, we also applied Bonferroni's method, although this procedure may be too conservative given that the tests were not independent of each other.

In addition, we monitored the distribution of CNVs in chromosomal regions 1q21, 2p16, 7q34-36, 15q11, 15q13, 16p11, 17p12, and 22q11 which have previously been reported to be associated with a variety of neuropsychiatric disorders (Table 2). Since the borders of these regions are known, we relaxed our filter criteria, i.e. we included all CNVs with LBF≥10 which were visually inspected by two independent investigators regardless of the number of affected SNPs. Association tests for these CNVs were performed using Fisher's exact test.

Verification of specific CNVs

All CNVs identified by QuantiSNP and PennCNV that had been found to be associated with neuropsychiatric disorders in previous studies, or which were located within chromosomal regions associated with BD, were visually inspected using Illumina's GenomeStudio.

The specific CNVs that were found to be associated with BD in the present study (6q27 and 10q11)

were verified by quantitative real-time PCR (qPCR) using TaqMan Copy Number Assays (Applied Biosystem, Foster City, CA, USA). We confirmed each CNV carrier by qPCR and also tested non-CNV carriers (as defined by QuantiSNP and PennCNV) to detect possible CNV carriers who had not been identified by QuantiSNP and PennCNV. The status of all CNV carriers was confirmed by qPCR, and no CNV carriers were detected among the putative non-CNV carriers tested. Copy numbers were calculated using the ΔΔCt method implemented in the CopyCaller Software (v1.0, http://www.appliedbiosystems.com/support/software/copycaller/).

Analysis of pathways and biological processes

We analyzed whether genes affected by CNVs were enriched in certain pathways or biological processes using the web-based the program Ingenuity Pathways Analysis platform (IPA, version 8.0, http://www.ingenuity.com). IPA is based on functional annotation and molecular interactions.

Gene lists were assembled using RefSeq genes that are affected by the CNVs identified in the burden analysis (microduplications, singleton deletions). Lists were uploaded into IPA and investigated using the "core analysis" function and default settings. In the functional analysis, biological functions were grouped into different categories from the Ingenuity Knowledge Base. To calculate the statistical significance of pathways and biological processes assigned to gene sets, Pvalues of the Fisher's exact test were corrected by the Benjamini-Hochberg method.

Results:

We systematically analyzed our samples for significant genome-wide differences in the distribution of all CNVs between patients and controls (genome-wide burden tests) as well as for a significant overrepresentation of specific CNVs in patients or controls.

General description of the CNV dataset

Following the baseline QC of SNP data, QuantiSNP identified a total of 124 146 putative CNVs in the initial sample of 957 patients and 880 controls (an average of 67.6 CNVs per individual).

Following the application of all QC filters, 1 044 potential CNVs remained in the filtered sample of 882 BD patients and 872 controls (an average of 0.59 CNVs per individual).

We examined the distribution of the number of CNVs per individual and identified a total of 27 extreme outliers in terms of CNV observations, with more than seven CNVs being detected in each individual. The samples of a total of 24 of these individuals were clustered in the outer rows of the same 96-well plate, and thus these findings are likely to represent plate effects. These individuals were excluded from the downstream analyses.

Association analysis for genome-wide CNV burden

Overall, 10 of the 72 genome-wide burden tests revealed nominally significant differences in CNV burden between patients and controls. In the following, we provide a detailed description of the most important findings, as outlined in Table 1.

All patients vs. controls. In the total sample of 882 BD patients and 872 controls, the genome-wide burden of singleton microdeletions showed nominally significant association for two out of the four tests performed. The average total length of all singleton microdeletions per individual was 487. the most significant of all of the burden analyses, withstood correction for multiple testing with the Bonferroni method which accounted for the number of all tests performed in this study (n=72, P adjusted =0.029). However, there were no significant differences in the proportion of individuals who carried either microdeletions or singleton microdeletions.

Patients with AO>21y vs. controls. In the third test group, we analyzed the burden of CNVs in the AO>21y subgroup (n=591) vs. all controls. These tests revealed no significant differences in CNV burden between patients and controls.

Association analyses of specific CNVs

We identified two common microduplications that were significantly overrepresented in patients compared to controls: (a) a 248 kb microduplication on chromosome 6q27, and (b) a 160 kb microduplication on chromosome 10q11 (Figure 1).

The 10q11 microduplication (Figure 1a) was observed in 53 patients (6.01%) and in 32 controls There were no significant differences in distribution when all BD patients were compared to controls. In CNV carriers, slight differences in CNV size were observed, with a shared overlap of around 248 kb (chr6:168.09-168.33 Mb, NCBI build 36) that was due to 110 adjacent SNPs. Three genes are affected by this microduplication: (a) kinesin family member 25 (KIF25), (b) FERM domain containing 1 (FRMD1), and (c) parts of the 3' terminus of mixed-lineage leukemia translocated to 4 (MLLT4). The gene dapper, antagonist of beta-catenin, homolog 2 (DACT2) lies around 115 kb downstream from this common microduplication. DACT2 participates in the WNT signaling pathway [START_REF] Machon | A dynamic gradient of Wnt signaling controls initiation of neurogenesis in the mammalian cortex and cellular specification in the hippocampus[END_REF][START_REF] Toledo | Wnt signaling in neuroprotection and stem cell differentiation[END_REF] which is known to regulate neurogenesis and neuroprotection.

(
In a follow-up analysis, we tested whether patients carrying either the common CNVs on 6q27 or the CNV on 10q11 showed differences in sex distribution or family history of psychiatric disorder compared to patients who were non-carriers. No significant associations were observed for either phenotypic item (data not shown).

Specific CNVs at loci previously associated with psychiatric disorder

We tested whether CNVs in one of six genomic regions that have previously been reported to be associated with neuropsychiatric disorders (1q21, 2p16, 7q34-36, 15q11, 15q13, 16p11, 17p12 and 22q11) were overrepresented in our BD patients. CNVs in these regions were not significantly overrepresented (Table 2). However, the power of our sample to find significant association with these rare CNVs was low.

Pathways and biological processes impacted by CNVs in early-onset patients

To further characterize the two top association findings of our burden analyses in AO≤21y patients (Table 1), we used IPA to explore possible functional relationships between genes covered by either microduplications or singleton microdeletions. In the analysis of the 46 genes hit by singleton microdeletions, the significant first five top hits for biological processes (i.e. "drug metabolism", "lipid metabolism", "molecular transport", "small molecular biochemistry" and "endocrine system development and function") were enriched by the presence of the genes SLCO1B1 and SLCO1B3, which are both located on 12p12. These were covered by a singleton microdeletion in one AO≤21y

patient. These categories were therefore omitted from data interpretation. Further top process categories that were significantly overrepresented due to several input genes being hit by singleton microdeletions included "endocrine system disorder", "genetic disorder", "metabolic disease", "immunological disease," and "infectious disease" (Table 3). The same functional categories were found to be enriched when analyzing the 287 genes affected by microduplications in AO≤21y patients (Table 4) although there was no overlap in affected genes between the two gene lists, except for MAD1L1.

No significant result was obtained for canonical pathways in IPA after correction for multiple testing (data not shown). This was probably due to the limited number of genes introduced into the analysis.

Discussion:

In the present study, we investigated a large sample of patients and controls of German descent for the presence of CNVs that may be involved in the development of BD. Our analyses included tests to monitor the genome-wide burden of CNVs (frequency and size) between patients and controls as well as tests to detect specific common and rare CNVs that were significantly overrepresented in either patients or controls. Since clinical and formal genetic studies have suggested that BD with an AO≤21y may be genetically distinct from BD with an AO >21 years, [START_REF] Weissman | Psychiatric disorders in the relatives of probands with affective disorders. The Yale University-National Institute of Mental Health Collaborative Study[END_REF][START_REF] Schurhoff | Early and late onset bipolar disorders: two different forms of manic-depressive illness?[END_REF][START_REF] Strober | A family study of bipolar I disorder in adolescence. Early onset of symptoms linked to increased familial loading and lithium resistance[END_REF][START_REF] Faraone | The genetics of pediatric-onset bipolar disorder[END_REF][START_REF] Grigoroiu-Serbanescu | Different familial transmission patterns in bipolar I disorder with onset before and after age 25[END_REF][START_REF] Baron | Age at onset in bipolar-related major affective illness: clinical and genetic implications[END_REF][START_REF] Leboyer | Age at onset and gender resemblance in bipolar siblings[END_REF] we also analyzed these subgroups separately. The basis for the statistical analyses was provided by high-quality CNV prediction (QuantiSNP) using the SNP intensity data of 1 754 individuals and verification (PennCNV, visual inspection, and TaqMan for specific CNVs). Parameters were specified to ensure that only relatively large CNVs (detected by ≥30 consecutive SNPs) with high statistical confidence (QuantiSNP: log Bayes factor ≥30; PennCNV: confidence value ≥30) passed QC. Using such stringent criteria is always a trade-off between sensitivity and quality of data. We have chosen our criteria to reduce type I error and are aware that this may lead to an increase of type II error at the same time.

Our genome-wide burden tests (RATE and PROP) provided no evidence that the overall number of singleton CNVs (both microdeletions and microduplications) overlapping with RefSeq genes was enriched in the total sample of BD patients compared with controls. The TOTKB and AVGKB tests demonstrated that the length of singleton microdeletions differed significantly between patients and controls (Table 1). Singleton microdeletions in patients were approximately twice the size (490 kb) of those in controls (260 kb; AVGKB P=0.014). Separate analyses of the AO≤21y subgroup and the AO>21y subgroup demonstrated that this effect was mainly attributable to the AO≤21y subgroup, in which the average size of microdeletions was approximately 680 kb (AVGKB P=0.0056), and to a lesser extent to the AO>21y subgroup, in which the average size of singleton microdeletions (348 kb) was not significantly larger than in controls (AVGKB P=0.17; data not shown). The most significant finding of our burden analyses was that the total number of patients with at least one microduplication (both common as well as singleton microduplications) was significantly higher in the AO≤21y group, but not in the AO>21y subgroup (PROP, P=0.0004). The burden results suggest that both a higher CNV load and a larger CNV size are associated with BD in patients with an AO≤21y. The effect in the early-onset group is was attributable to longer singleton deletions as well as to a higher frequency of microduplications (common and rare). Analysis of either the overall sample or the AO>21y sample alone produced only marginal evidence that the genome-wide burden of CNVs plays a role in disease development.

These results show a similar trend to those of Zhang et al. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] who reported a significant overrepresentation of singleton deletions of more than 100 kb in their BD cases (16.2%) compared to controls (12.3%; P=0.007). Interestingly, they also found that this effect was more pronounced in an early AO form of BD (age of mania onset ≤18 years). This is consistent with our observation of the presence of longer singleton deletions in early-onset BD patients. However, Zhang et al. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] found no evidence of a higher frequency of microduplications in their early-onset mania subsample.

Unfortunately, a more exact comparison of our data with that of Zhang et al. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] is limited by the methodological differences between the two studies: Zhang et al. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] used Affymetrix 6.0 arrays whereas we used Illumina HH550/H610Q arrays; Zhang et al. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] performed CNV detection with Birdsuite whereas we used QuantiSNP and PennCNV; Zhang et al. [START_REF] Zhang | Singleton deletions throughout the genome increase risk of bipolar disorder[END_REF] took all CNVs above a certain size threshold into account whereas we filtered for CNVs that overlapped with at least 30 consecutive SNPs and RefSeq genes.

Another recent study by Yang et al. [START_REF] Yang | Genomic landscape of a three-generation pedigree segregating affective disorder[END_REF] found no evidence for any significant association between the average number and size of CNVs and affective disorders (including BD and major depression).

However, their study investigated a single three-generation Old Order Amish family, and their main focus was upon the identification of CNVs that co-segregated with disease status across generations. The limited number of affected (n=19) and unaffected family members (n=32) and their genetic relatedness, as well as the broader phenotype definition used may have resulted in limited power to investigate the genome-wide burden aspects of CNVs in BD.

Two recent genome-wide studies of CNVs in BD have been published (Grozeva et al. [START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF] , WTCCC [START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] ), which investigated common and rare CNVs. They tested for association between BD and specific CNVs and CNVs at previously reported loci, as well as for genome-wide burden.

Neither of these two studies found evidence for the involvement of CNVs in disease development.

The majority of patients and controls investigated by Grozeva et al. [START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF] were also analyzed by the WTCCC [START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] , and thus the results of the two studies are not entirely independent of each other. Major methodological differences exist between these two studies and our study, which hamper any direct comparison of the results. Firstly, different arrays were used to screen for CNVs (Grozeva et al. [START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF] :

Affymetrix' Genechip Human Mapping 500K Array Set; WTCCC [START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] : Agilent Comparative Genomic Hybridization arrays). Secondly, Grozeva et al. [START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF] searched for rare CNVs (MAF <1%) only, whereas the WTCCC [START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] investigated common CNVs (MAF >5%) only. The present study did not have any restrictions with regard to CNV frequency, and the frequencies of the two specific BDassociated CNVs on 10q11 (3.68% in controls) and 6q27 (2.53% in controls) are within a frequency range that was not investigated by Grozeva et al. [START_REF] Grozeva | Rare Copy Number Variants A Point of Rarity in Genetic Risk for Bipolar Disorder and Schizophrenia[END_REF] or the WTCCC. [START_REF]Genome-wide association study of CNVs in 16,000 cases of eight common diseases and 3,000 shared controls[END_REF] Another important difference is the separate analyses of AO subgroups in the present study, which suggest that CNVs play a role in early onset, but not in later onset BD. Since the aforementioned studies did not take AO into account, it is unclear whether such an effect was present in their BD samples. Nonetheless, all studies, including the present study, are in agreement that CNVs do not appear to influence BD when AO is not taken into account.

In a subsequent step of our genome-wide burden analysis, we performed an IPA to investigate whether the genes affected by microduplications (n=287) or singleton microdeletions (n=46) in AO≤21y patients were significantly enriched for biological processes or pathways. When removing single gene-based enriched categories, the top five significantly overrepresented biological function categories were the same for both the microduplications and the microdeletions gene lists, i.e. the disorder and disease processes "endocrine system disorder", "genetic disorder", "metabolic disease", "immunological disease," and "infectious disease". Follow-up studies in independent samples are clearly necessary to confirm that biological functions within these disease categories are involved in the pathophysiology of (early-onset) BD. One interesting gene with prior suggestive evidence for association with BD, mitotic arrest deficient-like 1 (MAD1L1, 7p22.3), was affected by a singleton microdeletion as well as a singleton microduplication in early-onset BD patients. Support for the involvement of this locus was provided by a recent genome-wide association study (GWAS) of BD, in which two MAD1L1 SNPs (rs11764590 and rs10278591; r 2 =0.7) were the second and third best results in the meta-analysis step (P=1.28x10 -7 and P=1.81x10 -5 , respectively, unpublished data).

MAD1L1 is a mitotic spindle assembly checkpoint protein. Homozygous knockout of MAD1L1 in mice confers embryonic lethality, indicating that MAD1L1 plays an essential role during embryonic development. [START_REF] Iwanaga | Heterozygous deletion of mitotic arrest-deficient protein 1 (MAD1) increases the incidence of tumors in mice[END_REF] Two genes affected by microduplications in early AO patients, epidermal growth factor receptor (EGFR, 7p11.2, one patient) and nucleoredoxin (NXN, 17p13.3, one patient, two controls), were identified as susceptibility genes for BD by two recent GWAS. [START_REF] Sklar | Whole-genome association study of bipolar disorder[END_REF][START_REF] Baum | A genome-wide association study implicates diacylglycerol kinase eta (DGKH) and several other genes in the etiology of bipolar disorder[END_REF] EGFR rs17172438 (P=3.26x10 - 5 , OR=1.32) and EGFR rs729969 (P=3.30x10 -5 , OR=1.36) were ranked among the top 20 findings in the GWAS by Sklar et al. [START_REF] Sklar | Whole-genome association study of bipolar disorder[END_REF] EGFR and its ligands are cell signaling molecules that mediate diverse downstream cellular functions, including cell proliferation and differentiation. In a further GWAS, Baum et al. [START_REF] Baum | A genome-wide association study implicates diacylglycerol kinase eta (DGKH) and several other genes in the etiology of bipolar disorder[END_REF] found that NXN rs2360111 (P=0.0003, OR=1.23) were associated with BD.

The NXN gene encodes the protein nucleoredoxin which is involved in the inhibition of the WNT signaling pathway. [START_REF] Funatoy | The thioredoxinrelated redox-regulating protein nucleoredoxin inhibits Wnt-β-catenin signalling through dishevelled[END_REF] Another aim of the present study was to identify specific BD-associated CNVs. We found that the frequency of two relatively common CNVs on 6q27 and 10q11 differed significantly between patients and controls. Common microduplications located on chromosome 10q11 showed association in the overall sample (P=0.035, OR=1.53). Again, this effect was much stronger in the AO≤21y subgroup (P=0.00052, OR=2.79). The other microduplication on chromosome 6q27 was overrepresented in the AO≤?21y subgroup (P=0.0039, OR=2.40), but not in the overall sample.

Interestingly, the 6q27 CNV region was implicated in BD in a three-generation Old Order Amish pedigree. [START_REF] Yang | Genomic landscape of a three-generation pedigree segregating affective disorder[END_REF] The 6q27 region was one of four regions that were enriched in affected pedigree members and which had an effect on the expression of genes within or near the rearrangement. The identification of the 6q27 region in two independent studies is clearly of interest, although further independent studies are required to support this finding. One limitation which prevents stronger conclusions being drawn from the Yang et al. study is that there was only modest co-segregation of the 6q27 CNV with BD. Thus the possibility that that this relatively common variant segregated within the family independently of disease cannot be excluded.

In summary, the present study found evidence for a significant association between BD and microdeletions and microduplications. Although the frequency of microdeletions was not significantly higher in patients compared to controls, the size of singleton microdeletions was significantly larger. Our data also suggest that a higher frequency of microduplications is implicated in disease development. We found both the genome-wide burden of microduplications as well as common specific microduplications on 6q27 and 10q11 to be enriched in BD patients. A further important finding of our study is that CNVs were strongly associated in patients with an AO ≤21 years, but not in patients with an AO >21 years. In the overall BD sample, only a very weak association with CNVs was detected. Our results support for findings from previous clinical and formal genetic studies that early and later onset BD may represent genetically distinct forms of the disease. Future studies of CNVs in BD should therefore take the AO of their patient samples into account.

  P=0.0056) and spanned longer chromosomal regions per individual than in controls (679.6 kb in AO≤21y patients vs. 261.2 kb in controls, TOTKB: P=0.0084). Furthermore, we observed that the total proportion of individuals with at least

	one microduplication (44.3% in AO≤21y patients vs. 33.1% in controls, PROP: P=0.00040), at least
	one CNV (52.9% in AO≤21y patients vs. 42.3% in controls, PROP: P=0.00092), or at least one
	singleton microduplication (17.2% in AO≤21y patients vs. 11.9% in controls, PROP: P=0.017) was
	significantly higher in this BD subgroup. The PROP test P-value for microduplications, which was

6 kb in patients compared to 265.1 kb in controls (TOTKB: P=0.014). The average size per singleton microdeletion was 472.6 kb in patients and 249.3 kb in controls (AVGKB: P=0.014). The PROPand RATE-tests generated no significant P-values.

Patients with an early AO vs. controls. When comparing the 291 AO≤21y patients with all controls, we again found that singleton microdeletions in patients were, on average, larger (661.7 kb in AO≤21y patients vs. 249.3kb in controls, AVGKB:

Table 2

 2 Association results for specific CNVs in BD at loci with prior evidence for association with psychiatric disorders.Abbreviations: n, number of observed CNVs; Microdel, microdeletions of variable size; Microdup, microduplications of variables size, ASD, autism spectrum disorder; BD, bipolar disorder; EP, epilepsy; MR, mental retardation; SCZ, schizophrenia; P, P-value from Fisher's exact test; CI, confidence interval.

		Previous studies				Present study
						CNV (frequency)
	Locus	Type (n)	Associated disorder	Refere nce	Type (n)	In 882 patients In 872 controls	P (OR, 95%CI)
	1q21	Microd el	ASD, MR, SCZ	8,43	Microdel (1)	0 CNV (0%)	1 Microdel (0.11%)	0.49 (0.99 [0.01-77.56])
	2p16	Microd el	ASD, SCZ	44,45	Microdel (8)	5 Microdel (0.56%)	3 Microdel (0.34%)	0.73 (1.65 [0.32-10.67])
					Micordel (15),	11 Microdel	4 Microdel
	7q34-36	Microd el	SCZ	46	Microdu p(4)	(1.25%), 1 (0.11%) Microdup	(0.45%). 3 (0.34%) Microdup	0.36 (1.70 [0.62-5.13])
	15q11	Microd el	SCZ	8	Microdel (4), Microdu p (8)	1 Microdel (0.11%), 6 Microdup (0.56%)	3 Microdel (0.34%), (0.57%) 2 Microdup	0.77 (1.39 [0.38-5.56])
					Microdel	15 Microdel	20 Microdel
	15q13	Microd el	ASD, EP, MR, SCZ	8,47,48	(35), Microdu	(1.25%), 19 Microdup	(2.29%), 12 Microdup	0.90 (1.05 [0.62-1.78])
					p (31)	(1.70%)	(1.38%)
		Microd			Microdel	0 Microdel
	16p11	el, Microd	ASD	7	(0), Microdu	(0%), Microdup 3	0 CNV (0%)	0.25 (2.97 [0.24, 156.29])
		up			p (3)	(0,23%)
					Microdel	1 Microdel	0 Microdel
	17p12	Microd el	SCZ	49	(1), Microdu	(0.11%), 4 Microdup	(0%), Microdup 4	1 (1.24 [0.27-6.26])
					p (8)	(0.45%)	(0.45%)
					Microdel	10 Microdel	6 Microdel
	22q11	Microd el	ASD, BD, MR, SCZ	9,50,51	(16), Microdu	(1.13%), 18 Microdup	(0.69%), 18 Microdup	0.67 (1.16 [0.64-2.11])
					p (36)	(2.04%)	(2.18%)

Table 4

 4 Functional categories a significantly enriched by several genes affected by microduplications in the AO≤21y subgroup.

	Category	Genes in CNV	P-value a	Adjusted P-value b
	Immunological	ACACB, ALKBH1, ATP10A, CD209, CDC16,	1.01 x 10 -6 -2.65 x	1.34 x 10 -3 -1.38 x
	Disease	CDH12, CDH13, CLEC4M, CSMD1, CYFIP1,	10 -2	10 -1
		DEPDC6, DGKB, DOC2A, FAT1, GPC5,
		GPR158, HCP5, KDM4C, LYN, MAD1L1,
		MAMDC2, MICA, MTHFD1L, NEGR1, NXN,
		NXPH1, PARK2, PDE8A, PDE8B, RBP3,
		RSPO4, TAOK2
	Infectious	ACACB, CD209, CD36, CDH12, CHRNA7,	1.01 x 10 -6 -2.65 x	1.34 x 10 -3 -1.38 x
	Disease	CLEC4G, CLEC4M, CSMD1, CTSL1, CYFIP1,	10 -2	10 -1
		DEPDC6, DGKB, EGFR, F2, FAT1, GPC5,
		HCP5, KDM4C, LILRA3, LILRB5, LYN,
		MAD1L1, MPHOSPH6, NXN, PARK2, PDE8A,
		RSPO4, SEC61G, SNW1, STXBP2, TOP3B,
		TRIM58, XAB2, ZNF254, ZNF354A

Acknowledgements:

We thank all of the patients who participated in this study. We also thank all of the probands from the community-based cohort of PopGen as well as that of the Heinz Nixdorf Recall (HNR) study, which was established with the support of the Heinz Nixdorf Foundation. The present study was supported by the German Federal Ministry of Education and Research (BMBF) within the context of the National Genome Research Network plus (NGFNplus) and the MooDS-Net (grant 01GS08144 to S.C. and M.M.N., grant 01GS08147 to M.R.). M.M.N. also received support from the Alfried Krupp von Bohlen und Halbach-Stiftung. Some of the results of this study were obtained using the program package S.A.G.E., which is supported by a U.S. Public Health Service Resource Grant (RR03655) from the National Center for Research Resources.

We are grateful to J. Sebat (Department of Psychiatry, University of California, San Diego) for his critical review of the manuscript.

Conflict of Interest statement:

The authors declare they have no competing financial interests, as defined by Nature Publishing Group, or any other interests that might be perceived to influence the results and discussion reported in this paper.

Tables:

Table 1 Association results for genome-wide burden of CNVs in BD. CNVs in patients and controls were sorted according to their length. Schematic drawings were generated using UCSC Genome browser (NCBI Build 36.1).