{'fraction_non_alphanumeric': 0.03141424080251506, 'fraction_numerical': 0.015479133168542561, 'mean_word_length': 4.396969304494236, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical Depth based Normalization and Outlier Detection of Gene Expression Data

Normalization and outlier detection belong to the preprocessing of gene expression data. We propose a natural normalization procedure based on statistical data depth which normalizes to the distribution of gene expressions of the most representative gene expression of the group. This differ from the standard method of quantile normalization, based on the coordinate-wise median array that lacks of the well-known properties of the one-dimensional median. The statistical data depth maintains those good properties. Gene expression data are known for containing outliers. Although detecting outlier genes in a given gene expression dataset has been broadly studied, these methodologies do not apply for detecting outlier samples, given the difficulties posed by the high dimensionality but low sample size structure of the data. The standard procedures used for detecting outlier samples are visual and based on dimension reduction techniques; instances are multidimensional scaling and spectral map plots. For detecting outlier genes in a given gene expression dataset, we propose an analytical procedure and based on the Tukey's concept of outlier and the notion of statistical depth, as previous methodologies lead to unassertive and wrongful outliers. We reveal the outliers of four datasets; as a necessary step for further research.


Introduction
The study of gene expressions was first pursued through the use of DNA microarray and, although it is still in use, the technology is being replaced by the RNA sequencing (RNAseq). Both techniques offer the possibility to monitor the behavior patterns of, simultaneously, several thousands of genes. Thus, they allow us to study how the genes function and how they perform under different conditions. The raw gene expression data consists of an array of quantitative data extracted from one or more chips. For microarrays the data is obtained from scanned images of the chips where each value represents a spot or probe luminosity measured on some unknown scale. One of the early observations in gene expression research was that even among chips that were subject to the same processing under the same conditions by the same lab there could be big differences in the gene expressions that were produced by the experiment. These differences tend to be systematic and they reflect a lack of a unified chip to chip measurement scale: the measurement scales are homogeneous within a chip but not between chips. Technological advances in chip design and automation have mitigated these differences but the improvements have not succeeded to eliminate this problem. In RNAseq similar problems occur in the sense that the scale are not homogeneous. In addition another problem in RNAseq is that often the samples contain many zero values and it is standard to delete those samples with too many zero values and keep the rest for analysis. Therefore, in order to transform gene expression data into a format suitable for analysis, we have to preprocess it. In the preprocessing, the raw data is transformed into a scale adequate for analysis. In addition, the systematic sources of variation are removed and the outlier observations as well as the gene expressions that are outliers with respect to the given sample should be found. Due to the particular structure of gene expression data, small samples and high dimensional spaces, standard methods cannot be used. This paper is devoted to two aspects of the preprocessing, doing the so called normalization process and finding the outlier gene expressions. Section 1.2 provides the required background on the normalization process, which makes use of the component-wise median. We discuss in Section 2.1 the lack of appropriateness of the component-wise median as a robust measure of centrality and introduce an improved normalization process obtained by using the concept of statistical data depth. Gene expression data are known for containing outliers even after normalization that are difficult to detect. It is easy to detect severe outliers in individual variables but it is more difficult to detect outliers that may not show up in individual variables but are hidden in the high dimension. Thus, although seeking for outlier genes in a given gene expression is a topic that has been broadly studied, there has not been put much attention on finding gene expression samples that are outliers. In fact, there are only some graphical techniques based on principal components analysis or spectral maps that could give us an idea about the outliers, for instance, the multidimensional scaling graph and the biplot of spectral map analysis Amaratunga, Cabrera and Shkedy (2014); Wouters et al. (2003). None of these two approaches is exact as they are based on a reduction of the dimensionality and outliers are detected only if they appear in such projections as outliers. Section 2.2 presents the selection procedure of the outlier gene expressions in a gene expressions' sample.
Section 3 analyzes the below mentioned real data sets using the methodologies presented in Sections 2.1 and 2.2, for the preprocess of gene expression data. The proposed methodologies and employed datasets have been included in an R-package we have named fdaRNA.

Background on the Datasets
The procedures we propose will be applied to different datasets which are described in what follows. They are the Airway dataset Himes et al. (2014), two sets of the Slc17A5 experiment Raghavan et al. (2007) and Amaratunga, Cabrera and Shkedy (2014) and two known data sets: the called Khan Khan et al. (2001) and Tissue data sets Amaratunga and Cabrera (2004). The Airway dataset is from an RNA-Seq experiment on airway smooth muscle on cell lines that were treated with dexamethasone. The data was extracted from the airway R-package obtained from the Bioconductor repository. It consists of 8 samples, 4 treated and 4 control, and 29391 RNA's that have more than one nonzero values. About half of the RNA's of the original data had two or more samples with a zero value and were discarded. The Slc17A5 data set, which we will refer to as Sialin, compares the gene expression profiles of 6 knocked out mice, mice for which the S1c17A5 gene has been knocked out, versus 6 wild type mice. This is a microarray experiment that was repeated at two time points, one in day first of development during the embryo stage and another at 18 days of age during the adult stage of the mice. The biology suggests that there should be small differences after just a few hours of gestation and larger differences after 18 days. Therefore, the data available for each of the two time points consists of 12 microarrays divided in two groups of six, with 45,101 genes each. Note the large amount of genes and the small number of microarrays. The Khan data set is a well-known dataset consisting of gene expression measurements, obtained using cDNA microarrays, of four types of pediatric small round blue cell tumors: the Ewing family of tumors, Burkitt lymphomas, neuroblastoma and rhabdomyosarcoma. It has 6,567 genes that were filtered in Khan et al. (2001) to 2,308 genes. We make use of the filtered sample. Particularly, we use the 63 microarrays that form the training set, which contains 23 microarrays from the first class, 8 from the second, 12 from the third and 20 from the last one. The Tissue data set consists of the study of the tissue, on mice, using three different treatments. It has 3,487 genes and 41 microarrays. Treatment one and two contain each 11 microarrays and treatment three 19 microarrays. To summarize the information provided on the datasets, see Table1.

Background on Normalization
Normalization was proposed as a way to correct the measurement scale differences of gene expressions. It can be regarded as a sort of calibration process that improves the comparability among gene expressions treated alike. Initially, biologists applied normalization using simple methods. They were linear transformations that corrected the systematic measurement differences by equating either the mean, the sum of each sample, the component-wise median or the component-wise third quartile. The common feature of these normalization schemes was that they assume that the spot intensities on every pair of arrays being normalized are linearly related with no intercept, so that the lack of comparability can be corrected by adjusting every single spot intensity on any gene expression by the same amount, called the normalizing factor, regardless of its intensity level. However, the relationship between the spot intensities of two different gene expressions is usually nonlinear. Thus, an intensity dependent normalization method, i.e. a normalization scheme in which the normalization factor is a function of intensity level, should be preferable to a linear normalization method. In intensity dependent normalization, the transformed spot intensity data is normalized using a nonlinear normalization function. Nonlinear and quantile normalizations were introduced in Amaratunga and Cabrera (2001a) and Amaratunga and Cabrera (2001b), under the label of quantile and non- linear standardization, and in Bolstad et al. (2003), who coined the term "normalization". Recent publications that regard normalization of gene expression data include Wu et al. (2019) that presents an R-package for the evaluation of normalization methodologies, Fundel et al. (2008); Abbas-Aghababazadeh, Li and Fridley (2018) that compares different normalization methods and Amaratunga, Cabrera and Shkedy (2014) for an analysis of the different normalization approaches. Further publications are Libault et al. (2008); Shmulevich and Zhang (2002); Schadt et al. (2001). The different ways to define non-linear normalization have the following commonalities. Let X = {x ij } be a data matrix of gene expressions with i=1,...,G genes and j=1,..,n samples. We denote by X ·j the j−th sample. Suppose that each column of X is subject to slightly different measurement scales. This assumes that there is a data matrix Z = {z ij } such that x ij = h j (z ij ), with the h j s nonlinear monotonic transformations h j : R → R. There are two issues with the estimation of h j 's, the first one is that the h j functions are not identifiable up to a monotonic transformation and the second is the nature of the error term, which one of X or Z is observed with error. The identifiability issue on h j is usually overcomed by establishing an ideal scale or mock array that represents a universal standard. The process is as follows.  1. Let X * be obtained from X by sorting the columns in ascending order. 2. Let X M = (x M 1 , . . . , x M G ) be the vector of row medians of X * . 3. For any x ij , i = 1, . . . , G, j = 1, . . . , n, find the interval, [Q j k ,Q j k+1 ], to which it belongs to. 4. Obtain the normalized value of x ij , by linearly interpolating between the pair of points: (Q M k , Q j k ) and (Q M k+1 , Q j k+1 ). The algorithm for equating all quantiles of the columns of X to the quantiles of the component-wise median array is simpler Amaratunga and Cabrera (2001a): It suffices to substitute 3. and 4. above by: 3. For any value x ij , i = 1, . . . , G, j = 1, . . . , n, find its rank r with respect to the j−th column. 4. The normalized value of x ij is x M r .

Normalization using statistical functional depth
The main issue with quantile normalization is to determine X M . It is well know in multivariate analysis that the coordinate-wise median is not an adequate measure of centrality. To order a data set in R, it is reasonable to follow the decreasing order of the absolute values of the difference between the percentiles of the points in the set and the percentile 50. Thus, the data set central points are the median(s) of the dataset. In contraposition, the component-wise median of a sample of gene expressions is not necessarily similar to any of the gene expressions in the sample. This is because, if our data set is in R p , p > 1, the median has to be calculated coordinate by coordinate. It is well known that the component-wise median may not always give an idea of centrality, and in some cases even not fall inside the convex hull of the data. The following two examples illustrate these cases. Let us take in R 2 the vertices of an equilateral triangle and its center of mass, all with equal probability. By symmetry, the deepest point should be unique and coincide with the center of mass. However, that is not the case of the component-wise median, an area that depends on the selected coordinate axes. For the second example, in Figure 1 it is displayed the convex hull of the vectors (0, 1, 0), (0, 0, 0), (1, 0, 0), (1, 2, 5) and (3, 1, 5) in two different perspectives. The component-wise median of these vectors is (1, 1, 0) and, as it can be appreciate from the figure, if does not fall inside the convex hull of the data. We present an alternative normalization methodology that uses the notion of statistical depth. Tukey (1975) introduced this notion as a way to emulate the behavior of the one-dimensional quantiles, and median, to give a concept of order valid for spaces whose dimension is higher than one. Informally, the data is ordered such that, if a datum is moved toward the center of the data cloud, its statistical depth increases, and, if it is moved toward the outside, its statistical depth decreases. Thus, finding the deepest point in a set is the same as finding the point that is most inside the points cloud or, roughly speaking, the one with most points around it. Tukey's idea was for statistical depth to be a robust measure that would yield a definition of "multivariate median" (deepest element(s)) that was robust, as an alternative to earlier ideas like convex peeling and erosion peeling that lacked resistance to pattern of outliers and therefore had a breakdown point of zero. This work was followed by Donoho and Gasko (1992) and Stahel (1981) who simultaneously introduced the methodology of projection depth and they achieve a "multivariate median" with 50% breakdown point. Later, some requirements that every multidimensional statistical depth should fulfill were established in Zuo and Serfling (2000). Several definitions of multidimensional statistical depth have been proposed. Some of the best-known are Tukey's depth (or half-space depth), Oja's depth and simplicial depth (see, for instance, Liu, Parelius and Singh (1999)). Given a sample of data, to exemplify the Tukey's depth of an element of the sample with respect to the sample, we have represented in Figures 2 and 3 a set of data points by pawns and a king. We aim to compute the depth of the data point represented by the king. Such depth is computed with respect to the data points represented by the pawns and the king itself. Figure 2 exemplifies Tukey's depth in dimension one and Figure 3 in dimension two. Let us first focus on Figure 2. Computing Tukey's depth in dimension one is equivalent to first calculate the minimum between the number of pieces that have to be captured to capture the king if we approach it by the right and if we do so by the left. Then, this number is divided by the total number of pieces. In the plot, we obtain that this number is the minimum between three (from the left) and six (from the right). Then, the depth of the king with respect to the plotted pieces is 3/8. In higher dimensions, p > 1, Tukey's depth is computed by projecting the data in all possible vectors, of norm one, of the sphere of dimension p − 1 (it reduces to a finite amount of vectors when computed with respect to a sample) and computing the minimum of all these resulting one-dimensional depth values. Thus, in Figure 3 it is plotted a yellow line determined by one of the possible vectors and by the data point represented by the king. The procedure consists on projecting each pawn in the yellow line (the projected pawns are the green diamonds at the end of each red segment) and then computing the one-dimensional Tukey's depth. In the plotted line, the one-dimensional Tukey's depth of the king with respect to the plotted pieces is 5/16. A further example is the equilateral triangle and its center of mass commented at the top of this section. In that case, the center of mass obtains uniquely the highest Tukey's depth.  Representation of a data sample in dimension two by a king and pawns, a thick yellow line in which to project the pawns and the projected pawns (green diamonds at the end of the slim red segments). The statistical depth value, restricted to the thick yellow line, of the king with respect to the sample is illustrated with yellow spots. It takes value 5/16.
The objective of gene expression normalization is to equate the scales at which the various gene expressions have been measured. Each scale is a continuous curve that is approximated by a non-decreasing vector that is obtained by ordering the spot intensities of a gene expression. Thus, we propose to choose as X M the deepest gene expression among X * ·1 , . . . , X * ·n , which will preserve the idea of centrality and the simplicity of the procedure while lacking of the cited problematic of the component-wise median. Then, our proposal modifies the normalization procedure of Subsection 1.2.1 substituting item 2. by: 2.* Let X M be the deepest element of the set X * ·1 , . . . , X * ·n . The dimensionality of such vectors could easily reach over half a million because a typical gene expression has 45K-50K genes and for each gene we may observe eleven probes. Thus, classical multidimensional statistical depth functions are not applicable from a practical point of view because of their computational cost Mosler and Hoberg (2006). From a theoretical point of view, in Chakraborty and Chaudhuri (2014) it is studied the non appropriateness of using multidimensional statistical depth functions either, due to the high dependency structure among closed by observations that occur in this type of data. Using parallel coordinates Inselberg (1985), X * ·1 , . . . , X * ·n can be viewed as the partial observations of non-decreasing curves. Consequently, it is appropriate to use statistical functional depths (see Nieto-Reyes and Battey (2016) for a review on the matter), which are not time consuming and lack of the theoretical issues expressed in Chakraborty and Chaudhuri (2014).
Throughout this article, we use the definition of statistical functional depth introduced in Nieto-Reyes (2011). This definition is based on a distance. For its simplicity and interpretability, we propose to use throughout this paper the L 2 space and distance, which coincides with the euclidean because of the use of parallel coordinates. Therefore, what makes this statistical functional depth behave well is the fact that the regions in which the curves are close have less weight in computing the distances among curves than other regions. In addition, this definition of statistical functional depth is robust in the sense that it does not matter how far away is a curve, or a region of the curve, from the curves cloud as it is in the same layer as if it were closer.
The formal definition follows. Given a set of functions, X := {X ·1 , . . . , X ·n } ⊂ L 2 , we aim to compute the depth of X ·j with respect to X for each j = 1, . . . , n. This statistical functional depth computes the statistical functional deepness of each of these functions with respect to the set by constructing a sequence of layer sets, L k , and associated borders, B k . First, it is defined a layer and associated border formed by the most distant functions (or curves).
l=1 B l contains only one curve or is equal to B k0+1 . Thus, the statistical functional deepest curve/s among X ·1 , . . . , X ·n is/are in B k0+1 and B k is empty for k > k 0 + 1. While the statistical functional depth of any general element of L 2 requires the use of these layers, it suffices to make use of the borders to compute the depth of the sample curves. Thus, according to Nieto-Reyes (2011), for any j = 1, . . . , n, the statistical functional depth of X ·j with respect to X is D(X ·j , X) := {k ∈ {1, . . . , k 0 + 1} : X ·j ∈ B k }/n; the index of the border to which X ·j belongs divided by the sample size. Thus, the statistical functional depth of a sample curve is given by the border to which the curve belongs.
To illustrate this statistical functional depth we have included in Figure 4 a zoom, of the central genes, for the ordered genes in the process of normalization of the Khan dataset. With the colors, it is represented the statistical functional depth of the microarray curves, from blue for low statistical functional depth to red for high statistical functional depth, through green and yellow. Thus, it can be appreciated how the border sets are nested in. For a better appreciation, we have included as supplementary material, a movie where it is appreciated the construction of the borders involved in constructing this plot. More involved instances of statistical functional depth with these characteristics can be sought (see Nieto-Reyes and Battey (2016)). However, the intention here is to propose a straightforward, and easy, to understand and apply methodology that does not require of expert statisticians for its implementation and analysis.
Remark 2.0.1. In order to make the component-wise method more usable, the data must be preprocessed by individual linear transformations of each gene expression vector to make them as similar as possible, otherwise the coordinatewise median may yield unreasonable results. If we normalize to the deepest array then this problem may not occur and the linear preprocessing is not necessary. However, in practice, the linear preprocessing may yield a different and perhaps more suitable deepest array, so we believe that is not a bad idea.

Outlier gene expressions
Here we develop an analytical technique to detect outliers specifically design to this kind of data. It is quite simple and so, pleasant for either statisticians and biologist, who has to deal with gene expression data. However, it can be generally used. This technique is based on statistical functional depth and the well-known Tukey's method for outlier detection in R Hoaglin, Mosteller and Tukey (2000). According to Tukey, x ∈ R or y ∈ R is a potential outlier if where q 1 denotes the lower quartile and q 3 denotes the upper quartile of a given sample or distribution. Tukey proposed to take g equal to 1.5 for mild outliers and to 3 for severe outliers. However, later authors propose to take g depending on characteristics of the data such as sample size, see for instance Davies and Gather (1993). By subtracting the two inequalities in (1), the statement is rewritten as x − y > (2g + 1) · (q 3 − q 1 ).
Our methodology orders the data in pairs from the statistical least deep pairs to the statistical deepest. Let X := {X ·1 , . . . , X ·n } be a data matrix of gene expressions, with B k = {X ·k1 , X ·k2 } the k-th least deep pair of gene expressions with respect to X, for k = 2, ..., k 0 . Thus, we propose to generalize Inequality 2 from R to any general space with a metric d(·, ·), by considering the pair X ·k1 , X ·k2 as a potential outlier if d(X ·k1 , X ·k2 ) > G · d(Q 3 , Q 1 ). ( As before, for simplicity, we consider the L 2 (euclidean) distance. For the interquartile range, d(Q 3 , Q 1 ), we estimate it by computing the median{d(X ·k1 , X ·k2 ) : X ·k1 , X ·k2 ∈ B k , k ≤ k 0 + 1}.
If n is odd, d(X ·k0+1 1 , X ·k0+1 2 ) = 0. For an example, let us assume we have the data sample 1.3, 2.1, 2.8, 2.9, 3.2, 3.9, 4.1, 4.8, 4.9, 5.3 ∈ R of size 10. Then, we have to compute the median{4, 2.8, 2, 1.2, 0.7} that is 2, which is the value we would have obtained if computing directly the interquartile range in R of the given sample. To estimate G, we perform a Monte Carlo study in which we draw a multivariate normal distribution with the sample size, dimension and covariance structure of the given data. The latter can be estimated robustly. Here, to estimate G, it is estimated the quantile such that a .01 percent of the columns are flagged as outliers. The process is repeated a number of times, 100 times for the results in Section 3, and the median of those quantiles is the estimated value of G If a pair X ·k1 , X ·k2 satisfies (3), it is in general because one of the elements of the pair is outlying. Thus, we consider as potential outlier the element of the pair that is farther away from the statistical deepest gene expression of X, for instance in L 2 (euclidean) distance. However, a more involve procedure can be undertaken, if the practitioner consider it necessary, in which both elements of the pair would be potential outliers. Note that our approach is different from the one of the functional boxplot, see Dai and Genton (2018) and the references therein.

Result
Through this section we work with a variety of datasets, which we introduced in Section 1: the airway dataset, two sets we named Sialin and two known datasets, Khan and Tissue. Due to the reader might be unfamiliar with he Sialin datasets, we show a representation of them in the step of the normalization process where the columns have being sorted. Thus, in the left column plots of Figure 5 we have displayed on the top row the sorted data after 6 hours of gestation and on the bottom row after 18 days. The central column plots of the figure represent a zoom of the first part of the domain and the right column plots a zoom of a central part of the domain. Making use of the colors of the rainbow, we represent the statistical functional depth of the curves, from blue for low statistical functional depth to red for high statistical functional depth. Once the normalization based on statistical functional depth has being carried out, to show the adequacy of the procedure, we have done a boxplot of the log(· + 1) transform of the resulting datasets, which can be seen in Figure 7. To appreciate the effect of the normalization step, we have included in Figure 6 the boxplot of the log(· + 1) transform of the datasets previous to the normalization step.
In what follows, we compare empirically the proposed analytical procedure, described in Section 2.2, with the two mentioned graphical procedures.

Airway dataset
We apply the analytical proposed procedure on the Airway dataset and obtain that it contains an outlier. This can be observed in Table 4 where each pair of RNA-seq's is displayed following their statistical depth order. Thus, the pair (5, 2) is the statistical least deep one and the pair (8, 7) the statistical deepest one. Below each pair of RNA-seq's, it is shown the L 2 (euclidean) distance between the elements in the pair. A pair contains a potential outlier when this distance is larger than the benchmark for outliers, which is also displayed in the table. For instance, the distance between the elements of the statistical least deep pair of RNA-seq's, 484, 974.7, is just above the benchmark for outliers, which is 484, 923.1. Between the two elements of the statistical least deep pair, the RNA-seq with identification number 5 is the one consider as potential outlier because it has the larger L 2 (euclidean) distance to the average of statistical deepest elements, with identification numbers 8 and 7.  Tables 2 and 3. The y-axis represents the log(· + 1) transform of the counts, for the Airway dataset, and of the intensity levels for the Sialins, Kahn and Tissue datasets.   Table 4 The first two rows are the pairs of RNA-seq's of the Airway dataset ordered from statistical least deep (left) to deepest one (right). The third row contains the distance between the two elements of each corresponding pair, i.e., the distance between RNA-seq's with identification number 5 and 4 is 484, 974.7. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark to detect a pair of gene expressions as outliers. The fifth row is the value of G. The first intra-pair distance and the RNA-seq identifier 5 are in bold to represent it is a potential outlier.
The two graphical procedures are: one in the top left plot of Figure 8, multidimensional scaling; the other, in the top left plot Figure 2 from the Supplementary material, spectral map analysis. They show that the RNA-seq's 5 and 6 are distant from the others; but it is not clear from these graphical procedures where they are distant enough to be consider as outliers. In fact, these two RNA-seq's belong to the two statistical least deep pair of RNA-seq's and each of them belongs to a different class. According to the analytical procedure, RNA-seq 6 is an outlier when the computation of our interquartile distance is restricted to the class it belongs. This is observable from Table 5, from where it can also be seen that RNA-seq 5 is not a potential outlier when restricted to its class.  Table 5 Inside each class, the first two rows are the two pairs of RNA-seq's of the Airway dataset ordered by statistical depth: least deep (left) and deepest (right). The third row contains the distance between the two elements of each corresponding pair. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark to detect a pair of gene expressions as outliers. The fifth row is the value of G from Table 4. The first intra-pair distance of class 2 and the RNA-seq identifier 6 are in bold to represent it is a potential outlier with respect to its class.  Tables 2 and 3. According to this methodology, those far away from the cluster are potential outliers, for instance, in the Airway dataset those RNA-seq identified with number 5 and 6. The colors represent the class to which the identification value belongs: black (class1), red (class 2), green (class 3) and blue (class 4).

Sialin datasets
Using the analytical proposed procedure on each of the Sialin datasets, we obtain that the microarray with identification number 3 is consistently depicted as an outlier. This can be observed in Table 6 where the four less deepest pairs of microarrays of the study after 6 hours of gestation and of the study after 18 days are displayed. Additionally, when we restrict ourselves to each class to compute the distance corresponding to the interquartile range, microarray with identification number 3 is also depicted as an outlier in the study after 6 hours of gestation, see Table 7. Regarding the two graphical procedures, they point to microarray number 3 being a potential outlier in the study after 6 hours of gestation, taking the whole set of data or just the class to which it belongs. However, the graphical procedures do not indicate clearly that it is an outlier in the study after 18 days. Furthermore, the graphical procedures may also indicate that microarray number 8 is a potential outlier in the study after 6 hours of gestation while Table 6 shows that it belongs two the second statistical less deep pair of microarrays with a distance between its elements of 58.1 that is smaller than the benchmark for outliers of 67.6. The multidimensional scaling procedure is displayed in the top row of Figure 8; the central plot is for the study after 6 hours and the right plot for the study after 18 days. Analogously, the spectral map analysis is displayed in the top central and right panels of Figure  2 from the Supplementary material; the left plot is for the study after 6 hours and the right plot for the study after 18 days.  Table 6 The first two rows are the pairs of gene expressions identifiers of the Sialin datasets, where each dataset is ordered from statistical least deep (left) to deepest one (right). The third row contains the distance between the two elements of each corresponding pair. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark for each dataset to detect a pair of gene expressions as outliers. The fifth row is the value of G obtained for each dataset. The first intra-pair distance of each dataset and the gene expression identifiers 3 are in bold to represent they are potential outliers.

Khan dataset
The analytical procedure on the Khan data produces no outlier when the dataset is study as a whole. This can be observed in Table 8 where the eight statistical  Table 7 Inside each class and dataset, the first two rows are the two pairs of gene expressions of the Sialin datasets ordered by statistical depth: least deep (left) and deepest (right). The third row contains the distance between the two elements of each corresponding pair. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark to detect a pair of gene expressions as outliers. The fifth row is the value of G from Table 6. The first intra-pair distance of class 2 of the Sialin dataset after 6 hours of gestation and microarray number 3 are in bold to represent it is a potential outlier with respect to its class.
least deep pairs of microarrays are displayed. In addition, the table contains the benchmarks for outliers and the constant estimated to compute the benchmark. This can also be observed in the graphical procedures, see the bottom left plot of Figure 8 for the multidimensional scaling procedure and the bottom left plot of Figure 2 from the Supplementary material for the spectral map analysis. The 63 microarrays we work with are a training set formed by four classes. Applying the analytical procedure to the first class we obtain two outliers, microarrays number 16 and 4. If we apply the procedure to the second class, we obtain no outlier. However, applying the procedure to the third and the fourth class separately, we obtain one outlier in each, microarray number 32 and 51 respectively. This information is displayed in Table 9 where the two statistical least deep pairs of each of the four classes are displayed emphasizing in bold the outliers. Notice that the outliers are from the least deep pairs of microarrays in their classes. These findings can be seen as coherent with the graphical procedures. For instance, in the multidimensional scaling plot, microarrays 16 and 4 are in the up part of their class and microarrays 32 and 51 are plotted on the top of their respective class.

Tissue
The Tissue data contains one outlier when applied to the whole set of data. This can be seen in Table 10 where we have displayed the eight statistical least deep pairs of microarrays and the statistical least deep pair has a distance of 38.7 between its two elements, that is larger than the outlier's benchmark of 38.3. In the graphical procedures, we can observe that microarray 34 is the one with highest OY-axis in the multidimensional scaling plot, see the bottom right plot of Figure 8 for the multidimensional scaling procedure and the bottom right plot of Figure 3 of the Supplementary material for the spectral map analysis.  Table 9 Inside each class, the first two rows are the two pairs of gene expressions of the Khan dataset ordered by statistical depth with the least deep on the left and the second least deep on the right. The third row contains the distance between the two elements of each corresponding pair. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark to detect a pair of gene expressions as outliers. The fifth row is the value of G from Table 8.
There are four intra-pair distances and four corresponding gene expressions in bold to represent i potential outliers with respect to their class.
However, it would not be clear whether to flag it as a potential outlier. It would be more clear, however, to point to microarray number 20 as a potential outlier in its class. However, although the analytical procedure contains it in the least deep microarray of its class, it does not detect it as potential outlier. In fact, the analytical procedure, detects no outlier when the interquartile distance is computed for each class separately, see Tables 11.

A comparison with Robust Multi-array Average
Robust multi-array average, generally known as RMA, is a conventional preprocessing technique that makes use of quantile normalization. We develop here a corresponding technique that makes use of our statistical functional depth based normalization, and refer to it as Functional Depth Normalization (FDN).
We compare both through simulations in what follows.

Table 10
The first two rows are the eight statistical less deep pairs of gene expressions of the Tissie dataset, with the statistical least deep on the left. The third row contains the distance between the two elements of each corresponding pair. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark of the dataset to detect whether each pair of gene expressions are potential outliers. The fifth row is the value of G obtained for this dataset. In bold it is the first intra-pair distance and the corresponding gene expression that is farther away in L 2 (euclidean) distance from the deepest gene expression.

Table 11
Inside each class, the first two rows are the two pairs of gene expressions of the Tissue dataset ordered by statistical depthwith the least deep on the left and the second least deep on the right. The third row contains the distance between the two elements of each corresponding pair. The fourth row is the value G · d(Q 3 , Q 1 ) of Inequality 3, the benchmark to detect a pair of gene expressions as outliers. The fifth row is the value of G from Table  10.

Simulations
In order to compare FDN with RMA, we performed a simulation involving 12 samples with 1000 genes and 11000 probes per sample. Each probe was generated from a t-distribution center at 3. We performed the procedure three times for the t-distribution with 10, 5 and 2 degrees of freedom. The few negative values were replaced by a small positive number. Then, the first 11000 probes of the first six samples were shifted by a positive δ taking values δ = 0, 0.25, 0.5, 1, 2 to create a difference in 100 genes. The other 900 genes should show no difference. Finally the resulting probe values for each of the 12 samples were elevated to the power of 3 + where the 12 epsilon values where sampled from a uniform(0,2) distribution. This last step applies a non-linear scale change to each sample separately, to create the need for normalization. We generated 100 datasets following this method and we calculated gene expression sets using RMA including full median polish, not one-step median polish as was performed in the original RMA Amaratunga and Cabrera (2001a); Irizarry et al. (2003). The normalization step was performed using quantile normalization (RMA) or with statistical functional depth normalization (FDN). To study the performance of both methods we applied a conditional t-test Amaratunga and Cabrera (2009) to compare the gene expressions between the two groups. Recall that the first 100 genes were differentially expressed by δ and the other 900 were not. The results are given in Table 12 columns 1 and 2 showing the number of true discoveries out the 100, which is equivalent to the power. It appears that for this example the FDN technique produces slightly better power with less false discoveries (see columns 4 and 5 of Table 12), but they are not very different. In Table 12 column 3, and 6, we present the result of the simulation using FDN with Tukey's biweight M-estimator replacing median polish, for comparing FDN with Median Polish. Table 12 shows that for this simulation the power of the conditional t-test for FDN with the M-estimator is better than the power for FDN with median polish. In conclusion this simulation shows situations where FDN together with Tukey's biweight M-estimator are good choices for implementing an alternative RMA algorithm.

Table 12
Power results and false discoveries out of 900 of the conditional t-test for RMA (columns 1 and 4), FDN with median polish (column 2 and 5) and FDN with M-Estimator (column 3 and 6). The best result of the three in each case is in bold blue and the second best in blue.

Discussion and conclusion
The propose methodology allows to perform a robust normalization procedure that satisfy the appropriate mathematical properties and an outlier detection methodology that is unambiguous, no depending on the expert pursuing it.
Our research has found the following pressure findings: • The Airway dataset contains one outlier and when restricting two each class separately, another outlier is found. The analytical procedures would have suggested, however, both RNA-seq's as outliers of the whole dataset. • The Sialin dataset after 6 hours of gestation detects the same outlier when studying the dataset as a whole and when studying it by classes. The graphical procedures suggest the same outlier and one additional outlier. • The Sialin dataset after 18 days contains one outlier when the data is study as a whole which is not spotted by the graphical procedures. Like the Tissue data, below, it contains no outlier when the data is studied by classes. • The Khan dataset has no outlier when studying the dataset as a whole.
However, it contains outliers in three of its four clases when studying them separately. In fact, one of the clases contains two outliers. • In contraposition with the Khan dataset, the Tissue dataset contains an outlier when studying the dataset as a whole but no outliers when studying the classes separately. The analytical procedure differs from the suggestion of the graphical procedures. For instance, the multidimensional scaling producere would suggest an outlier in the second class. See