Bioinformatics, 37(16), 2021, 2365–2373 doi: 10.1093/bioinformatics/btab105 Advance Access Publication Date: 20 February 2021 Original Paper Gene expression CuBlock: a cross-platform normalization method for gene-expression microarrays 1,2 Valentin Junet , Judith Farrés1, José M. Mas1,* and Xavier Daura 2,3, * 1 Anaxomics Biotech SL, Barcelona 08008, Spain, 2Institute of Biotechnology and Biomedicine, Universitat Autònoma de Barcelona, Barcelona 08193, Spain and 3Catalan Institution for Research and Advanced Studies (ICREA), Barcelona 08010, Spain *To whom correspondence should be addressed. Associate Editor: Jonathan Wren Received on July 14, 2020; revised on February 4, 2021; editorial decision on February 10, 2021; accepted on February 16, 2021 Abstract Motivation: Cross-(multi)platform normalization of gene-expression microarray data remains an unresolved issue. Despite the existence of several algorithms, they are either constrained by the need to normalize all samples of all platforms together, compromising scalability and reuse, by adherence to the platforms of a specific provider, or sim- ply by poor performance. In addition, many of the methods presented in the literature have not been specifically tested against multi-platform data and/or other methods applicable in this context. Thus, we set out to develop a nor- malization algorithm appropriate for gene-expression studies based on multiple, potentially large microarray sets collected along multiple platforms and at different times, applicable in systematic studies aimed at extracting know- ledge from the wealth of microarray data available in public repositories; for example, for the extraction of Real- World Data to complement data from Randomized Controlled Trials. Our main focus or criterion for performance was on the capacity of the algorithm to properly separate samples from different biological groups. Results: We present CuBlock, an algorithm addressing this objective, together with a strategy to validate cross- platform normalization methods. To validate the algorithm and benchmark it against existing methods, we used two distinct datasets, one specifically generated for testing and standardization purposes and one from an actual experi- mental study. Using these datasets, we benchmarked CuBlock against ComBat (Johnson et al., 2007), UPC (Piccolo et al., 2013), YuGene (Lê Cao et al., 2014), DBNorm (Meng et al., 2017), Shambhala (Borisov et al., 2019) and a simple log2 transform as reference. We note that many other popular normalization methods are not applicable in this con- text. CuBlock was the only algorithm in this group that could always and clearly differentiate the underlying biologic- al groups after mixing the data, from up to six different platforms in this study. Availability and implementation: CuBlock can be downloaded from https://www.mathworks.com/matlabcentral/fil eexchange/77882-cublock. Contact: xouse@anaxomics.com or xavier.daura@uab.cat Supplementary information: Supplementary data are available at Bioinformatics online. 1 Introduction (RWD)—data derived from a number of sources, outside the context of Randomized Controlled Trials (RCTs), and associated with out- Since the first whole-genome microarray study of gene expression comes in an heterogeneous patient population (Berger et al., was published in 1997 (Lashkari et al., 1997; Schena et al., 1995), 2017)—may in fact boost the use of the wealth of available micro- high-throughput gene-expression microarrays have been a standard array data in the near future. The importance of RWD as a comple- in many experimental designs in biological and biomedical research. mentary information source in drug-evaluation studies is based on Although their use is being replaced by next-generation sequencing the observation that data from RCTs does not always match results techniques such as RNA-Seq (Nagalakshmi et al., 2008), the large from observational studies (Trotta, 2012), mostly owing to the lim- amounts of microarray data relevant to an equally large variety of ited number of RCT patients, their over-monitoring and the limited biological and biomedical problems and available in public data- follow-up time. Thus, certain adverse drug reactions or lack-of- bases constitutes a valuable resource that will remain in use for efficacy problems are hidden until RWD studies are performed, and many years. The potentiality of resources such as the Gene drug administrations have come to encourage the extraction of in- Expression Omnibus (GEO) as sources of Real-World Data formation from sources complementary to RCTs to increase C The Author(s) 2021. Published by Oxford University Press. V 2365 This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com 2366 V.Junet et al. evidence around treatments (U.S. Food and Drug Administration, normalization methods for the integration of gene-expression 2018; Sherman et al., 2017). A problem of RWD is that it tends to microarrays. be highly heterogeneous, thus requiring careful analysis and statis- One should also note that the methods mentioned above were tical treatment (Bartlett et al., 2019; Berger et al., 2017). Translated not necessarily developed with a same purpose. For example, to the context of this study, microarray data relevant to a particular ComBat was developed for the adjustment of batch effects and is problem will often originate from different laboratories and experi- often used in combination with other methods in cross-platform ments, possibly using different microarray platforms (Bumgarner, normalization procedures; although it does not normalize 2013) and almost certainly obtained in different batches. In order to platforms separately, we introduced it in this study because of its make a sensible use of such a heterogeneously sourced data, a data- broad use. Thus, the SCAN/UPC package offers the possibility to normalization step is required before data analysis. Normalization apply ComBat after SCAN normalization and summarization at can be relatively straightforward when dealing with different gene-level—and before transformation to UPCs if so chosen. Thus, batches of a same experiment using the same biological samples, methods like ComBat are often called integration methods, in so platform and operator, but gets increasingly complex as different that they integrate previously normalized data. Nevertheless, operators, platforms and sample sources are introduced. This often ComBat has become, on its own, a popular first choice for cross- leads to studies discarding part of the available data, which could platform normalization and a frequent benchmark standard for otherwise be used to increase the chances of discovery of meaningful other methods (Irigoyen et al., 2018; Walsh et al., 2015). patterns or improve their statistics. Shambhala, on the other hand, is classified as a harmonization A main source for sample differences arising from systematic method because it uses a golden standard as reference, while biases is the mixing of data from different microarray platforms. YuGene is a transformation and XPN, DBNorm and SCAN are referred to as normalization methods. In this study, like in Rudy and Unfortunately, most standard and widely used normalization meth- Valafar (2011), we will refer to all these methods as cross-platform ods are applicable to or have been developed for the single- normalization methods so far as they are being used in the literature microarray-platform context (Rudy and Valafar, 2011), making to make data across different platforms comparable for the purpose them generally inappropriate for the cross-study analysis of existing of analysis. datasets. On the other hand, most existing cross-platform normal- Although the number of normalization methods proposed in the ization methods, such as ComBat (Johnson et al., 2007), XPN literature is large, to our knowledge there are no other major cross- (Shabalin et al., 2008) or DWD (Benito et al., 2004), require the platform normalization methods that can be applied to gene- data from different platforms to be normalized together—XPN and expression microarrays in a platform agnostic way and that have DWD were in fact developed for pairwise cross-platform normaliza- been tested and validated as such. To enable systematic studies tion. For large datasets, normalizing platforms together can be re- involving the download of microarray data from databases (possibly strictive. In addition to involving the normalization of a large joined at different times) and its normalization and storage for later re- dataset, the eventual addition of new microarray data requires glo- trieval, allowing a non-linear use of the data—for example, in suc- bal renormalization. This led to the more recent development of cessive analyses incorporating different amounts of data as available sample-wise, cross-platform normalization methods such as SCAN or necessary, it is essential that a downloaded microarray set need (Piccolo et al., 2012) and UPC (Piccolo et al., 2013), YuGene (Lê not be normalized more than once. Here, we introduce a novel Cao et al., 2014), DBNorm (Meng et al., 2017)—which can operate cross-platform normalization method fulfilling all these conditions. sample or platform wise—and Shambhala (Borisov et al., 2019). The algorithm is called CuBlock, which stands for Cubic approxi- SCAN performs a sample-wise normalization assuming a double mation by Block. We validate its performance using various metrics Gaussian mixture distribution. It was, however, specifically and compare it to six methods that can be used in a cross-platform designed for Affymetrix and two-channel Agilent platforms, thereby context, namely, the log2 transform of raw data, ComBat, YuGene, restricting its general use. Developed by the same authors, the DBNorm, Shambhala and UPC. Overall, CuBlock shows the best Universal exPression Code (UPC) builds on SCAN to generate stand- performance in this group. ardized estimates of expression that have a consistent interpretation across platforms, measuring how much the expression of the gene deviates from model-estimated background levels within the sample. 2 Materials and methods Although, the derivation of UPCs is platform specific and, to our In this section, we introduce the datasets used for the validation of knowledge, it is currently available only for Affymetrix and Agilent CuBlock and describe the data preprocessing approach and the microarrays through the package SCAN/UPC, the program offers methods used for benchmarking and validation. also a generic UPC function applicable to any microarray platform, by making assumptions on the background and background-plus- signal distributions. The other distribution-based normalization 2.1 The datasets method, DBNorm, scales the data distributions from the individual We selected two benchmark datasets previously used in similar stud- microarrays to a common form, which does not need to be predeter- ies (Borisov et al., 2019; Rudy and Valafar, 2011). The first set (here mined (e.g. the distribution from a reference microarray). As a called the reference dataset) originates from projects MAQC downside, it is very slow. On the other end, YuGene uses a simple (MAQC-I) (MAQC Consortium, 2006) and SEQC/MAQC-III transform that assigns a modified cumulative proportion value to (SEQC/MAQC-III Consortium, 2014), which made use of reference each measurement, making the normalization very fast. Finally, RNA samples to assess repeatability of gene-expression microarray Shambhala uses a harmonization method that transforms each pro- data within a specific site, reproducibility across multiple sites and file so that it approaches the output of a chosen golden-standard comparability across multiple platforms. The second set (here called platform. the experimental dataset) originates from a study trying to assess Some methods like UPC, Shambhala and MatchMixeR (Zhang profile differences of human spermatozoal transcripts from fertile et al., 2020) have specifically included in their design the possibility and teratozoospermic males (Platts et al., 2007). The use of these to integrate data from both microarray and RNA-seq sources. As a two datasets allows us to assess, independently, effects from tech- matter of fact, any method that can be applied at the gene level nical replicates (same biosample, analyzed in different labs with could be adapted for such studies. It should be kept in mind, how- repetition) and biological replicates (different biosamples corre- ever, that this requires source-specific preprocessing steps taking sponding to a same condition). into account the fundamental differences between these two types of data. In fact, the conceptual differences between microarray and 2.1.1 The reference dataset RNA-seq gene-expression measurements are so significant that they The data of this set are accessible in GEO with accession numbers may require distinct normalization procedures (Rapaport et al., GSE5350 (MAQC-I) and GSE56457 (MAQC-III), respectively. The 2013). In this study, we will solely focus on the application of set contains microarray gene-expression data corresponding to four CuBlock 2367 titration pools from two distinct reference RNA samples: (A) 2.2 Data processing Stratagene’s Universal Human Reference RNA pool; (B) Ambion’s To make the analysis as platform agnostic as possible, we took the Human Brain Reference RNA pool; (C) pool with an A:B ratio of image-processed raw intensities for all non-control probes and disre- 3:1; (D) pool with an A:B ratio of 1:3. These biosamples had been garded any platform-dependent background-signal correction such analyzed using different platforms and in different sites, as described as that provided by mismatch probes in Affymetrix platforms. CEL (MAQC Consortium, 2006; SEQC/MAQC-III Consortium, 2014). files for Affymetrix and txt files for the other platforms were used. Following the work from Rudy and Valafar (2011) and Borisov Probes with invalid intensities (NaN) in datasets HT12 and HUG et al. (2019), we selected data from six of the platforms (between were ignored. For Affymetrix microarrays, the intensities of probes parentheses, data-set identifier in this study, GEO platform ID and constituting a probe set were averaged. We note that in the context project of origin): of this study preprocessing of Affymetrix probe sets by simple aver- aging performed just as well as more complex treatments including • Affymetrix Human Genome U133 Plus 2.0 Array (AFX, background correction and RMA median polish (Irizarry et al., GPL570, MAQC-I): three experiments (sites) (AFX_1 to 2003), which was not completely unexpected—for example, AFX_3), with four biosamples (A–D) per experiment and five Hubbell et al. (2002) found that simple averaging performed com- replicates per biosample (60 samples) parably to more robust approaches from the Affymetrix Micro • Agilent-012391 Whole Human Genome Oligo Microarray Array Suite under low-noise conditions. From this point onward, the Affymetrix probe-set-average intensities were treated the same G4112A (AG1, GPL1708, MAQC-I): three experiments (AG1_1 way as the raw probe intensities from the other platforms (the num- to AG1_3), with four biosamples (A–D) per experiment and five ber of Affymetrix probe sets per gene being on the same order as the replicates per biosample (60 samples) number of probes per gene in other platforms). Note that we could • Illumina Sentrix Human-6 Expression BeadChip (ILM, have skipped the probe-set-averaging step and worked directly with GPL2507, MAQC-I): three experiments (ILM_1 to ILM_3), all Affymetrix probes. While we tried this, it increased significantly with four biosamples (A–D) per experiment and five replicates the computational cost of the normalization procedure (due to the per biosample (59 valid samples) several-fold increase in number of probes) at a marginal gain. At this • point, probe-set average intensities (Affymetrix) and raw probe Illumina HumanHT-12 V4.0 Expression Beadchip (HT12, intensities (other platforms) were log2 transformed—we also tried GPL10558, MAQC-III): two experiments (ILM_COH and applying quantile normalization in the preprocessing without it ILM_UTS), with four biosamples (A–D) per experiment and improving significantly the results. CuBlock normalization was then three replicates per biosample (24 samples) applied for each platform separately. Since probes vary among the • TM different microarray platforms, the normalized datasets were then GeneChipV R PrimeView Human Gene Expression Array (PRV, GPL16043, MAQC-III): one experiment (AFX_USF_PRV), with transformed from the probe level to the protein level by mapping probes to UniProtKB accession numbers (ACs) and keeping only four biosamples (A–D) and four replicates per biosample (16 those probes that map to an AC present in all platforms. The map- samples) ping to proteins was performed by taking the gene identifiers from • Affymetrix Human Gene 2.0 ST Array (HUG, GPL17930, the GEO tables containing the microarray data. Each selected AC MAQC-III): one experiment (AFX_USF_HUG), with four was then assigned an intensity equal to the average of the normal- biosamples (A–D) and four replicates per biosample (16 ized intensities of associated probes in the given microarray. The samples) choice of UniProt ACs, rather than gene identifiers, was made to fa- cilitate streamlining with protein-level post-normalization analysis in studies where microarray data is used to infer protein expression Note that in the MAQC-I study the following microarrays (i.e. probes matching CDS regions). We note, however, that from AG1 were discarded as outliers after the Agilent’s Feature CuBlock delivers normalized data at the probe level, meaning that Extraction QC Report: AG1_1_A1, AG1_2_A3, AG1_2_D2, the user is free to summarize the data at the gene level if appropriate, AG1_3_B3. Since the data for these microarrays is nevertheless and that potential information regarding, for example, probes deposited and we wanted our analysis to be as independent as matching non-coding exon regions, will remain available. possible of platform-dependent data-preprocessing steps, we To benchmark CuBlock against established normalization meth- considered also their inclusion. To this end, we evaluated the ods applicable in a generic cross-platform context, we compared it correlation of the data between all AG1 samples and observed to a simple log2 transform and to the methods ComBat (Johnson that the ‘outliers’ are highly correlated to the non-outliers of et al., 2007), YuGene (Lê Cao et al., 2014), DBNorm (Meng et al., the same experiment and of the other two experiments (about 2017), Shambhala (Borisov et al., 2019) and UPC (Piccolo et al., 0.97 in both cases). A dimension reduction of the raw data 2013). YuGene, DBNorm and UPC were applied following the same showed also no outliers. We therefore decided to include these procedure used for CuBlock, i.e. normalization of the log2 transform four microarrays in the dataset. of the probe intensities and successive mapping to ACs. ComBat requires all microarrays to be normalized together, which implies 2.1.2 The experimental dataset their merging before normalization. Therefore, in this case the map- This dataset contains spermatozoal RNA samples from ping to UniProtKB ACs and selection of ACs present in the different normally fertile (N) and heterogeneously teratozoospermic (T) platforms was performed after log2 transform and before ComBat subjects and is accessible in GEO with accession number GSE6969. normalization. We note that DBNorm allows normalization per The samples had been analyzed on three different platforms sample and per platform. We performed both, but show only the (between parentheses, data-set identifier in this study and GEO results obtained with sample-wise normalization since they are bet- platform ID): ter. Comparison to Shambhala was done only for the datasets AFX, AG1 and ILM from the reference dataset, since Shambhala- • Affymetrix Human Genome U133 Plus 2.0 Array (AFF, normalized data for these sets has been already reported by the GPL570): 13 independent biosamples of type N and 8 of type T authors as supplementary data to Borisov et al. (2019). DBNorm • Illumina Sentrix Human-6 Expression BeadChip (ILL1, was only used on the experimental dataset, as the calculations GPL2507): five independent biosamples of type N and eight of turned out to be forbiddingly slow. To perform the calculations we used the R package sva for ComBat (https://bioconductor.org/pack type T. All ILL1 biosamples are replicates of AFF biosamples. ages/release/bioc/html/sva.html) and the packages provided by • Illumina Sentrix HumanRef-8 Expression BeadChip (ILL2, YuGene (https://cran.r-project.org/web/packages/YuGene/index. GPL2700): four independent biosamples of type N and 6 of type T html) and DBNorm (https://github.com/mengqinxue/dbnorm) authors in the respective papers. Calculations with these programs 2368 V.Junet et al. were performed with default settings. For DBNorm, in order to re- the number of neighbors is called perplexity, and is typically given produce the general case (e.g. this study), in which a reference values between 5 and 50. Due to its stochastic nature and the de- microarray cannot be straightforwardly selected, we used the pendence on the chosen perplexity parameter, the algorithm may option of normalization into a normal distribution. For UPC converge to irrelevant solutions. We thus performed 10 runs for (https://www.bioconductor.org/packages/release/bioc/html/SCAN. each of a number of perplexity values and selected the one produc- UPC.html), we used its generic function for expression set, which ing the most consistent biological partitioning according to the aver- takes expression  samples data, with default parameters from the age silhouette values. For the reference dataset we used perplexity package. values from 5 to 50, in increments of 5, and selected the representa- tion giving the best clustering relative to sets A, B, C and D. For the experimental dataset, we used perplexity values 5, 10 and 15 and 2.3 Comparison and validation methods selected the representation giving the best clustering relative to sets To validate and compare the cross-platform normalization methods T and N. The MATLAB function tsne was used to perform the t- evaluated in this study we used the methodology described below. SNE dimension reduction. The objective was to increase the sensitivity, i.e. the identification of true biological differences, while minimizing platform and various kinds of replica effects. All validation methods were applied on a 2.3.3 Dendrogram subset of 500 proteins that best distinguish two given biological We performed a hierarchical clustering analysis using the Euclidean groups. distance as metric and the arithmetic mean as linkage criterion, and To select the 500 proteins we first performed a differential ana- represented the resulting cluster hierarchy as a dendrogram. To as- lysis on the normalized data, for all platforms in the reference or ex- sess the significance of the clusters, we applied multiscale bootstrap perimental dataset. To this end, we performed Welch’s t-test to resampling as provided in the R package pvclust (Suzuki and evaluate, for each protein, the difference between the associated Shimodaira, 2006). By default, this package considers 10 relative mean intensities in units of uncertainty (the t-statistic) in two bio- bootstrap sample sizes (bootstrap sample size divided by total sam- logical groups, A and B from the reference dataset (total of 16 624 ple size), from 0.5 to 1.4, with 1000 resamplings per sample size, proteins) or N and T from the experimental dataset (total of 16 937 leading to a total of 10 000 bootstrap resamples. The package pro- proteins). Note that we deliberately avoid considerations on vides two statistics to estimate the significance of the obtained clus- whether the datasets meet the requirements of the t-test, since we ters: the bootstrap probability (BP) or frequency (expressed as used the test simply to identify the 500 proteins with largest separ- percentage) of observation of a given cluster in the bootstrap resam- ation of group means per uncertainty unit, that is, with lowest asso- ples, and the approximately unbiased P-value (AU), an unbiased ver- ciated P-values, irrespective of the error in the P-value and, sion of BP. More details on multiscale bootstrap resampling can be therefore, of its valid interpretation as a probability. Although for found in Shimodaira (2004). We plot the dendrograms using the R such purpose we do not require the calculation of FDR-adjusted P- package dendextend. values (q-values) (Storey, 2002), since they conserve P-value rank- ing, we did obtain them and show corresponding ROC-like curves 2.3.4 SVM classification (the cumulative distribution function of the q-values) in The goal of this analysis was to assert whether relevant patterns can Supplementary Figure S1 (Supplementary Information). We decided be found using the data from only one platform. Support vector to select a fixed number of proteins, rather than proteins with a P- machines (SVM) (Cortes and Vapnik, 1995) are binary classifiers or q-value below a given arbitrary threshold, to enable the compari- applicable to problems that are reducible to a binary outcome, such son of methods using datasets of equal and reasonably large dimen- as the T and N phenotypes in our experimental dataset. We trained sionality. We note, nevertheless, that the 500-protein cut a linear support vector machine model for each platform using the corresponds to an FDR well below 10–2 (Supplementary Fig. S1). following approach. To reduce feature-vector dimensionality, where The differential analysis was performed with the MATLAB function dimensions are proteins (more specifically their microarray-derived mattest. intensities), while retaining the capacity to asses how well the 500 proteins separate the T and N populations, the training was per- 2.3.1 Silhouette plot formed six times, starting with dimension 5 and increasing it up to Silhouette plots are graphical displays of data partitions dimension 10. For each of 1000 runs with a given dimensionality, (Rousseeuw, 1987), where clusters are represented by so-called sil- we selected randomly from the 500 protein set as many proteins as houettes generated by comparison of cluster tightness and separ- dimensions, extracted the corresponding data from sets T and N, ation. The method assigns a silhouette value between -1 and 1 to trained a linear SVM model for each platform, separately, and tested each element of a cluster, indicating if the element is well clustered it on the other two platforms. This led to a total of 6000 models per (value close to 1), lies between two or more clusters (close to 0) or is platform. For each platform, we calculated the mean and standard likely misclassified (close to -1). The silhouette plot is then generated deviation of different classification scores over the 6000 models, by representing the values for all elements as bars, for the different namely, Accuracy, Matthews Correlation Coefficient (MCC) cluster partitions. We computed three silhouette plots for the refer- (Boughorbel et al., 2017), Balanced Accuracy and Area Under the ence dataset: one identifying clusters with platforms, one where the ROC Curve (AUC) (Fawcett, 2006). The MATLAB function svm data was assigned to groups A [ C and B [ D and one where the was used to train the SVM models. partitioning was represented by sets A, B, C and D. For the experi- mental dataset, silhouette plots based on platform partitioning and T versus N partitioning were computed. The MATLAB function sil- 3 Algorithm houette was used to compute the silhouette plots. The CuBlock algorithm relies on the simple and widely used as- sumption that most genes are neither over- nor under-expressed 2.3.2 t-SNE dimension reduction (Yang et al., 2002). Thus, a transformation following the cubic poly- t-SNE (Maaten and Hinton, 2008) is a stochastic dimension- nomial x3 will leave most of the genes around 0 and slowly differen- reduction method aimed to preserve the local structure of data tiate the extremes, i.e. the under- and over-expressed genes. The (keeping the low-dimensional representation of very similar data complementary idea in CuBlock is the use of a block-wise trans- points close together) while retaining essential traits of global struc- formation, which had been already implemented successfully in ture. It analyzes the neighborhood of the data points by calculating XPN (Shabalin et al., 2008). To this end, CuBlock partitions probes pairwise conditional probabilities representing their similarity. The into clusters and, for each sample and probe cluster (i.e. for each method then tries to find a low-dimensional representation that min- data block), transforms the data by a procedure that involves its imizes the difference between the high-dimensional and low- mapping to objective values between -1 and 1 (with density increas- dimensional conditional probabilities. The parameter controlling ing toward 0) and the fitting of a cubic polynomial to the resulting CuBlock 2369 Fig. 2. Example histograms for samples from different platforms. (A–C) Histograms before CuBlock normalization and after log2 transform. (D–F) Histograms after CuBlock normalization. A, D: biosample A, platform AFX; B, E: biosample B, plat- form AG1; C, F: biosample A, platform AG1 Although we initially determine the probe clusters using all samples, we then normalize sample by sample to reduce the dependence of the normalization on the full sample collection. The strategy of nor- malization by block is similar to that used by the cross-platform nor- malization method XPN (Shabalin et al., 2008). However, XPN defines probe clusters and sample clusters with two independent ap- plication of k-means (one on the input matrix and the other on its Fig. 1. Pseudocode describing the CuBlock algorithm (see description in Section 3) transpose) and blocks are then constituted by all possible combina- tions of one sample cluster with one probe cluster. distribution (see below). By using data blocks, different profiles pre- For each block, and each of the N repetitions of the k-means sent (mixed) in the full dataset are considered, and as many different clustering, CuBlock fits a cubic polynomial to a mapped set of cubic polynomials are fitted to them, underlying different shapes points symmetrically distributed between -1 and 1 with density contained within the original distribution. A pseudo code of the increasing toward zero. This is performed in four steps, as shown in CuBlock algorithm is described in Figure 1. It calls two additional Figure 1. First, the block data is linearly transformed to z-scores algorithms with pseudo codes provided in Supplementary Figures S2 (zero mean and unit standard deviation). These are then used as in- and S3 (Supplementary Information). put values of a mapping function whose output values will be used The input to CuBlock is a matrix X containing the log2 trans- to fit the cubic polynomial, as described in the pseudocode shown in form of the gene-expression microarray intensities, where columns Supplementary Figure S2. The mapping associates the sorted values are samples and rows are probes. As discussed in Section 2.2 it is up present in the block to an equal number of equidistant points be- to the user to decide any level of preprocessing of the input log2- tween -1 and 1, and takes these new points to an uneven power in transformed intensities, for example, probe-set summarization for order to have their distance decrease as they approach zero from ei- Affymetrix microarrays. ther side (Supplementary Fig. S4). The exact uneven power will de- CuBlock makes use of the k-means clustering algorithm (Lloyd, termine how slow is the growth of the points around zero, and is 1982) to partition probes in the space defined by the samples—a selected such that, on average, the values of the block that are within probe data point is a vector of probe intensities of dimension equal standard deviation, i.e. the block values between -1 and 1, are to the number of samples—and is applied per platform, i.e. the k- mapped to a value smaller than 0.1 (Supplementary Fig. S5). The al- means clustering is performed for all samples of a given platform. k- gorithm tries uneven powers between 3 and 21 and the first one that means is an iterative algorithm that tries to partition the data into a fulfills the criterion is selected. Next, the algorithm finds the coeffi- predefined number k of non-overlapping clusters, starting from a cients of a cubic polynomial that, when evaluated on the sorted random initialization of their centroids. Because of the random ini- block data (input values), best fits the output values from the map- tialization, clusters from different runs may differ and the core part ping function (Supplementary Fig. S4). We chose to fit a cubic poly- of the CuBlock algorithm is repeated several times for different solu- nomial instead of a higher degree one to avoid overfitting. tions of k-means. Thus, input parameters k and N in Figure 1 refer Polynomial coefficients were obtained with the MATLAB function to the chosen number of k-means clusters and repetitions, which in polyfit, with degree 3. this work took values of 5 and 30, respectively. The number of clus- If the block data is not symmetric or contains many outliers, a ters was chosen to be low enough that the k-means algorithm will cubic polynomial will produce a poor fit. Thus, the polynomial will not, for some partitions, converge always to the same solution and increase along the symmetric part of the block and decrease as it high enough that blocks with different distributions will be reaches the outliers (Supplementary Fig. S6). Despite leading to a obtained. The CuBlock algorithm finds first a probe-cluster partition poor fit, this feature can be used to identify asymmetry issues and in the space defined by the samples and then applies its normaliza- outliers. When this is the case and decreasing values are identified tion scheme to data blocks defined as those (log2) probe-intensity after evaluating the polynomial on the block data, the decreasing values from a sample that belong to a given cluster. Therefore, for k values are corrected in order to preserve data sorting upon normal- clusters and m samples we have a total of k  m blocks. The advan- ization. Roughly, the correction equates the decreasing values to the tage of the normalization by block is that it decomposes the distribu- last increasing value (after an increasing section) or to the last tion of probe intensities of a sample into its different block decreasing value (before an increasing section). The precise correc- distributions, according to similarities between probes found by the tions are described in Supplementary Figure S6, and the cases where clustering algorithm in the space of all samples. These different dis- the cubic polynomial might decrease are considered in tributions will enable the emergence of different patterns present in Supplementary Figure S3. the data. Instead, if the blocks were selected at random or the whole The output of CuBlock is a matrix of normalized gene- sample was used, the normalization method would estimate parame- expression values, where columns are samples and rows are probes. ters based on a unique distribution, masking these different patterns. As discussed in Section 2.2 it is up to the user to decide at which 2370 V.Junet et al. level and using which database codes for the mapping to that level, number GSE65212 and involving five biological groups (Maire the probe values should be summarized. et al., 2013a,b; Maubant et al., 2015), demonstrating that results Figure 2 shows the histograms of different samples before and from CuBlock are very consistent with those from RMA after normalization. While before normalization the samples normalization. follow clearly different distributions, after normalization the distributions are much more homogeneous. We note that before nor- 4.1 Reference dataset: six platforms malization the distributions are clearly platform dependent (com- Figure 3 and Supplementary Figures S8–S11 show the results pare A and C, which correspond to the same biosample but different obtained with the different normalization methods using the dendro- platforms, and B and C, which correspond to different biosamples gram, silhouette and t-SNE analyzes, respectively. The three valid- and the same platform). This effect is remarkably corrected after ation methods show that CuBlock and ComBat separate very clearly normalization. the biological groups A, B, C and D (except for a couple of A points in ComBat’s case). ComBat tends to produce tighter but less cleanly separated clusters for these four groups, as illustrated by both the t- 4 Results and discussion SNE (Supplementary Fig. S9C) and silhouette (Supplementary Fig. The algorithm described in the previous section was applied to the S8C) plots. CuBlock is the only method that clusters the biological data introduced in Section 2.1 after preprocessing (see Section 2.2), groups A and C, and B and D together in the dendrogram plot and the results were compared to those obtained with other normal- (Fig. 3A), and this is also underlined in the corresponding silhouette ization methods as explained in Section 2.3. In line with the objec- plot in Supplementary Figure S10A, showing high and homogeneous tives stated in the introduction, the discussion of the results evolves silhouette values. On the contrary, log2, ComBat, YuGene and to a around the ability of the methods to highlight biological patterns in smaller extent UPC tend to cluster C with D (Fig. 3B–E). In fact, a multi-platform context. We want to note, however, that although log2 and YuGene have difficulties to separate these two groups at CuBlock was not developed for single-platform normalization—this all, while UPC has serious difficulties to separate the groups C and being already well covered by other methods, it can be also used for D from the parent groups (A and B; see also Supplementary Fig. S9). this purpose. Thus, as a side example, Supplementary Figure S7 Supplementary Figures S10 and S11 show silhouette plots using the shows results for a third dataset, available in GEO with accession groups A [ C and B [ D and the platforms as given clusters, respect- ively. We note that even though CuBlock is shown to emphasize the biological differences and Supplementary Figure S11A indicates weak platform clusters, both the t-SNE (Supplementary Fig. S9A) and dendrogram (Fig. 3A) plots show that, within each of the A, B, C, D clusters, the samples are subclustered by platform. As can be seen in these Figures, ComBat mixes the data from the different plat- forms best, while YuGene, UPC and log2 are, approximately in this order, worst at mixing platform data. In Figure 3 and throughout this study, the log2 transform plays the role of control method. As a second potential control, we also used a common approach consisting in platform-specific normaliza- tion of the samples followed by a centering transformation. Supplementary Figure S12 shows the dendrogram plot of the RMA method (background correction, quantile normalization, summar- ization at the probe-set level for Affymetrix platforms and log2 transformation) followed by Z-score transformation (subtracting the sample’s mean and dividing by its standard deviation). It can be observed that the results are only slightly better than those from a log2 transformation, for which reason we kept the latter as the sim- plest approach. 4.2 Reference dataset: three platforms To compare the results from CuBlock and Shambhala [the latter reported by Borisov et al. (2019) for the same dataset], we also per- formed the analysis for the three-platform subset used by the authors of Shambhala, namely AFX, ILM and AG1. They had con- cluded that Shambhala separates well A [ C from B [ D but not A from C or B from D. Using our selection of 500 proteins that best distinguish A from B, when looking at the results for Shambhala in Supplementary Figure S13B, D we observe that, while A [ C forms a relatively clear cluster, all the B [ D points from AG1 samples are clustered with A [ C, making B [ D a well defined cluster only for AFX and ILM. As illustrated by the t-SNE and dendrogram plots and by the negative silhouette values in Supplementary Figure S13F, Shambhala does also not distinguish A, B, C and D from each other well. The results for CuBlock in Supplementary Figure S13A, C, E show the same features already discussed in Section 4.1 using the data for six platforms. 4.3 Experimental dataset Figure 4 and Supplementary Figures S14–S16 show the results Fig. 3. Dendrogram analysis of the reference dataset (six platforms) after normaliza- tion with CuBlock (A), log2 (B), ComBat (C), YuGene (D) and UPC (E). Color bars obtained for the human sperm dataset, after normalization with under the dendrograms indicate the biological group and platform corresponding to CuBlock, log2, ComBat, YuGene, DBNorm and UPC. CuBlock is each leaf; the BP (green) and AU (red) values (see Section 2.3.3) for some selected the only normalization method that significantly distinguishes the clusters are indicated at the origin of the branches two biological groups, T and N. The dendrogram plots in Figure 4 CuBlock 2371 platforms to highlight existing patterns in the samples. In this sense, CuBlock will never lead to a high correlation between replicates (which might be desirable in some studies), but will on the other hand enable comparison across platforms using a selected small set of differentially expressed proteins. 4.4 Reference dataset: missing biological groups CuBlock is a method that works with the actual distribution of the data, without making any assumption on its shape. It is, in that sense, dependent on the actual differences found in the input data at the platform level—which are in turn highlighted by the block treat- ment enabling the uncovering of the different distributions present in the data. To test this dependence, we analyzed again the reference dataset using CuBlock and ComBat while removing the groups B and D from the platforms HUG and AG1. In other words, these two platforms were normalized only with A and C samples. The bio- logical difference between A and C is that 25% of C is made of B RNA samples. The other four platforms were normalized with all four biological groups. As it can be seen in Figure 5A, C, CuBlock results in the clustering of the C samples of HUG and AG1 separate- ly and closer to the D cluster of the other platforms than to their C cluster. However, the A cluster remains a well-defined cluster for all platforms. This suggests that, in the two platforms with missing groups, CuBlock emphasizes the difference between the available data, as predicted above. For HUG and AG1, this means emphasiz- ing the differences between A and C (the only groups it sees), thus bringing C closer to the D cluster formed by the other platforms, since, as C itself, D is also a mixture of A and B. ComBat does how- ever do similarly in this regard (Fig. 5B, D), with the C samples of HUG and AG1 being even more mixed with the D samples, and the A samples of the two platforms getting closer (see t-SNE plot) to the C samples of the four other platforms, some of them ending up clus- tered (see dendrogram) in this group. 5 Conclusion Fig. 4. Dendrogram analysis of the experimental dataset after normalization with CuBlock (A), log2 (B), ComBat (C), YuGene (D), DBNorm (E) and UPC (F). Color We have introduced an algorithm for cross-platform normalization bars under the dendrograms indicate the biological group and platform correspond- of gene-expression microarray data as well as a strategy to validate ing to each leaf; the BP (green) and AU (red) values (see Section 2.3.3) for some cross-platform normalization methods, with a focus on the capacity selected clusters are indicated at the origin of the branches of the algorithm to properly separate samples from different bio- logical groups after normalization and across multiple platforms. and t-SNE plots in Supplementary Figure S15 show that ComBat Overall, CuBlock showed good results on the two datasets used in has troubles to establish a clear boundary between the two groups, this evaluation, a dataset specifically generated for testing and stand- particularly for the ILL2 platform, while the other methods tend to ardization purposes and a dataset from an actual experimental misclassify the ILL2 samples corresponding to the N group. study. CuBlock could always differentiate, clearly, the underlying Similarly to the results for the reference dataset (Section 4.1), biological groups after mixing data from up to 6 different platforms. CuBlock tends to sort the samples by platform within the clusters T Nevertheless, we observed that within each biological group the al- and N (except for one T sample from ILL2). To investigate whether gorithm tends to subcluster samples by platform, indicating a patterns that are found using one platform can be extrapolated to remaining, yet comparatively small, platform effect. The ComBat al- the other platforms, we performed a SVM classification test as gorithm (Johnson et al., 2007) showed also good performance on described in Section 2.3.4. The results are shown in Table 1. In all the reference dataset, with better mixing of data from different plat- cases, CuBlock outperforms the other methods. It is also worth not- forms than all the other methods tested. However, on the experi- ing that no matter which platform is used for the training based on mental dataset, where samples are from different individuals and the CuBlock, the results are always very similar. To a lesser extent, this difference between biological groups might become less obvious is also true for DBNorm. However, using log2, ComBat, YuGene than in the reference dataset, ComBat did not perform as well. and UPC, training with ILL2 gives worse results than with the other Platform mixing was still good but the distinction between the two platforms, probably due to the fact that this platform constitutes a biological groups was not clear. To rationalize the differences in better defined cluster, as shown in the silhouette plots in platform-mixing properties between CuBlock and ComBat, we per- Supplementary Figure S16. formed a rank correlation analysis and show that methods normaliz- In Figure 4, it can be seen that ComBat clusters more strongly ing all platforms together, like ComBat, tend to correlate the data AFF and ILL1 (which contains replicates of AFF) than CuBlock from the different platforms, as opposed to methods normalizing does. In Supplementary Table S1 we provide, for the six normaliza- platforms separately, which tend to decorrelate them. While high tion methods, average pairwise Kendall rank correlations between correlation among samples from different platforms might be desir- the replicates shared by the platforms AFF and ILL1, between the able in some set-ups, there is no guarantee that this correlation will non-replicates of these two platforms and between the other two be meaningful, as shown for ComBat, which in our test induced cor- platform pairs, as well as a discussion on the correlation data. In relation among platforms almost independently of the actual level of summary, we show that methods like ComBat will tend to correlate expected correlation between them. As already mentioned, ComBat samples from different platforms, even when they are not from the also requires the platforms to be normalized together, making it a same biological group (and should be therefore less correlated), less convenient method for systematic application to multiple data- while methods like CuBlock tend to decorrelate the different sets. On the other hand, DBNorm—used only with the experimental 2372 V.Junet et al. Table 1. SVM classification scores (see Section 2.3.4) Training platform Accuracy MCC Balanced accuracy AUC CuBlock AFF 0.88 6 0.08 0.76 6 0.17 0.86 6 0.10 0.98 6 0.03 ILL1 0.92 6 0.06 0.85 6 0.12 0.91 6 0.07 0.99 6 0.02 ILL2 0.86 6 0.16 0.74 6 0.32 0.86 6 0.17 0.99 6 0.04 log2 AFF 0.77 6 0.07 0.52 6 0.16 0.71 6 0.08 0.70 6 0.11 ILL1 0.80 6 0.06 0.64 6 0.11 0.81 6 0.06 0.81 6 0.04 ILL2 0.44 6 0.16 –0.13 6 0.35 0.46 6 0.16 0.24 6 0.31 ComBat AFF 0.73 6 0.11 0.53 6 0.20 0.76 6 0.10 0.90 6 0.08 ILL1 0.76 6 0.08 0.57 6 0.15 0.77 6 0.08 0.90 6 0.06 ILL2 0.61 6 0.30 0.22 6 0.61 0.61 6 0.30 0.63 6 0.36 YuGene AFF 0.81 6 0.07 0.63 6 0.15 0.77 6 0.08 0.91 6 0.06 ILL1 0.87 6 0.08 0.74 6 0.16 0.87 6 0.09 0.96 6 0.04 ILL2 0.55 6 0.07 0.07 6 0.17 0.53 6 0.07 0.97 6 0.11 DBNorm AFF 0.81 6 0.07 0.62 6 0.16 0.77 6 0.09 0.91 6 0.07 ILL1 0.86 6 0.06 0.76 6 0.09 0.87 6 0.05 0.95 6 0.04 ILL2 0.79 6 0.13 0.61 6 0.24 0.77 6 0.13 0.96 6 0.11 UPC AFF 0.74 6 0.08 0.45 6 0.22 0.68 6 0.10 0.82 6 0.08 ILL1 0.82 6 0.06 0.68 6 0.09 0.83 6 0.05 0.90 6 0.05 ILL2 0.57 6 0.10 0.10 6 0.24 0.54 6 0.11 0.75 6 0.34 Note: Mean and standard deviation over the 6000 models per platform and method. Finally, we also showed, by training SVMs with single-platform data from the experimental set, that when normalizing the data with CuBlock the patterns that are found using one platform can be extrapolated to the other platforms significantly better than when the normalization is done with any other of the methods. Because no assumptions are made on the distributions underly- ing the input data, CuBlock can be thought of as a transformation which result remains close to the input. CuBlock fits cubic polyno- mials to data blocks that are found by k-means clustering, thereby trying to best fit the different distributions found in the data corre- sponding to a sample (different blocks need not have the same distri- bution) and it does so without assuming a shape for these distributions. As a consequence, CuBlock emphasizes the differences within the input data at the platform level, making it sample- composition dependent despite the only step in the algorithm where the microarray samples from a given platform are considered to- gether is when applying the k-means algorithm (afterwards, each sample is considered separately). This sample dependence was high- lighted in this study when analyzing the reference dataset after removing biological groups in some platforms. Nevertheless, the sample dependence is even more prominent for algorithms normaliz- ing platforms together, such as ComBat. We note that, while in this study the results were presented at the protein level, CuBlock is applied and returns the normalized data at the probe level, as an appropriate level for gene-expression microar- Fig. 5. t-SNE dimension reduction and dendrogram plots for the reference dataset, with exclusion of the B and D samples from platforms HUG and AG1, after normal- rays. In fact, like most other methods, it could technically be applied ization with CuBlock and ComBat. Point color and shape indicate biological group at any expression  samples level. Studying the effect of the applica- and platform, respectively (right-hand legend). (A) t-SNE for CuBlock normalized tion of CuBlock and other cross-platform normalization methods at data; perplexity (Prp) and mean silhouette index (SI) values (see Section 2.3.2): Prp different levels with appropriate (and possibly study-specific) pre- ¼ 15, SI ¼ 0.80. (B) Corresponding analysis for ComBat-normalized data; Prp ¼ processing steps could be an interesting development which would 10, SI ¼ 0.67. (C, D) Dendrograms for CuBlock (C) and ComBat (D) normalized ideally lead to an appropriate systematic integration of RNA-seq data; color bars under the dendrograms indicate the biological group and platform with gene expression microarrays. This however goes beyond the corresponding to each leaf scope of this study. In summary, we have shown that CuBlock can be applied to data from multiple microarrays in a platform agnostic way and pre- dataset, as it proved computationally much more time demanding serves the biological grouping of the samples, demonstrating a good than the rest, YuGene and UPC performed only slightly better than performance for different types of samples. It is therefore a tool ap- log2. For the set evaluated, Shambhala lagged clearly behind the propriate for gene-expression studies based on multiple microarray other methods, arguably including a simple log2 transformation. sets collected along multiple platforms and at different times, thus CuBlock 2373 facilitating the extraction of knowledge from the wealth of micro- Maaten,L.v.d. and Hinton,G. (2008) Visualizing data using t-SNE. J. Mach. array data available in public repositories and enabling the use of Learn. Res., 9, 2579–2605. these repositories as sources of Real-World Data. Maire,V. et al. (2013a) Polo-like Kinase 1: a potential therapeutic option in combination with conventional chemotherapy for the management of patients with triple-negative breast cancer. Cancer Res., 73, 813–823. Maire,V. et al. (2013b) TTK/hMPS1 is an attractive therapeutic target for Acknowledgements triple-negative breast cancer. PLoS One, 8, e63712–15. The authors thank Malu Calle for useful discussions during the preparation MAQC Consortium (2006) The MicroArray Quality Control (MAQC) pro- of the manuscript. ject shows inter- and intraplatform reproducibility of gene expression meas- urements. Nat. Biotechnol., 24, 1151–1161. Maubant,S. et al. (2015) Transcriptome analysis of Wnt3a-treated triple-nega- Funding tive breast cancer cells. PLoS One, 10, e0122333–26. Meng,Q. et al. (2017) DBNorm: normalizing high-density oligonucleotide This project received funding from the European Union’s Horizon 2020 re- microarray data based on distributions. BMC Bioinformatics, 18, 527. search and innovation programme under the Marie Skłodowska-Curie grant Nagalakshmi,U. et al. (2008) The transcriptional landscape of the yeast gen- agreement No. 765158 (COSMIC; www.cosmic-h2020.eu). ome defined by RNA sequencing. Science, 320, 1344–1349. Piccolo,S.R. et al. (2012) A single-sample microarray normalization method to Conflict of Interest: none declared. facilitate personalized-medicine workflows. Genomics, 100, 337–344. Piccolo,S.R. et al. (2013) Multiplatform single-sample estimates of transcrip- tional activation. Proc. Natl. Acad. Sci. USA, 110, 17778–17783. References Platts,A.E. et al. (2007) Success and failure in human spermatogenesis as Bartlett,V.L. et al. (2019) Feasibility of using real-world data to replicate clin- revealed by teratozoospermic RNAs. Hum. Mol. Genet., 16, 763–773. ical trial evidence. JAMA Netw. Open, 2, e1912869. Rapaport,F. et al. (2013) Comprehensive evaluation of differential gene ex- Benito,M. et al. (2004) Adjustment of systematic microarray data biases. pression analysis methods for RNA-seq data. Genome Biol., 14, R95. Bioinformatics, 20, 105–114. Rousseeuw,P.J. (1987) Silhouettes: a graphical aid to the interpretation and Berger,M.L. et al. (2017) Good practices for real-world data studies of treat- validation of cluster analysis. J. Comput. Appl. Math., 20, 53–65. ment and/or comparative effectiveness: recommendations from the joint Rudy,J. and Valafar,F. (2011) Empirical comparison of cross-platform nor- ISPOR-ISPE Special Task Force on real-world evidence in health care deci- malization methods for gene expression data. BMC Bioinformatics, 12, sion making. Pharmacoepidemiol. Drug Saf., 26, 1033–1039. 467. Borisov,N. et al. (2019) Shambhala: a platform-agnostic data harmonizer for Schena,M. et al. (1995) Quantitative monitoring of gene expression patterns gene expression data. BMC Bioinformatics, 20, 66. with a complementary DNA microarray. Science, 270, 467–470. Boughorbel,S. et al. (2017) Optimal classifier for imbalanced data using SEQC/MAQC-III Consortium. (2014) A comprehensive assessment of Matthews Correlation Coefficient metric. PLoS One, 12, e0177678. RNA-seq accuracy, reproducibility and information content by the Bumgarner,R. (2013) Overview of DNA microarrays: types, applications, and Sequencing Quality Control Consortium. Nat. Biotechnol., 32, 903–914. their future. Curr. Protoc. Mol. Biol., 101, 22.1.1–22.1.11. Shabalin,A.A. et al. (2008) Merging two gene-expression studies via Cortes,C. and Vapnik,V. (1995) Support-vector networks. Mach. Learn., 20, cross-platform normalization. Bioinformatics, 24, 1154–1160. 273–297. Sherman,R.E. et al. (2017) Accelerating development of scientific evidence for Fawcett,T. (2006) An introduction to ROC analysis. Pattern Recogn. Lett., medical products within the existing US regulatory framework. Nat. Rev. 27, 861–874. Drug Discov., 16, 297–298. U.S. Food and Drug Administration. (2018) Framework for FDA’s Real-World Shimodaira,H. (2004) Approximately unbiased tests of regions using Evidence Program. https://www.fda.gov/media/120060/download multistep-multiscale bootstrap resampling. Ann. Stat., 32, 2616–2641. Hubbell,E. et al. (2002) Robust estimators for expression analysis. Storey,J.D. (2002) A direct approach to false discovery rates. J. R. Stat. Soc. B, Bioinformatics, 18, 1585–1592. 64, 479–498. Irigoyen,A. et al. (2018) Integrative multi-platform meta-analysis of gene ex- Suzuki,R. and Shimodaira,H. (2006) Pvclust: an R package for assessing the pression profiles in pancreatic ductal adenocarcinoma patients for identify- uncertainty in hierarchical clustering. Bioinformatics, 22, 1540–1542. ing novel diagnostic biomarkers. PLoS One, 13, e0194844. Trotta,F. (2012) Discrepancies between observational studies and randomized Irizarry,R.A. et al. (2003) Exploration, normalization, and summaries of high controlled trials. Focus Farmacovigilanza, 73, 1. density oligonucleotide array probe level data. Biostatistics, 4, 249–264. Walsh,C.J. et al. (2015) Microarray meta-analysis and cross-platform normal- Johnson,W.E. et al. (2007) Adjusting batch effects in microarray expression ization: integrative genomics for robust biomarker discovery. Microarrays, data using empirical Bayes methods. Biostatistics, 8, 118–127. 4, 389–406. Lashkari,D.A. et al. (1997) Yeast microarrays for genome wide parallel genetic Yang,Y.H. et al. (2002) Normalization for cDNA microarray data: a robust and gene expression analysis. Proc. Natl. Acad. Sci. USA, 94, 13057–13062. composite method addressing single and multiple slide systematic variation. Lê Cao,K.-A. et al. (2014) YuGene: a simple approach to scale gene expression Nucleic Acids Res., 30, e15. data derived from different platforms for integrated analyses. Genomics, Zhang,S. et al. (2020) MatchMixeR: a cross-platform normalization method 103, 239–251. for gene expression data integration. Bioinformatics, 36, 2486–2491. Lloyd,S. (1982) Least squares quantization in PCM. IEEE Trans. Inf. Theory, 28, 129–137.