Tomáš Flouri

Costas S Iliopoulos

Kunsoo Park

Solon P Pissis

GapMis-OMP: pairwise short-read alignment on multi-core architectures

Pairwise sequence alignment has received a new motivation due to the advent of next-generation sequencing technologies, particularly so for the application of re-sequencing-the assembly of a genome directed by a reference sequence. After the fast alignment between a factor of the reference sequence and a high-quality fragment of a short read by a short-read alignment programme, an important problem is to find the alignment between a relatively short succeeding factor of the reference sequence and the remaining low-quality fragment of the read allowing a number of mismatches and the insertion of a single gap in the alignment. In this article, we present GapMis-OMP, a tool for pairwise short-read alignment that works on multi-core architectures. It is designed to compute the alignments between all the sequences in a first set of sequences and all those from a second one in parallel. The presented experimental results demonstrate that GapMis-OMP is more efficient than most popular tools.

Introduction

A gap is a consecutive sequence of holes in an alignment. When applying alignment algorithms to biological sequence data, it is sometimes desirable to globally penalise the formation of long gaps instead of locally penalising the individual deletion or insertion of letters. A gap in a biological sequence can be regarded as the absence (resp. presence) of a region, which is (resp. is not) present in another sequence, because of the natural diversity among individuals. The gap concept in sequence alignments is therefore important in many biological applications because the insertion or deletion of an entire region (particularly in DNA) often occurs as a single event. Many of these single mutational events can create gaps of varying sizes with almost equal likelihood-within a wide, but bounded, range of sizes.

Our work is motivated by the well-known and challenging problem of re-sequencing-the assembly of a genome that is directed by using a reference genome. Continuous advances in sequencing technology are turning whole-genome sequencing into a routine procedure, resulting in massive amounts of sequence data that need to be processed. Tens of gigabytes of data, in the form of short sequences (reads), need to be mapped (aligned) back to reference sequences, a few gigabases long, to infer the read from which the genomic location derived. This is a challenging task because of the high data volume and the large genome sizes. In addition, the performance, in terms of speed, sensitivity, and accuracy, deteriorates in the presence of inherent genomic variability and sequencing errors, particularly so, for relatively short consecutive sequences of deletions or insertions in the reference sequence or in the reads.

A broad variety of short-read alignment programmes, (e.g. Bowtie [START_REF] Langmead | Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[END_REF], SOAP2 [START_REF] Li | SOAP2: an improved ultrafast tool for short read alignment[END_REF], REAL [START_REF] Frousios | REAL: an efficient REad ALigner for next generation sequencing reads[END_REF]) has been released recently to address the task of mapping millions of short reads to a genome, placing emphasis on various aspects of the challenge (e.g. time and memory efficiency, sensitivity, and accuracy). Although all programmes allow for a small number of mismatches in the alignment, some of them either perform poorly when allowing the insertion of gaps or do not allow it at all.

The seed-and-extend strategy is applied in most current short-read alignment programmes [START_REF] Altschul | Basic Local Alignment Search Tool[END_REF]. After a fast alignment between a fragment of the reference sequence and a high-quality fragment of a short read (positions 1-3 in square brackets in Figure 1) by a short-read alignment programme, an important problem is to find the alignment between a relatively short succeeding fragment of the reference sequence and the remaining low-quality fragment of the read (positions 4-9 in Figure 1); allowing a number of mismatches (position 8 in Figure 1) and the insertion of a single gap (positions 5-6 in Figure 1) either in the fragment of the reference sequence or in the read.

One main observation from Figure 1 is that a gap might need to be inserted in the leftmost position of the alignment (position 4 in Figure 1), and that we are unable to know the exact length of the succeeding fragment of the reference sequence to be aligned a priori. Hence, we need an intermediate approach between the global (e.g. Needleman-Wunsch algorithm [START_REF] Saul | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF]) and the local alignment (e.g. Smith-Waterman algorithm [START_REF] Smith | Identification of common molecular subsequences[END_REF]), Concerning the length of the gaps, a broad range of lengths is possible. In practice, however, the length of reads is too small to confidently and directly detect a large gap. In Figure 2, the distribution of lengths of gaps in exome sequencing is demonstrated1 . The shape of the gap length distribution is consistent with other studies (cf. [START_REF] Ng | Targeted capture and massively parallel sequencing of 12 human exomes[END_REF]). The presented data reflect a gap occurrence frequency of approximately 5.7 × 10 -6 across the exome-the part of the genome formed by exons that codes portions of genes in the genome which are expressed.

Reference G C G A C G T C C G A A | | | | | . | G C G A ⋆ ⋆ T A C
The main observation from Figure 2 is the exponential decrease of gaps with increasing length, and a preference for lengths which are multiples of 3 base pairs (bp)-loss or gain of codons. For short reads in the range of 25-150 bp, the presence of multiple gaps is unlikely given the above gap occurrence frequency, and can greatly reduce the mapping confidence (accuracy) of those reads. Hence, applying a traditional dynamic programming approach, which allows for multiple mismatches, deletions, and/or insertions in the alignment, can deteriorate mapping confidence.

In this article, motivated by the aforementioned observations, we consider the problem of approximate string matching with k-mismatches and a single gap. We present GapMis-OMP, a tool for pairwise short-read alignment that works on multi-core architectures. It is based on algorithm GapMis, introduced in [START_REF] Flouri | Approximate string-matching with a single gap for sequence alignment[END_REF], for computing a modified version of the traditional dynamic programming matrix for sequence alignment to solve the above problem. GapMis-OMP is designed to compute the alignments between all the sequences in a first set of sequences and all those from a second one in parallel.

Definitions and notation

An alphabet Σ is a finite non-empty set whose elements are called letters. A string on an alphabet Σ is a finite, possibly empty, sequence of elements of Σ. The zero-letter sequence is called the empty string, and is denoted by ε. The length of a string x is defined as the length of the sequence associated with the string x, and is denoted by |x|. We denote by x[i], for all 1 ≤ i ≤ |x|, the letter at index i of x. Each index i, for all 1 ≤ i ≤ |x|, is a position in x when x = ε. It follows that the ith letter of x is the letter at position i in x, and that

x = x[1 . . |x|]
A string x is a factor of a string y if there exist two strings u and v, such that y = uxv. Let x, y, u, and v be strings, such that y = uxv holds. If

u = ε, then x is a prefix of y. If v = ε, then x is a suffix of y.
Let x be a non-empty string and y be a string. We say that there exists an occurrence of x in y, or, more simply, that x occurs in y, when x is a factor of y. Every occurrence of x can be characterised by a position in y. Thus we say that x occurs at the starting position i in y when y[i . . i + |x| -1] = x. It is sometimes more suitable to consider the ending position i + |x| -1. The Hamming distance δ H for two strings of the same length, is defined as the number of positions where the two strings possess different letters.

A don't care letter is a special letter that does not belong to alphabet Σ, and matches with itself as well as with any letter of Σ. It is denoted by ⋆. A gap is a finite sequence of such don't care letters. A gap string is a finite, possibly empty, sequence of elements of the alphabet Σ ∪ {⋆}. Two letters a and b of alphabet Σ ∪ {⋆} are said to correspond, denoted by a ≈ b, if they are equal, or, if at least one of them is the don't care letter. The G-distance, denoted by δ G , for two gap strings of the same length is defined as the number of positions in which the two strings possess letters that do not correspond. A gap string x is called single-gap string if there exist two strings u and v on alphabet Σ and a gap g, such that x = ugv. Let conc(y ′ ) be an operation that, given a gap string

y ′ = y 0 g 0 y 1 g 1 . . . y n-2 g n-2 y n-1
where y i ∈ Σ * , for all 0 ≤ i < n, and g j ∈ {⋆} * , for all 0 ≤ j < n -1, returns the string y = y 0 y 1 . . . y n-1 , such that y ∈ Σ * .

The approximate string-matching with k-mismatches and a single gap problem can now be formally defined as follows.

Problem 1. Given a text t of length n, a pattern x of length m ≤ n, an integer k, such that 0 ≤ k < m, and integers α and β, such that 0 ≤ α ≤ β and β < n, find all prefixes of t, such that for each prefix y either there exists a single-gap string y ′ , with a gap g, such that y = conc(y ′ ), δ G (x, y ′ ) ≤ k, and α ≤ |g| ≤ β or there exists a single-gap string x ′ , with a gap g, such that x = conc(x ′ ), δ G (x ′ , y) ≤ k, and α ≤ |g| ≤ β or δ H (x, y) ≤ k and α = 0 Example 1. Let t = AGCAGAGGAGCAGGCGTTCCGTGGT, x = ACCGT, k = 2, α = 6, and β = 7. The solution to this problem instance is the set {11, 17, 22} of ending positions. For example, the solution contains 11 since there exists a single-gap string x ′ = ACC⋆ ⋆ ⋆ ⋆ ⋆⋆GT, with a gap g = ⋆ ⋆ ⋆ ⋆ ⋆⋆, such that x = conc(ACC⋆ ⋆ ⋆ ⋆ ⋆⋆GT), δ G (x ′ , t[1 . . 11]) = 2, and |g| = 6.

Let G[0 . . n, 0 . . m] be a matrix, where G[i, j] contains the minimum number of mismatches of the alignment between factor t[1 . . i] of t and factor x[1 . . j] of x allowing the insertion of a single gap either in t[1 . . i] or in x[1 . . j], for all 1 ≤ i ≤ n, 1 ≤ j ≤ m. We say that x[1 . . j] matches t[1 . . i] with at most k-mismatches and a single gap if and only if

G[i, j] ≤ k, for all 1 ≤ j ≤ m, 1 ≤ i ≤ n.
In order to compute the exact location of the inserted gap, either in the text or in the pattern, we also need to compute a matrix H[0 . . n, 0 . . m],

such that

H[i, j] =    j -i if G[i, j] = G[i, i] and i ≤ j i -j if G[i, j] = G[j, j] and i > j 0 otherwise
Example 2. Let t = AGGTCAT, x = GGGTA, and β = 2. Table 1a and Table 1b illustrate matrix G and matrix H, respectively.

0 1 2 3 4 5 ǫ G G G T A 0 ǫ 0 0 0 1 A 0 1 1 1 2 G 0 0 1 1 1 3 G 0 0 1 1 1 4 T 1 1 1 1 5 C 1 1 2 6 A 1 1 7 T 2 (a) Matrix G 0 1 2 3 4 5 ǫ G G G T A 0 ǫ 0 1 2 1 A 1 0 0 0 2 G 2 0 0 0 2 3 G 0 0 0 1 2 4 T 0 0 0 1 5 C 2 1 0 6 A 2 0 7 T 0 (b) Matrix H Table 1
. Matrix G and matrix H for t = AGGTCAT, x = GGGTA, and β = 2

Algorithm GapMis

Algorithm GapMis, introduced in [START_REF] Flouri | Approximate string-matching with a single gap for sequence alignment[END_REF], computes matrices G and H. It takes as input the text t of length n, the pattern x of length m, and the threshold β.

Proposition 1 ([8]

). There exist at most 2β + 1 cells of matrix G that yield a solution for Problem 1.

Proposition 2 ([8]

). Problem 1 can be solved by algorithm GapMis in time O(mβ).

As soon as we have computed matrices G and H, we can apply a simple alignment scoring scheme, depending on the application of the algorithm, to compute the maximum score among all possible alignments of t and x in time Θ(β) by Proposition 1. The total score of each alignment is obtained by adding these two scores, and the optimal alignment is the one with the maximum total score. The same alignment score scheme is applied in package EM-BOSS [START_REF] Rice | EMBOSS: The European Molecular Biology Open Software Suite[END_REF].

The programme takes as input arguments two files with the two sets of sequences (DNA or protein) in FASTA format. It consists of an outer for loop, going through the sequences of the second set (target sequences), and an inner for loop, going through the sequences of the first set (query sequences), to compute the optimal semi-global alignment for each pair. We used the omp parallel for directive to effectively define the inner for loop as a parallel region, as there exists no dependency between those tasks (pairs). The number of threads to be used can be set by the user. The programme produces an EMBOSS-like text file with all the optimal semi-global alignments as output.

Experimental results

To the best of our knowledge, GapMis-OMP is the first tool for pairwise short-read alignment between two sets of sequences. In [START_REF] Flouri | GapMis: a tool for pairwise sequence alignment with a single gap[END_REF], in order to highlight the suitability of algorithm GapMis, we compared it against EMBOSS needle, which implements the Needleman-Wunsch algorithm for global alignment, and EMBOSS water, which implements the Smith-Waterman algorithm for local alignment. At present, these two programmes are two of the most widely used pairwise sequence alignment programmes. We demonstrated that needle and water can, by design, not guarantee the insertion of at most one gap in the alignment.

In this article, in order to evaluate the efficiency of GapMis-OMP, we compared its performance to the respective performance of needle and water. We simulated 10, 000 pairs of 100 bp-long DNA sequences, such that each pair consists of two identical sequences; and then we inserted, in one of the two sequences, a single gap of random length ranging from 1 to 30 and a random number of mismatches ranging from 1 to 10. We repeated the same experiment by simulating 150, 200, and 250 bp-long sequences. Two versions of GapMis-OMP were used; one with the modifier -m 30 to set β = 30 and use the Θ(mβ) algorithm, and one without it to use the Θ(mn) algorithm.

As it is demonstrated by the results in Figure 3 GapMis-OMP was able to complete the assignment much faster than needle and water. The version with the modifier -m 30 was always the fastest confirming our theoretical results.

In order to evaluate the parallel efficiency of GapMis-OMP, we simulated a set of 100, 000 75 bp-long query sequences and 100 100 bp-long target sequences. Six versions of GapMis-OMP were used: with 1 thread (-t 1); with 4 threads (-t 4); with 8 threads (-t 8); with 16 threads (-t 16); with 32 threads (-t 32); and with 48 threads (-t 48). We repeated the same experiment with 150 and 200 bp-long target sequences. As depicted in Figure 4, the -t 4 version is able to complete the assignment

Fig. 1 .

 1 Fig.1. Alignment between the fragment of the reference sequence, starting at position 1 and ending at position 9, and the read with one mismatch at position 8 and a single gap of length two inserted in the read after position 4

Fig. 2 .

 2 Fig. 2. Distribution of gap lengths in exome sequencing

Fig. 3 .

 3 Fig. 3. Processing times of needle, water, and GapMis-OMP for aligning 10, 000 pairs of sequences

Fig. 4 .

 4 Fig. 4. Processing times of GapMis-OMP for aligning 100, 000 query sequences and 100 target sequences

Data generated by the Exome Sequencing Programme at the NIHR Biomedical Research Centre at Guy's and St Thomas' NHS Foundation Trust in partnership with King's College London.

http://www.exelixis-lab.org/gapmis ⋆ This work was supported in part by the NSF-funded iPlant Collaborative (NSF grant #DBI-0735191)

up to 4× faster, the -t 8 version up to 7.8× faster, the -t 16 version up to 15.5× faster, the -t 32 version up to 30.6× faster, and the -t 48 version up to 41.7× faster than the -t 1 version.

The experiments were conducted on 1 node of a cluster architecture using 1 to 48 2.4 GHz AMD 6136 processors and 125 GB of main memory, and running the GNU/Linux operating system. GapMis-OMP is distributed under the GNU General Public License (GPL).


Martin Eklund

Ulf Norinder

Scott Boyer

Lars Carlsson

Application of conformal prediction in QSAR

QSAR modeling is a method for predicting properties, e.g. the solubility or toxicity, of chemical compounds using statistical learning techniques. QSAR is in widespread use within the pharmaceutical industry to prioritize compounds for experimental testing or to alert for potential toxicity. However, predictions from a QSAR model are difficult to assess if their prediction intervals are unknown. In this paper we introduce conformal prediction into the QSAR field to address this issue. We apply support vector machine regression in combination with two nonconformity measures to five datasets of different sizes to demonstrate the usefulness of conformal prediction in QSAR modeling. One of the nonconformity measures provides prediction intervals with almost the same width as the size of the QSAR models' prediction errors, showing that the prediction intervals obtained by conformal prediction are efficient and useful.

Introduction

Elucidating the structural properties of chemical compounds required to elicit a desired pharmacological effect (cure or alleviate a disease) and to concomitantly avoid toxicity is a fundamental issue in drug development. Quantitative Structure-Activity Relationships (QSAR) is a framework for employing statistical learning methods to predict the pharmacological effect and the toxicity (as well as other desirable properties) from chemical compounds' structures.

For making informed decisions based on predictions from a QSAR model, we are interested in the confidence in the predictions. Substantial efforts have been devoted to research on this topic within the QSAR community over the last decade and a number of methods have been suggested for estimating the confidence of QSAR predictions (see e.g. [START_REF] Netzeva | Current status of methods for defining the applicability domain of (quantitative) structure-activity relationships. The report and recommendations of ECVAM Workshop 52[END_REF][START_REF] Dragos | Predicting the predictability: a unified approach to the applicability domain problem of QSAR models[END_REF][START_REF] Jaworska | Towards optimization of chemical testing under REACH: a Bayesian network approach to Integrated Testing Strategies[END_REF] and references therein). These confidence estimates are typically based on the very loosely defined concept of a QSAR model's "applicability domain" (AD), which in [START_REF] Bassan | Computational Tools for Regulatory Needs[END_REF] was described as "the response and chemical structure space in which the model makes predictions with a given reliability". The assumption is that the further away a molecule is from a QSAR model's AD (according to some measure), the less reliable the prediction. The problem with the current approaches to estimating the prediction confidence martin.eklund@farmbio.uu.se in QSAR models is that their interpretation is cumbersome. For example, what does it mean that a chemical compound is outside the QSAR model's AD by a certain amount according to a certain measure? How does the plethora of different metrics used to define the AD relate to each other? What we ideally would like to know is in fact that a prediction is within a given prediction interval with a certain confidence (e.g. 80% or 95%).

In this paper we introduce conformal prediction [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF][START_REF] Shafer | A Tutorial on Conformal Prediction[END_REF] into the QSAR field. Conformal prediction uses previous data to determine prediction regions for new predictions. Given a confidence level (of say 80% or 95%), it produces a valid prediction region (i.e. a region that contains the true value with a probability equal to or higher than the given confidence level) under the assumptions that the observed data is i.i.d. [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. The conformal prediction framework provides a unified view of the different approaches to estimating a QSAR model's AD. Moreover, conformal prediction gives a natural and intuitive way of interpreting the AD estimates as prediction intervals with a given confidence.

The paper is organized as follows: We give a short introduction to QSAR modeling and AD estimation in Section 2. Continuing with Section 3 that contains a description of the methods (learner, nonconformity measure) used in the paper and how we apply them to QSAR modeling. The results are presented in Section 4 and a discussion of the results in Section 5.

QSAR

A molecule can be represented by an undirected labeled graph G = (V, E, A) (possibly embedded in IR 3 if the three dimensional structure of the molecule is taken into account), where the vertices V are the atoms and the edges E are the bonds. The atoms of a molecular graph are labeled by A, a set of atom types, which for instance can be the set of elements of the periodic table, the set of physicochemical properties of atoms, or any set of atom types provided by a molecular force field (see e.g. [START_REF] Halgren | Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94[END_REF]).

A descriptor function may be defined as a random variable

d : G → IR p . (1)
A realization x i of d applied to a molecular structure G i is called the descriptor vector (or simply descriptor ) of G i (NB: we will in this paper use "descriptor" and "attribute" interchangeably to refer to an element in x i ). Given a training set (x 1 , y 1 ), ..., (x n , y n ), where each y i , i = 1, ..., n, is the known activity of the molecule with descriptor vector x i , we can apply statistical learning methods to estimate a QSAR. A QSAR model fitted to (x 1 , y 1 ), ..., (x n , y n ) permits us to predict the activity of new (possibly not yet synthesized) chemical compounds.

Applicability domain

It is difficult for researchers to use the predictions from a QSAR model if there is no information available on whether the predictions are reliable or not. A large number of approaches have been suggested in the QSAR literature to assess whether a give prediction is reliable, typically by defining a QSAR's applicability domain. For example:

-Let μ and Ŝ be the respective estimates of the mean vector and the covariance matrix of the distribution P from which x 1 , ..., x n were sampled. A prediction of a new compound with descriptor vector x new is regarded less reliable the larger the distance D(x new ) = (x new -μ) T Ŝ-1 (x new -μ) and may be defined as "unreliable" if D(x new ) exceeds some predefined limit c. The hypersphere D(x) ≤ c is thus the AD.

-Estimate the density of P . This is often done by assuming independence between descriptors, or by projecting (x 1 , y 1 ), ..., (x n , y n ) to a lower dimensional space using e.g. principal component analysis (PCA). The density value at x new is determined and the lower this value is, the less reliable the prediction is. If it is lower than a predefined limit, the prediction may be deemed "unreliable". -A large number of methods are based on resampling (e.g. bootstrapping) from (x 1 , y 1 ), ..., (x n , y n ) to estimate the variability of the prediction of a new compound x new (see e.g. [START_REF] Bosnić | Comparison of approaches for estimating reliability of individual regression predictions[END_REF] for a survey of eight different resampling based methods). The AD is then implicitly defined as the region in descriptor space where the variability of predictions is lower than a given cutoff value.

A reader familiar with conformal prediction will realize that all these measures can be viewed as different nonconformity measures.

3 Methods

The signature descriptor

Let w be an atom of G, the signature of height h of w, σ h w , is a canonical representation of the subgraph of G containing all atoms that are at distance h bonds from w. It is usually encoded as a string, see Figure 1 and Table 1. The signatures can be viewed as letters of a finite alphabet, σ h w ∈ Σ and the occurrence of each letter is used as a descriptor value (i.e. each letter in Σ is represented by a given position in a descriptor vector and the value on that position is the number of occurrences of the letter in a molecule). In this study Σ was defined by the set of signatures defined at training time of each particular QSAR model. Typically the number of letters in Σ is much greater than the number of examples in a training set and the maximum number of letters present for a molecule will be the number of heights times the number of atoms. Detailed explanations on the algorithm of the signature descriptor and the usage of the signature descriptor can be found in [START_REF] Faulon | The signature molecular descriptor. 1. using extended valence sequences in QSAR and QSPR studies[END_REF][START_REF] Faulon | The signature molecular descriptor. 4. Canonizing molecules using extended valence sequences[END_REF].

The motivation for using the signature descriptor is that it has performed very well in benchmarking experiments on in-house AstraZeneca datasets. The limitation with the signature descriptor is that it only is a two-dimensional description of a molecule (as opposed to a potentially more information rich three

w h = 0 h = 1 h = 2 1 [C] [C]([C]=[C]) [C]([C](=[C])=[C]([C][N])) 2 [C] [C]([C]=[C]) [C]([C](=[C])=[C]([C][O])) 3 [C] [C]([C]=[C][N]) [C]([C](=[C])=[C]([C])[N]([C])) 4 [C] [C]([C]=[C][O]) [C]([C](=[C])=[C]([C])[O]) 5 [C] [C]([C]=[C]) [C]([C](=[C][O])=[C]([C])) 6 [C] [C]([C]=[C]) [C]([C](=[C][N])=[C]([C])) 7 [C] [C]([C]) [C]([C]([N]=[O])) 8 [C] [C]([C][N]=[O]) [C]([C][N]([C])=[O]) 9 [O] [O](=[C]) [O](=[C]([C][N])) 10 [N] [N]([C][C]) [N]([C]([C]=[O])[C]([C]=[C])) 11 [O] [O]([C]) [O]([C]([C]=[C]))
Table 1. The signatures of three different heights for paracetamol. w denotes the atom number and h the height.

dimensional description). However, valid three-dimensional representations are both costly to compute and there is also a great deal of uncertainty regarding what three-dimensional conformations a molecule would have for a particular problem, which introduces errors within a machine-learning based model.

Learning method

We used the support vector regression machine (SVM) [START_REF] Vapnik | Statistical learning theory[END_REF] implemented in the Java library LIBSVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] with a Gaussian radial basis kernel function

K(x, x ) = exp (-γ x -x 2 ). (2)
γ and the cost parameter C were either chosen in a cross-validated grid search or preset such that γ = 2 -9 and C = 100. During the grid search γ could take on vales 2 -s with s = 3, . . . , 10 and C values 10 t/2 with t = 0, . . . , 5.

Nonconformity measures

We will in this paper use inductive conformal prediction (ICP) due to its substantially smaller computational cost compared with transductive conformal prediction (TCP). We will therefore use nonconformity measures that are suitable for the inductive framework. These measures will essentially follow the nonconformity measures (31) and (32) defined in [START_REF] Papadopoulos | Regression conformal prediction with nearest neighbours[END_REF]. We will however make some changes compared to [START_REF] Papadopoulos | Regression conformal prediction with nearest neighbours[END_REF], which we point out below.

Let T be the training set (x 1 , y 1 ), ..., (x n , y n ) and let (x i1 , y i1 ), . . . , (x i k , y i k ) denote the k nearest neighbors of x i in T . For an example x i , we estimate the standard deviation of the labels of its k nearest neighbors as

s k i = 1 k k j=1 (y ij -ȳi 1,...,k ) 2 , (3)
where

ȳi 1,...,k = 1 k k j=1 y ij . (4)
The standard deviation of the labels is normalized with the median standard deviation of the k nearest neighbor labels of all training examples

ξ k i = s k i median({s k j : (x j , y j ) ∈ T }) . (5)
Let

d k i = k j=1 K(x i , x i ) + K(x ij , x ij ) -2K(x i , x ij ) = 2 k j=1 1 -K(x i , x ij ) (6)
be the sum of the distances in feature space between x i and its k nearest neighbors, where K(•, •) is the kernel function in Equation ( 2). This distance is different to the one used in [START_REF] Papadopoulos | Regression conformal prediction with nearest neighbours[END_REF] (where the Euclidean distance was used), however is natural to use in this work considering our choice of learning method. Similarly to (5), we divide [START_REF] Shafer | A Tutorial on Conformal Prediction[END_REF] with the median of the distances of all training examples from their k nearest neighbors:

λ k i = d k i median({d k j : (x j , y j ) ∈ T }) . ( 7
)
We now introduce the two nonconformity measures that we will use:

α i = y i -ŷi exp (γλ k i ) + exp (ρξ k i ) (8)
and

α i = y i -ŷi κ + ηλ k i + (1 -η)ξ k i , (9)
where ŷ is the SVM prediction. The measure ( 8) is identical to measure (32) in [START_REF] Papadopoulos | Regression conformal prediction with nearest neighbours[END_REF] (apart from the difference in distance function). Measure ( 9) is similar to measure (31) in [START_REF] Papadopoulos | Regression conformal prediction with nearest neighbours[END_REF], however we have introduced the parameter η to control the relative influence of λ k i and ξ k i . In this work, we used κ = η = 0.5 in all tests. The rationale behind the denominator in measures ( 8) and ( 9) is done to normalize with the expected accuracy of the underlying method (see [START_REF] Papadopoulos | Regression conformal prediction with nearest neighbours[END_REF] for a more detailed discussion). We ask for 95% prediction intervals in all tests.

Data description

We used the signature descriptor of height zero to three to compute descriptor vectors for the chemical compounds in the datasets described below (two of the datasets are publicly available, whereas three datasets are proprietary in-house AstraZeneca datasets). We then applied SVM and ICP as outlined in Sections 3.2 and 3.3 to derive QSAR models with prediction intervals. All datasets were split into training-, calibration-, and test sets. The training sets were used for model fitting and parameter optimization, the calibration set for computing the prediction region, and the test set for evaluating the model fits and the prediction regions. The calibration examples were chosen by random sampling (repeated between 1 and 10 times for different datasets). The testset was chosen by selecting the last (most recent) q examples from the dataset. The rationale behind this method for selecting the testset is that when designing new compounds there is usually distinct changes to the composition of the chemical compounds, and it can thus be argued whether new examples are i.i.d.

Datasets

The publicly available datasets contain the aqueous solubilities for a diverse set of 1297 organic compounds [START_REF] Huuskonen | Estimation of aqueous solubility for a diverse set of organic compounds based on molecular topology[END_REF] and the US Food and Drug Administration (FDA) recommended daily doses for 1215 pharmaceuticals. These datasets are available at http://cheminformatics.org/ and http://www.epa. gov/ncct/dsstox/sdf_fdamdd.html, respectively. We are not permitted to disclose the biopharmaceutical endpoints or the structures of the molecules in the AstraZeneca in-house datasets due to legal reasons, however summary statistics of these datasets are shown in Table 2 along with the same statistics for the public datasets.

Results

The main results are shown in Table 3 and Figure 2. Table 3 shows comparisons of the predictive region tightness for the different nonconformity measures outlined in Section 3.3 applied to the datasets described in Section 3.4, respectively. Figure 2 shows violin plots of the labels of the different datasets, the absolute value of the prediction errors and the widths of nonconformity measure (9) applied to them (a violin plot is a combination of a box plot and a rotated kernel density plot on each side of the box plot, see [START_REF] Hintze | Violin plots: A box plot-density trace synergism[END_REF][START_REF] Adler | vioplot: Violin plot[END_REF]).

Discussion

Conformal prediction formalizes the AD-based approaches for estimating the reliability of QSAR predictions. Different AD measures can simply be regarded as special cases of nonconformity measures. In this paper we used the nonconformity measures ( 8) and ( 9) with the 3 and 10 nearest neighbors, respectively. Measure (9) using the 10 nearest neighbors appears to produce the best results (giving the overall tightest predictive regions) on the tested datasets. The reason for this is that the variance of the labels among the near neighbors can be very high for a subset of the molecules in the datasets (in fact, this is a well-known phenomenon in QSAR modeling and has been termed "the QSAR paradox" or the "the Kubinyi paradox" in the literature [START_REF] Van Drie | Pharmacophore discovery-lessons learned[END_REF]). This leads to unnecessarily large prediction regions for measure [START_REF] Bosnić | Comparison of approaches for estimating reliability of individual regression predictions[END_REF] due to the exponentiation of the estimated variance among the near neighbors. Using 10 near neighbors instead of 3 alleviates this problem, since the effect on the variance of a problematic molecule is "diluted" to a larger extent. Further, using measure (9) avoids this problem and -for the large datasets where reasonably good predictions could be achieved (the AZ2 and AZ3 datasets) -generates informative prediction regions. The improvement of measure (9) over measure [START_REF] Bosnić | Comparison of approaches for estimating reliability of individual regression predictions[END_REF] suggest that it may be interesting to optimize the parameters in nonconformity measures during the training of QSAR models (analogously to how e.g. γ and C in the SVM are chosen using cross-validation). It may also be interesting to try different AD measures suggested in the QSAR literature and evaluate their efficiency (tightness) compared to measures [START_REF] Bosnić | Comparison of approaches for estimating reliability of individual regression predictions[END_REF] Table 3. The tightness and reliability results of nonconformity measure ( 8) and ( 9) on the test data. 3NN and 10NN denote the number of near neighbors used. Range is the range spanned by the labels in the different datasets (a constant was added to all labels for each dataset to make the range start at 0) and R 2 t is the R 2 statistic of the SVM predictions on the test datasets. Measure indicates which nonconformity measure that was used and number of runs shows the number of repetitions used when selecting the calibration set. Median, mean, and max width show statistics of the widths of the obtained prediction intervals. % outside shows the number of predictions outside the 95% prediction interval.

The number of predictions outside the 95% prediction interval is too high on two of the tested datasets (AZ2 and AZ3, see Table 3). This indicates that the i.i.d. assumption is not valid for these datasets; a quantile-quantile plot of the α i from the calibration set and the α i from the test set (plot not shown due to space constraints) and testing for equality of the distribution of α i from the calibration set and the α i from the test set further strengthen this conclusion (Kolmogorov-Smirnov test; p-value < 10 -8 ). QSAR models are often described as "global" or "local", where global models are derived on dataset containing structurally diverse compounds reflecting a range of different mechanistic actions, whereas local models are more focused and usually based on smaller sets of chemically similar compounds. In this paper we used ICP due to its computational efficiency, which is likely the only realistic approach for global QSAR modeling. However, for the smaller sets of compounds in local models, TCP makes more sense.

Fig. 1 .

 1 Fig. 1. Paracetamol is a drug molecule for treatment of pain and fever. The atom numbers are displayed in this figure.

Table 2 .

 2 Description of the datasets. The table shows the number of attributes for each dataset, as well as the number of examples (molecules) used for training, calibration, and testing. Whether γ and cost parameter C were preset or estimated using crossvalidation is also indicated in the "Preset" column, whereas the column "Random" indicates whether the calibration examples were chosen at random or if the q most recent ones were used.

	Dataset Attributes Training Calibration Test Preset Random
	Solu.	6862	926	99	257 No	Yes
	FDA	11295	901	99	215 No	Yes
	AZ1	18313	3331	99	391 Yes	No
	AZ2	52470	18320	499	871 Yes	No
	AZ3	82432	38763	499	2862 Yes	No

  and (9).

	Dataset Range	R 2 t	Measure Number Median Mean	Max	%
				of runs width width	width	outside
	Solu.	0 to 10.3 0.90 (8)/3NN 10	14	9.8 × 10 15 7.5 × 10 18 1.8
		0 to 10.3 0.90 (8)/10NN 10	140	3.8 × 10 12 2.2 × 10 15 2.1
		0 to 10.3 0.90 (9)/3NN 10	4.4	7.2	54	1.8
		0 to 10.3 0.90 (9)/10NN 10	2.0	2.4	12	3.9
	FDA	0 to 6.7 0.45 (8)/3NN 10	2.8	1.6 × 10 12 8.5 × 10 14 1.1
		0 to 6.7 0.46 (8)/10NN 10	2.2	4.1 × 10 8 2.2 × 10 11 4.0
		0 to 6.7 0.44 (9)/3NN 10	3.1	4.8	81	1.2
		0 to 6.7 0.45 (9)/10NN 10	2.0	2.5	22	2.4
	AZ1	0 to 3.5 -0.10 (8)/3NN 1	2.0	3.8	369	0.8
		0 to 3.5 -0.07 (9)/10NN 5	1.0	1.1	3.4	6.5
	AZ2	0 to 5.4 0.75 (8)/3NN 1	1.0	4.5 × 10 8 2.0 × 10 11 4.0
		0 to 5.4 0.70 (9)/10NN 5	0.7	0.7	2.7	17
	AZ3	0 to 5.4 0.74 (8)/3NN 1	0.7	4.8 × 10 2 1.2 × 10 6 5.7
		0 to 5.4 0.70 (9)/10NN 1	0.6	0.7	3.1	14

In conclusion, conformal prediction provides a unified view of previous attempts to estimate the reliability of QSAR prediction. By regarding AD estimates as nonconformity measures, conformal prediction permits an intuitive interpretation (p-values and prediction intervals) of AD measures, which are otherwise often very difficult to interpret.


Christos A Frantzidis
email: christos.frantzidis@gmail.com

Maria D Diamantoudi

Eirini Grigoriadou

Anastasia Semertzidou

Antonis Billis

Evdokimos Konstantinidis

Manousos A Klados

Ana B Vivas

Charalampos Bratsas

Magda Tsolaki

Constantinos Pappas

Panagiotis D Bamidis

A Mahalanobis Distance Based Approach towards the Reliable Detection of Geriatric Depression Symptoms Co-existing with Cognitive Decline

Keywords: Geriatric Depression, Electroencephalography, Mild Cognitive Impairment, Mahalanobis Distance Classifier, Wavelet Entropy, Neuropsychological Estimation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The overall prevalence rate of depressive disorders among the elderly varies between 10% and 20% [START_REF] Barua | Prevalence of depressive disorders in the elderly[END_REF]. According to projections by the World Health Organization (WHO) depression will be the leading cause of Disability Adjusted Life Years (DALYs) lost in high income countries in 2030 [START_REF] Kok | Can reporting heterogeneity explain differences in depressive symptoms across Europe?[END_REF]. Contrary to popular belief, depression is not a natural part of aging and often reversible with prompt and appropriate treatment [START_REF] Kurlowicz | The Geriatric Depression Scale (GDS). Best Practices in Nursing[END_REF]. However, if left untreated, it leads not only to considerable influences on functional impairment and disability [START_REF] Veer-Tazelaar | Prevention of anxiety and depression in the age group of 75 years and over: a randomized controlled trial testing the feasibility and effectiveness of a generic stepped care program among elderly community residents at high risk of developing anxiety and depression versus usual care[END_REF], [START_REF] Tam | Cognitive and functional impairment in chinese elderly with late-onset depression[END_REF], but also to increased cardiovascular diseases, mortality risks and a decreased quality of life, along with increased health care utilization and medical costs. Hence, early detection and treatment is important in chronically ill elderly persons, thereby preventing a downward spiral [START_REF] Lamers | Effectiveness and cost-effectiveness of a minimal psychological intervention to reduce non-severe depression in chronically ill elderly patients: the design of a randomised controlled trial[END_REF], [START_REF] Wint | Depression: a shared risk factor for cardiovascular and Alzheimer disease[END_REF].

Depression is highly prevalent among individuals with dementia of any degree and any type [START_REF] Barca | Factors associated with a depressive disorder in Alzheimer's disease are different from those found from other dementia disorders[END_REF]. Depressive symptoms are also common in mild cognitive impairment (MCI) and appear to hasten conversion from MCI to clinical Alzheimer 's disease (AD) [START_REF] Rapp | Cognitive decline in patients with dementia as a function of depression[END_REF], [START_REF] Royall | Depressive symptoms predict longitudinal change in executive control but not memory[END_REF]. Depression, subcortical dementia, and normal aging may all have similar neurobehavioral manifestations, a fact that makes the differentiation between depression and dementia in the elderly extremely complicated [START_REF] Barua | Prevalence of depressive disorders in the elderly[END_REF]. Still controversial is whether depression prior to dementia represents a risk factor for dementia or a prodromal feature of dementia and the most likely scenario seems that depression can be both an early or midlife risk factor for dementia and also an early sign of incipient dementia. For instance, in a late stage of dementia depression could be due to specific structural or biochemical changes in the brain [START_REF] Barca | Factors associated with a depressive disorder in Alzheimer's disease are different from those found from other dementia disorders[END_REF]. Furthermore, maturational changes in the brain may increase vulnerability for the development of late-life depression [START_REF] Wilkins | Latelife depression with comorbid cognitive impairment and disability: non pharmacological interventions[END_REF] while in late-onset depression, structural imaging studies point not only to changes in deep frontal white matter but also to structural changes in the hippocampus and prefrontal cortex [START_REF] Rapp | Cognitive decline in patients with dementia as a function of depression[END_REF]. These alterations including connections with the limbic system have been hypothesized to increase vulnerability for depression and contribute to the manifestation of declines in cognition and subsequent executive dysfunction associated with late-life depression [START_REF] Wilkins | Latelife depression with comorbid cognitive impairment and disability: non pharmacological interventions[END_REF]. The widespread structural brain changes mentioned include structures which are affected by AD pathology too [START_REF] Royall | Depressive symptoms predict longitudinal change in executive control but not memory[END_REF] rendering the diagnosis of depression in the context of dementia complicated [START_REF] Fiske | Depression in Older Adults[END_REF].

A closer look at the relationship between depression and dementia development illustrates complexity of this relationship. For instance, cognitively intact elderly individuals who develop depression are at increased risk of subsequent MCI [START_REF] Geda | Depression, Apolipoprotein E genotype, and the incidence of mild cognitive impairment: a prospective cohort study[END_REF]. In fact in cognitively intact subjects, symptoms of depressive mood were related to AD development [START_REF] Palmer | Predictors of progression from mild cognitive impairment to Alzheimer disease[END_REF]. In a different study it is demonstrated that although not all patients with MCI develop dementia, in those who do, those who present with depression as well, are at more than twice the risk of developing AD (as those without depression) [START_REF] Modrego | Depression in patients with mild cognitive impairment increases the risk of dementia of Alzheimer type[END_REF]. Though it is not clear yet whether depression is a cause or a consequence of cognitive impairment in older populations, it has been found nevertheless that there is indeed an association between MCI in specific and depression, which in turn leads to a worse clinical outcome of depression (and potentially faster cognitive decline) [START_REF] Defrancesco | Association of Mild Cognitive Impairment (MCI) and depression[END_REF]. Depression, however, has also been found to be a risk factor for AD in individuals without cognitive impairment [START_REF] Wilson | Depressive symptoms, cognitive decline, and Risk of AD in older persons[END_REF]. In a review of relevant studies, it was concluded that depressive symptoms are not related to progression from MCI to dementia [START_REF] Palmer | Predictors of progression from mild cognitive impairment to Alzheimer disease[END_REF], while it has not been precluded that depression may contribute to cognitive deficits in the elderly. If depressive symptoms in certain MCI cases are indeed related to a different etiology, separate from a neurodegenerative mechanism [START_REF] Palmer | Predictors of progression from mild cognitive impairment to Alzheimer disease[END_REF] that remains to be examined.

The role of depression in dementia is without doubt complicated. Depression in older individuals may mimic or exacerbate symptoms of other progressive dementing conditions, thus making differential diagnosis difficult [START_REF] Derix | Neuropsychological abnormalities in depression: Relation between brain and behaviour[END_REF]. Despite the large literature and findings on the connections between depression and dementia (whether that be preclinical, MCI or AD), there is no integrated method to strengthen the already existing approaches and lead to more solid differential diagnoses. Bearing that in mind, the current study sought to set the basis for a stronger diagnosis method. Therefore, it is clear that improved diagnostic tools are essential to further clarify the relationship between depressive symptomatology and dementia, and mainly MCI, in older populations. Depressive symptomatology needs to be examined through clinical diagnosis, supported by neuropsychological assessment, adding to data from neurophysiological measures. Neuropsychological assessment of an individual with clinical depression is expected to underline impairments typically seen in cognitive faculties such as attention, memory components such as short-term recognition and recall or long-term episodic and semantic memory, abstraction and executive functions, language, motor skills as well as visuospatial functions [START_REF] Derix | Neuropsychological abnormalities in depression: Relation between brain and behaviour[END_REF], or show slowed mental processing [START_REF] Lezak | Neuropsychological Assessment[END_REF]. Neuropsychological assessment will allow for the brain-behavior relationship to be explored and to further develop inferences on the neuroanatomical substrates of the disease [START_REF] Derix | Neuropsychological abnormalities in depression: Relation between brain and behaviour[END_REF].

Several research questions arise regarding geriatric depression and pathological aging due to neurodegeneration. Can MCI develop in all patients with depressive symptomatology? And if yes, how soon, and to what extend? What changes in the brain due to aging may facilitate this process of MCI development, with or without depressive symptoms present? Can MCI along with depression be viewed as a stronger predictor of AD development than MCI or depression alone? Which neuropsychological measures in particular may be useful in the differential diagnosis of depression from early dementia types such as MCI? If they are used parallel to neurophysiological measures such as EEG, will their predictive value increase, and will this combination of measures assist in differentiation and further treatment?

Our hypothesis is that a large portion of MCI patients suffer simultaneously from geriatric depression. The coexistence of depression reduces their commitment to interventions regarding dementia and prompts them to self-occlusion and further deterioration. Despite depression diagnosis is of vital importance, there is not a plethora of objective and reliable diagnostic tools. So, it is mainly supported on subjective questionnaires and to the therapist's experience. Aiming to provide a robust and simple computerized tool that may be used to clinical practice , we propose the use of resting state EEG data and a pattern recognition technique towards the detection of the depressive symptoms.

Methodology

The study employed 66 (six male) elderly participants recruited from day care centers and the Greek Association of Alzheimer's Disease and Related Disorders. They were subdivided equally in two groups of 33 participants / group. The first group consisted of elderly participants diagnosed with geriatric depression and memory loss or minor cognitive deficits like Mild Cognitive Impairment. The second group consisted of elderly participants free of depressive symptoms. Both groups were matched regarding their age and cognitive status. Inclusion criteria were age of 58 and older, approval to participate from their personal doctor and signing the informed consent. Exclusion criteria were dementia diagnosis, drug or substance abuse, recent usage of cholinesterase inhibitors or anti-depressive drugs (less than 3 months), severe neurological disorder due to ischemic attack and severe mobility problems. This study was part of the Long Lasting Memories EU FP7 funded program which aimed to promote healthy and independent aging through cognitive and/or physical training. Prior to their participation they were informed about the research procedure and aims. Then, they had the chance to ask anything they wanted regarding the study and signed an informed consent form. The study was approved by the local ethical committee. Group details are given in Table 1. The neuropsychological examination focused on the presence of cognitive impairment relevant to MCI and to depression. The Short Form Geriatric Depression Scale (GDS), which consists of 15 questions, was used for the depression identification. The GDS is a tool with validity and reliability supported through both clinical practice and research. The participant's cognitive status was assessed through the Mini Mental State Examination (MMSE) and the Montreal Cognitive Assessment (MoCA) test. The MMSE is commonly used for memory complaints of memory problem, since it offers a general screening tool for a number of different mental abilities including memory, attention and language. It is often combined with the MoCA examination which is a rapid screening instrument for mild cognitive dysfunction assessing different cognitive domains. Verbal memory was assessed through the California Verbal Learning Task (CVLT). Conclusions regarding psychomotor speed and executive control function were derived by Trail-Making Test Part A and B and working memory was assessed with Digit Span Task. The ability to perform daily self-care activities was measured by the Instrumental Activities of Daily Living (IADL) test. Finally, the WHOQOL-BREF, a shorter version of the World Health Organization Quality Of Life Instrument, was used to estimate the subjective quality of life including physical, psychological, social and environmental domains.

Neurophysiological data acquisition was employed in order to provide a direct window of mature brain functioning enabling thus the reliable detection and quantification of cognitive decline and brain deviations due to depression symptoms. Participants were sat in a comfortable armed chair during the neurophysiological examination, performed through a Nihon Kohden JE-207A device equipped with 57 active electrodes attached on a cap fitted to the participant's scalp. Electrode impedances were kept lower than 2 ΚΩ. Reference electrodes were positioned at the mastoids. The sampling rate was set at the 500 Hz. The data acquisition procedure is visualized in Figure 1. The experimental procedure involved resting state EEG data acquisition with eyes closed. The data acquisition lasted for approximately four minutes and was part of a longer neurophysiological examination within the context of the Long Lasting Memories screening procedure [START_REF]The Long Lasting Memories (LLM) project[END_REF].

The pre-processing was performed through the MATLAB programming environment and the EEGLAB graphic user interface. The EEG signals serve as an input to the Independent Component Analysis (ICA) algorithm. The ICA algorithm computes the source estimations that are hypothesized to produce the input EEG signals. Thanks to this kind of analysis, noise sources (e.g. eye blinks, muscle artifacts, heart modulation) are recognized and rejected. After the removal of the noisy source components, the whole EEG data were visually inspected by two experts in order to identify data segments that are temporally contaminated with noise. These short time segments are manually rejected. Then, a continuous interval of 20 seconds of artifact-free EEG data was selected for the feature extraction phase.

Figure 1 Visualization of the data acquisition equipment

The 20 seconds EEG data were divided in 150 windows of 128 milliseconds duration. The Discrete Wavelet transform (DWT) at periodization mode was then employed in order to compute the wavelet coefficients for each time window. The duration of the time window was set at 128 ms, which was the minimum time duration that needed for the computation of one wavelet coefficient correspond to the delta frequency band. The DWT scheme performs iterative time-frequency decomposition (first layer of Figure 2) in order to extract the activity of each frequency band with optimal resolution through recursive filtering (second layer of Figure 2) [START_REF] Frantzidis | Toward Emotion Aware Computing: An Integrated Approach Using Multichannel Neurophysiological Recordings and Affective Visual Stimuli[END_REF]. The wavelet coefficients that correspond to the delta (0.5-4 Hz), theta (4-8 Hz), alpha (8-16 Hz), beta (16-32 Hz) and gamma (32-62 Hz) band are used in order to compute the energy (Ej, j=1…5) of each one of the aforementioned frequency bands, through the following type:

(1). Then, the relative energy of each frequency band is computed (band's contribution to the total EEG energy). Finally, the synchronization degree between the energy distributions of each electrode pair is extracted (third layer of Figure 2) according to the following type:

(2), where p and q the energy distributions of the two electrodes forming the electrode pair under consideration.

The electrodes were divided according to their topological distributions in 12 categories: prefrontal, frontal, central, frontocentral, frontotemporal, temporal, anteriofrontal, parietal, parietotemporal, occipital, parietooccipital and centroparietal. So, except from the synchronization features, mean energy values for each one topological cluster were also computed and served as additional features. Therefore, the pool of features consists of both synchronization and energy characteristics.

The huge amount of features extracted from the relative wavelet entropy analysis (57 × 57 = 3249 features) pose the issue of feature selection. Therefore, a strict threshold was defined. The threshold absolute value was set at 0.06, allowing both synchronization and desynchronization features to be selected. The strict threshold resulted in the selection of three synchronization features and seven desynchronization features. The Mahalanobis Distance (MD) classifier was employed [START_REF] Frantzidis | Toward Emotion Aware Computing: An Integrated Approach Using Multichannel Neurophysiological Recordings and Affective Visual Stimuli[END_REF]. It is a simple yet robust distance based classifier, which is defined by the following formula:

(3)

Let assume that we have N features that would be employed for the classification scheme. The average value (centroid) of each feature forms the N-dimensional μ vector. S denotes the covariance matrix and the x variable denotes each instance (one value from each feature corresponding to a single participant) as defined in ( 4):

(4)

Results

The huge amount of the features extracted from both the rhythmic activity computation and the synchronization evaluation were submitted to statistical analysis through Student t-test. The statistical significance threshold was set at p <0.05. Bonferroni correction was also used. This procedure served as a pre-processing feature evaluation in order to detect the most important features that could differentiate the two groups. Regarding the synchronization features, the most prominent differences were found among the following electrode pairs: F8-Afz, Fz-FC3, P2-F7, Af4-C2, CP1-Af4, FC2-Af4, C2-Af4, CPz-FC3 and F7-P2.

Rhythmic activity energy features were also statistically significant. High frequency oscillations occurring within the gamma band differentiated the two groups on anteriofrontal and frontocentral locations. Beta band features were significant on central and parietal positions. No feature was highly statistical significant within the alpha frequency band. Theta energy features from parietooccipital and prefrontal areas were also used. Features extracted from the delta frequency band proved to best differentiate the two groups on several areas like centroparietal, parietooccipital, occipital, central, frontocentral and temporal locations. These features (10 from the synchronization analysis and 11 from rhythmic activity computation) were used in order to discriminate neurophysiological alterations between depressed and nondepressed patients that also suffer from subjective or mild cognitive impairment.

The classification rates for the two groups (depressive & controls) and the overall accuracy rate is depicted in Figure 3. The number of features used was 21. The classification accuracy for the patients suffering from depression was 93.94% (31/33) and 84.85% (28/33) for the age and cognitive status matched control group. So, the overall accuracy rate reached the 89.39%. According to these results the number of true positives is 31/33, false positives are 5/33, true negatives are 28/33 and false negatives are 2/33.

Discussion

The proposed Mahalanobis Distance (MD) based tool is one of the early attempts aiming to detect depression symptoms that co-exist with cognitive deficits characterizing pathological aging. Therefore, the problem's difficulty is much more complex since brain alterations occur both due to depression and neurodegeneration. Its clinical utility would be of major importance, since it would be used as a screening tool prior to non-pharmacological interventions aiming to delay or prevent dementia onset. Therapy dropout ratings are mainly attributed to depressive patients due to less commitment and reduced socialization. Therefore, such a system would recognize depressive participants at high risk of not finalizing the program (therapy wise) and special care would be given towards their encouragement to complete it and towards providing them with positive feedback. Towards this direction, it is important to introduce a multi-modal screening tool that processes neurophysiological and neuropsychological data, self reports and clinical diagnosis in order to provide a definitive depression diagnosis through a data mining tool.

Despite the promising results, there are still study limitations such as the lack of gender balanced participants. The sampled population of this study did not have a clear cut diagnosis of clinical depression (only depressive symptomatology), making further comparisons of groups more difficult. For instance, the results of the short form GDS are not necessarily suggestive of depressive symptomatology, an issue of concern in the current study. Moreover, type of depressive symptomatology was not examined, such as anxiety versus depressed mood types of symptoms, which appear to present differently in the presence of cognitive difficulties. Moreover, the sample size is relatively small in order to extract normative ratings that could establish a standardized diagnostic tool.

In future research, clinical depression should be diagnosed by means of a full clinical interview, fortified with the inclusion of more detailed self report measures. Furthermore, more detailed neuropsychological assessment should be performed to assist in the differentiation between dementia and depression, and examine further the etiology. Repeat testing will provide valuable information on depression and dementia features alike. To summarize, this study proposes a novel recognition scheme aiming to form a screening tool in the under-studied field of geriatric depression occurring simultaneously with cognitive neurodegeneration. Potential applications except from the definitive diagnosis of depression would be the personalized selection of cognitive and physical dementia-related intervention treatment and selection of an appropriate way of providing emotional feedback to elderly in order to increase their self esteem, socialization and independent living.

Figure 2

 2 Figure 2 Visualization of the iterative five-step frequency decomposition scheme adopted by the Discrete Wavelet Transform (DWT) (layer 1); the recursive filtering framework used for the extraction of the activity of each frequency band under consideration (delta, theta, alpha, beta and gamma) (layer 2); and the Relative Wavelet Entropy (RWE) synchronization matrix (layer 3).

Figure 3

 3 Figure 3 Classification accuracy of the Mahalanobis Distance (MD) based classifier.

Table 1 Group details regarding mean age, depression and cognitive status

 1

	Group	Age	GDS	MMSE	MoCA
	Depression	69.818	7 ± 1.62	25.455 ± 2.762	22.485 ± 3.589
	Control	70.333	1.667 ± 1.574	25.727 ± 1.892	22.273 ± 2.798

Acknowledgments. This work was granted by the USEFIL ICT FP7 EU funded program.


Khuong Nguyen
email: khuong@cantab.net

Zhiyuan Luo
email: zhiyuan@cs.rhul.ac.uk

Conformal Prediction for Indoor Localisation with Fingerprinting Method

Keywords: Conformal prediction, indoor localisation, Fingerprinting

Indoor localisation is the state-of-the-art to identify and observe a moving human or object inside a building. Location Fingerprinting is a cost-effective software-based solution utilising the built-in wireless signal of the building to estimate the most probable position of a real-time signal data. In this paper, we apply the Conformal Prediction (CP) algorithm to further enhance the Fingerprinting method. We design a new nonconformity measure with the Weighted K-nearest neighbours (W-KNN) as the underlying algorithm. Empirical results show good performance of the CP algorithm.

Introduction

The purpose of indoor localisation is to identify and observe a user inside a building. Global Positioning System (GPS) has long been an optimal solution for outdoor localisation, yet the indoor counterpart remains an open research problem, because the complex building infrastructure hinders the GPS signal. In this paper, we applied the Conformal Prediction (CP) algorithm to enhance an effective indoor tracking solution known as "Fingerprinting method" [START_REF] Bahl | RADAR: an in-building RF-based user location and tracking system[END_REF]. We designed a new nonconformity measure with the W-KNN as the underlying algorithm. To the best of our knowledge, we are the first to apply CP for the Localisation purpose in general, and for the Indoor Localisation context in particular. The logical progression of the paper is graphically depicted in Figure 1.

Indoor Localisation Problem

A user can be coarsely identified at room-level or precisely localised at sub-room level. The finest-grained systems offer up to 3cm accuracy, 95% of the time, using intensive hardware implementation [START_REF] Want | The Active Badge location system. ACM Transaction on Information Systems[END_REF][START_REF] Ward | A New Location Technique for the Active Office[END_REF]. However, they demand massive investment to install and maintain. Balancing between fine-grained tracking and affordable deployment is desirable for recent approaches, which concentrates purely on software solutions. This section outlines the recent approaches with the indoor localisation, as a background to apply Conformal Prediction (CP) in the next section.

Coarse-grained Tracking

The most basic indoor localisation comes from the proximity-based tracking idea. Whenever a user and a base station can communicate, the user's location is determined as the location of the station. This method is further enhanced by dividing the tracking space into grids. The stations are strategically placed in such a way that each grid block is overlapped by the signal from as many different stations as possible (Figure 2). Thus, instead of coarsely predicting the user's location to be somewhere within the station's broadcasting range, the accuracy is improved by interpreting the user's location to be the overlapped portion of the stations the user sees. However, many stations must be deployed to have good signal coverage. Overall, coarse-grained tracking idea is great for localising user at the roomlevel resolution. The next section discusses how to identify the user at much finer-grained sub-room resolution.

Location Fingerprinting

Coarse-grained tracking is improved further by analysing the wireless signal strength, based on the fact that the radio signal attenuates and gets weaker as it travels in the air. There are two popular measurements to roughly represent the distance between a user and a station: the Received Signal Strength Indication (RSSI) and the Link Quality (LQ). However, two distinct locations might not have a linear relationship in terms of RSSI/LQ and the distance between them. This phenomenon is caused by the signal being blocked by the indoor objects, known as the multipath issue. An elegant solution -the "Location Fingerprinting" method has been widely adopted for its simplicity and efficiency [START_REF] Bahl | RADAR: an in-building RF-based user location and tracking system[END_REF]. This method utilises the built-in wireless signal of the building to survey a signal-to-position mapping database beforehand, which is known as the off-line stage. By surveying the whole signal variation at each position, Location Fingerprinting does not calculate the distance from the user to the station based on the signal strength. Instead, the system applies pattern-matching algorithms of the real-time signal and the database record to estimate the most probable position during the on-line stage.

In comparison to proximity-based tracking, this solution offers much higher fine-grained tracking even with a few stations. The overall accuracy of the system relies on two factors. First, if the database is poorly constructed during the offline stage, the correct estimation will degrade. Second, the prediction algorithm chosen in the on-line stage will determine how good the estimated position is. The Fingerprinting method is depicted in Figure 3.

Modelling the Indoor Localisation Problem

Without any loss of generality, let us formally model the Location Fingerprinting method. A general example scheme can be found in Figure 4, where there are N = 6 base stations (A-F) deployed in a tracking zone. A single location L is modelled in a 3-dimensional space as L = (d x L , d y L , d z L ). The signal strength RSSI between the user at this location L and all N stations is modelled as an N-tuple RSSI L = (s 1 , s 2 , . . . , s N ), where s i is the signal strength received from the station i (1 ≤ i ≤ N ). A database B organises the surveyed locations in terms of RSSI and the corresponding physical co-ordinate as shown in Table 1. The task is: given an RSSI tuple RSSI unknown = (s 1 , s 2 , . . . , s N ) of an unknown location inside the tracking zone, the system compares this tuple with all recorded fingerprints in the database B to estimate a closest match, which will tell the correct co-ordinate (d x , d y , d z ) of the unknown location. The main challenge is that RSSI measurement at the same location can vary from time to time, and distinct locations might have a similar RSSI combination due to human movements, humidity, furniture re-arrangement, as well as the multipath fading of the indoor environment.

Conformal Prediction for Indoor Localisation

Having discussed the problem of Indoor Localisation, this section explains and applies the Conformal Prediction (CP) algorithm to address the specific indoor problems outlined previously. Particularly, we show why CP is a suitable algorithm for our purpose. The next section will compare its performance with other traditional algorithms.

Conformal Prediction

Conformal Prediction is a recently developed machine learning algorithm, which uses experiences in the past to confidently and precisely predict the outcome of a new sample [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF][START_REF] Shafer | A Tutorial on Conformal Prediction[END_REF]. However, what differentiates CP from other similar machine learning algorithms is the ability to produce prediction region for the given confidence level parameter. It has been mathematically proved that the prediction region generated by CP is valid in online setting [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. In other words, for a confidence level at 95%, the correct estimated position is expected to be included in the prediction region at least 95% of the time. Achieving the prediction result in such a powerful sense, however, CP demands a relatively weak assumption that the training database and the new sample to be classified are generated from the same distribution independently. CP has been successfully used in many applications such as medical diagnosis and network traffic classification [START_REF] Bellotti | Qualified Predictions for Microarray and Proteomics Pattern Diagnostics with Confidence Machines[END_REF][START_REF] Dashevskiy | Reliable Probabilistic Classification of Internet Traffic[END_REF].

In general, CP is determined by some nonconformity measure [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. The nonconformity measure is a real-valued function A(B, z) measuring how different a sample z is to the training database B. Whenever a new sample needs to be classified, we exhaustedly test every possible label recorded in the training data, and CP tells us whether to accept that label or not by a test of randomness for a given significance level. CP uses existing machine learning algorithms such as nearest neighbours, SVM and neural network as the underling algorithms to compute a "nonconformity measure" and generate the prediction region [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. Regardless of the chosen nonconformity measure, the set of locations predicted by CP is always valid in the online setting. However, the efficiency, in other words, the tightness of the prediction region is affected. Ideally, it is more preferred to have as few predictions in the prediction region as possible, without sacrificing the confidence level too much. Each problem requires a customised nonconformity measure to fit the purpose, thus opening more research opportunities. In the next part, we design a new nonconformity measure for the indoor localisation problem.

Applying Conformal Prediction into Fingerprinting method

Giving a training database B mapping physical location co-ordinate to correct Bluetooth signal pattern, and a signal fingerprint at an unknown location, CP precisely predicts a set of locations in the database, which likely matches this new signal fingerprint. This can be formulated as a classification problem, because we divide the physical location into grid points, and the label set is finite. The signal strengths are regarded as the object set X, and the physical locations are regarded as the label set Y. We will apply CP using both the old examplesthe training database B = (z 1 , z 2 , . . . , z n-1 ), and the signal fingerprint of the unknown location (as a new object for z n ). Each sample z i is a combination of the signal strength RSSI i = (s i 1 , s i 2 , . . . , s i N ) and the co-ordinate

L i = (d x i , d y i , d z i ). A prediction region of K examples is R ε (L 1 , L 2 , . . . , L K ) ⊂ Y.
First, CP demands that the training database and the new real-time sample must be generated from the same distribution. This assumption holds because all signals are generated by the same set of base stations independently. Second, to evaluate the difference amongst the samples, the 'Weighted K-nearest neighbours' (W-KNN) [START_REF] Duda | Pattern Classification[END_REF] is chosen as the underlying algorithm to compute the 'nonconformity score' α. The reason we opted for the 'Weighted version' of the traditional K-nearest neighbours algorithm is because of the multipath problem of the indoor signals, which cause many locations further apart to have a similar combination of signal strength. If we treat each of the K neighbours equally, the final estimated position is shifted towards the wrong neighbours further away. By considering the "weight", corresponding to the inverse "signal distance" between each neighbour and the unknown position, the closer neighbours are prioritised over further away ones. The outline of the Conformal Prediction algorithm with W-KNN is presented in Algorithm 1.

Given a training database B, and a new example with a known signal strength pattern but unknown co-ordinate, CP assumes the correct position to be each of the recorded locations in Y. First, we use the Euclidean approach to calculate the "distance" between two signal strength sequences RSSI P and RSSI Q .

dist(RSSI P , RSSI Q ) = (s P 1 -s Q 1 ) 2 + . . . + (s P N -s Q N ) 2 . ( 1
)
Using Equation ( 1), we can find K entries (z 1 , . . . , z K ) in the database B which have the smallest distances dist(RSSI i , RSSI u ) to the unknown location U , as well as having the same co-ordinate label

L i = (d x i , d y i , d z i ) with U (1 ≤ i ≤ K).
Then, a weighted average location is calculated from these K entries, giving an estimated position e same = (d x s , d y s , d z s ). The below equation is repeated for each co-ordinate of e same , where δ is a very small constant to prevent division by zero.

d x s = K i=1 1 dist(RSSI i , RSSI u ) + δ d x i K i=1 1 dist(RSSI i , RSSI u ) + δ . ( 2
)
Similarly, using Equation ( 1), we find another K entries (z 1 , . . . , z K ) in the database B which have smallest distances dist(RSSI i , RSSI u ) to the unknown location U , but having different co-ordinate label

L i = (d x i , d y i , d z i ) to U (1 ≤ i ≤ K).
-d x d ) 2 + (d y s -d y d ) 2 + (d z s -d z d ) 2 . ( 3
)
Using Equation (3), we calculate the nonconformity score α i for every entry in the database B. Then, we count the number of α i which is larger than the nonconformity score α u of the new sample location U , and divide the total number of entries in the database B to have the p-value number [START_REF] Vovk | Algorithmic Learning in a Random World[END_REF]. Given a significance level ε beforehand (such as ε = 0.05), the current assumed coordinate label is accepted as the label for the new sample, if and only if the p-value ≥ ε.

Performance Evaluation

Having shown how to apply CP into the Indoor Localisation context, we perform experiments to evaluate the performance of CP with K-NN, and that of the W-KNN alone.

Conformal Prediction with W-KNN

For off-line evaluation of CP with W-KNN, we use the Bluetooth signal location database described in [START_REF] Nguyen | Robot-based Evaluation of Bluetooth Fingerprinting[END_REF]. There are 1,052 examples in the training database, and 260 examples in the testing set. Figure 5(a) shows that CP using W-KNN as the underlying algorithm outperforms the W-KNN algorithm alone by 8%-13%. Since CP returns a set of locations, while traditional W-KNN returns just one estimated location for each test sample, we averaged the prediction region returned by CP as a single location for comparison. We compared the performances of CP with the W-KNN algorithm using a Cumulative Distribution Function (CDF) plot (Figure 5(a)). At a significance level ε = 0.13 and K = 16, the CP with W-KNN has fewer than 1.5m location error in more than 90% of the time; or it has fewer than 1m location error, 40% of the time. Comparing to the un-Weighted version of K-NN, CP with W-KNN outperforms CP with K-NN by 25% (Figure 5(b)). This enhancement is similar to what reported for W-KNN and K-NN without CP [START_REF] Nguyen | Robot-based Evaluation of Bluetooth Fingerprinting[END_REF]. Varying the ε and K does affect the accuracy, as explained separately in the next parts.

To verify the validity of the proposed nonconformity measure, we performed leave-out-one on-line cross validation with the above training data only. For L = 1, 052 samples, we used L-1 samples for training and the remaining one for testing. At the same ε = 0.13, the error rate was achieved at 10.8%, which remains valid.

Confidence Level vs. Prediction Region

To study the credibility of the prediction region, we evaluate different confidence levels (see Table 2). The CP Error Rate is the percentage in which CP does not produce a prediction region containing the exact location. For simplicity, we considered a 2-dimensional database, where the prediction region is presented as a set of 2D locations concentrated around a peak position. We observed that the size of the prediction region decreases as the significance level ε increases. When the ε went past 0.38, the prediction region was empty for all testing points. We had precisely 1 prediction label when ε was around 0.29 for 70% of the test samples. By dividing the training database into 2 equal sets of samples, and used the second set as a test set and the first set as a training set, we worked out that a reasonable confidence level was around 87% with a prediction range of 1.5m -2m, which is acceptable for tracking users at sub-room resolution. There was only 3.5% of samples with more than 3.5m error. Increasing the confidence level to above 90% allows too many samples getting in the prediction region.

Calibrating the K Parameter

Finding an optimal K parameter for the W-KNN is challenging. Each environment and each signal property require the K parameter to be calibrated exclusively. In our system, we experimentally found K = 16 to be an optimal value across all testing points. However, some locations had smaller location error with different K value. Generally, starting from K = 1, which is equivalent to 1-nearest neighbour only, the accuracy tends to increase as K increases, up to a certain threshold (K = 16 in our case). Then, it begins to decrease as K continues to increase. In our experience, from K = 16 onwards, the performance is only 7%-10% affected by increasing K.

Conclusion and Further Work

This section outlines the main contributions of this paper, and the scheme to apply CP for tracking a mobile user.

Contributions

We have demonstrated the application of Conformal Prediction into the Indoor Localisation context. To the best of our knowledge, we are the first to implement Conformal Prediction algorithm into the Localisation problem in general, and the Indoor Localisation in particular. We also designed a new nonconformity measure with the Weighted K-nearest neighbours as the underlying algorithm.

Our initial results showed very promising location accuracy, compared to other traditional algorithms such as KNN and W-KNN algorithms.

Tracking a Mobile User

Another advantage of CP is that it can work at the on-line learning setting, in which CP learns from its previous predictions and updates its training database to improve the accuracy of the next prediction. This scheme perfectly fits the purpose of tracking an indoor mobile user, since the user cannot jump a long distance in a short period of time, his next movement is within a certain radius of his current position.

A challenge to apply CP, however, is the low latency of many tracking systems. The environment does not always response immediately. Especially with Bluetooth tracking, the expected delay can be as high as 1.28 seconds, and the user might have moved away within that time-frame. This is the case of slow teacher, where immediate feedback cannot be guaranteed. A possible solution is to evaluate the trustworthy of the received signal at a particular moment. For example, a signal received in less than 0.5 second would have higher reliability than that of 2 seconds. We only add the samples with reliable signal measurement into the training database, and discard the others. It is our on-going work to inspect the improvement of CP for tracking a mobile user, considering the user movement speeds.

Fig. 1 .

 1 Fig. 1. Paper progression

Fig. 2 .

 2 Fig. 2. Overlapping signal of 3 stations

Fig. 3 .

 3 Fig. 3. Location Fingerprinting Method

Fig. 4 .

 4 Fig. 4. Location fingerprinting deployment example

Table 1 .

 1 Off-line database snapshot 3-D location (Y) RSSI measurement (X) ... ...... (51, 136, 70) (-57, -41, -62, -59, -86, -81) (51, 146, 70) (-59, -36, -65, -58, -82, -75) (51, 156, 70) (-56, -38, -69, -58, -84, -73) (51, 166, 70) (-63, -32, -73, -57, -87, -68) (51, 176, 70) (-67, -33, -68, -60, -91, -69) ... ......

  Then, another weighted average location is calculated from these K entries, giving an estimated position e dif f = (d x d , d y d , d z d ). Our nonconformity measure is the normal Euclidean distance of the two physical 3-Dimensional locations e same and e dif f α = dist(e same , e dif f ) = (d x s

Fig. 5 .

 5 Fig. 5. Performance of CP with W-KNN and K-NN

Table 2 .

 2 System accuracy at different confidence levels

	Confidence Significance Prediction range Prediction size CP Error
	level	level ε	(Average)	(Average)	(Average)
	90%	0.10	>2.1m	62	7.2%
	87%	0.13	1.5m	29	11.5%
	71%	0.29	3m	1	28.7%
	61%	0.39	null	0	null

Algorithm 1 Conformal Predictor with W-KNN

Input: Training database B= {z1, . . . , z l }, significance level , the new example z l+1 with known RSSI and N is the number of stations. Output: Prediction region R. function dist RSSI(P , Q)