MINT: a multivariate integrative method to identify reproducible molecular signatures across independent experiments and platforms

Abstract

Background

Molecular signatures identified from high-throughput transcriptomic studies often have poor reliability and fail to reproduce across studies. One solution is to combine independent studies into a single integrative analysis, additionally increasing sample size. However, the different protocols and technological platforms across transcriptomic studies produce unwanted systematic variation that strongly confounds the integrative analysis results. When studies aim to discriminate an outcome of interest, the common approach is a sequential two-step procedure; unwanted systematic variation removal techniques are applied prior to classification methods.

Results

To limit the risk of overfitting and over-optimistic results of a two-step procedure, we developed a novel multivariate integration method, MINT, that simultaneously accounts for unwanted systematic variation and identifies predictive gene signatures with greater reproducibility and accuracy. In two biological examples on the classification of three human cell types and four subtypes of breast cancer, we combined high-dimensional microarray and RNA-seq data sets and MINT identified highly reproducible and relevant gene signatures predictive of a given phenotype. MINT led to superior classification and prediction accuracy compared to the existing sequential two-step procedures.

Conclusions

MINT is a powerful approach and the first of its kind to solve the integrative classification framework in a single step by combining multiple independent studies. MINT is computationally fast as part of the mixOmics R CRAN package, available at http://www.mixOmics.org/mixMINT/and http://cran.r-project.org/web/packages/mixOmics/.

Background

High-throughput technologies, based on microarray and RNA-sequencing, are now being used to identify biomarkers or gene signatures that distinguish disease subgroups, predict cell phenotypes or classify responses to therapeutic drugs. However, few of these findings are reproduced when assessed in subsequent studies and even fewer lead to clinical applications [1, 2]. The poor reproducibility of identified gene signatures is most likely a consequence of high-dimensional data, in which the number of genes or transcripts being analysed is very high (often several thousands) relative to a comparatively small sample size being used (<20).

One way to increase sample size is to combine raw data from independent experiments in an integrative analysis. This would improve both the statistical power of the analysis and the reproducibility of the gene signatures that are identified [3]. However, integrating transcriptomic studies with the aim of classifying biological samples based on an outcome of interest (integrative classification) has a number of challenges. Transcriptomic studies often differ from each other in a number of ways, such as in their experimental protocols or in the technological platform used. These differences can lead to so-called ‘batch-effects’, or systematic variation across studies, which is an important source of confounding [4]. Technological platform, in particular, has been shown to be an important confounder that affects the reproducibility of transcriptomic studies [5]. In the MicroArray Quality Control (MAQC) project, poor overlap of differentially expressed genes was observed across different microarray platforms (∼ 60%), with low concordance observed between microarray and RNA-seq technologies specifically [6]. Therefore, these confounding factors and sources of systematic variation must be accounted for, when combining independent studies, to enable genuine biological variation to be identified.

The common approach to integrative classification is sequential. A first step consists of removing batch-effect by applying for instance ComBat [7], FAbatch [8], Batch Mean-Centering [9], LMM-EH-PS [10], RUV-2 [4] or YuGene [11]. A second step fits a statistical model to classify biological samples and predict the class membership of new samples. A range of classification methods also exists for these purposes, including machine learning approaches (e.g. random forests [12, 13] or Support Vector Machine [14–16]) as well as multivariate linear approaches (Linear Discriminant Analysis LDA, Partial Least Square Discriminant Analysis PLSDA [17], or sparse PLSDA [18]).

The major pitfall of the sequential approach is a risk of over-optimistic results from overfitting of the training set. This leads to signatures that cannot be reproduced on test sets. Moreover, most proposed classification models have not been objectively validated on an external and independent test set. Thus, spurious conclusions can be generated when using these methods, leading to limited potential for translating results into reliable clinical tools [2]. For instance, most classification methods require the choice of a parameter (e.g. sparsity), which is usually optimised with cross-validation (data are divided into k subsets or ‘folds’ and each fold is used once as an internal test set). Unless the removal of batch-effects is performed independently on each fold, the folds are not independent and this leads to over-optimistic classification accuracy on the internal test sets. Hence, batch removal methods must be used with caution. For instance, ComBat can not remove unwanted variation in an independent test set alone as it requires the test set to be normalised with the learning set in a transductive rather than inductive approach [19]. This is a clear example where over-fitting and over-optimistic results can be an issue, even when a test set is considered.

To address existing limitations of current data integration approaches and the poor reproducibility of results, we propose a novel Multivariate INTegrative method, MINT. MINT is the first approach of its kind that integrates independent data sets while simultaneously, accounting for unwanted (study) variation, classifying samples and identifying key discriminant variables. MINT predicts the class of new samples from external studies, which enables a direct assessment of its performance. It also provides insightful graphical outputs to improve interpretation and inspect each study during the integration process.

We validated MINT in a subset of the MAQC project, which was carefully designed to enable assessment of unwanted systematic variation. We then combined microarray and RNA-seq experiments to classify samples from three human cell types (human Fibroblasts (Fib), human Embryonic Stem Cells (hESC) and human induced Pluripotent Stem Cells (hiPSC)) and from four classes of breast cancer (subtype Basal, HER2, Luminal A and Luminal B). We use these datasets to demonstrate the reproducibility of gene signatures identified by MINT.

Methods

We use the following notations. Let X denote a data matrix of size N observations (rows) ×P variables (e.g. gene expression levels, in columns) and Y a dummy matrix indicating each sample class membership of size N observations (rows) ×K categories outcome (columns). We assume that the data are partitioned into M groups corresponding to each independent study m: {(X ⁽¹⁾,Y ⁽¹⁾),…,(X ^(M),Y ^(M))} so that $\sum _{m=1}^{M} n_{m}=N$, where n _m is the number of samples in group m, see Additional file 1: Figure S1. Each variable from the data set X ^(m) and Y ^(m) is centered and has unit variance. We write X and Y the concatenation of all X ^(m) and Y ^(m), respectively. Note that if an internal known batch effect is present in a study, this study should be split according to that batch effect factor into several sub-studies considered as independent. For $n\in \mathbb {N}$, we denote for all $a\in \mathbb {R}^{n}$ its ℓ ¹ norm $||a||_{1}=\sum _{1}^{n}|a_{j}|$ and its ℓ ² norm $||a||_{2}=\left (\sum _{1}^{n}a_{j}^{2}\right)^{1/2}$ and |a|₊ the positive part of a. For any matrix we denote by ^⊤ its transpose.

PLS-based classification methods to combine independent studies

PLS approaches have been extended to classify samples Y from a data matrix X by maximising a formula based on their covariance. Specifically, latent components are built based on the original X variables to summarise the information and reduce the dimension of the data while discriminating the Y outcome. Samples are then projected into a smaller space spanned by the latent component. We first detail the classical PLS-DA approach and then describe mgPLS, a PLS-based model we previously developed to model a group (study) structure in X.

PLS-DA Partial Least Squares Discriminant Analysis [17] is an extension of PLS for a classification frameworks where Y is a dummy matrix indicating sample class membership. In our study, we applied PLS-DA as an integrative approach by naively concatenating all studies. Briefly, PLS-DA is an iterative method that constructs H successive artificial (latent) components t _h =X _h a _h and u _h =Y _h b _h for h=1,..,H, where the h ^th component t _h (respectively u _h ) is a linear combination of the X (Y) variables. H denotes the dimension of the PLS-DA model. The weight coefficient vector a _h (b _h ) is the loading vector that indicates the importance of each variable to define the component. For each dimension h=1,…,H PLS-DA seeks to maximize

$$ \underset{||a_{h}||_{2} = ||b_{h}||_{2} =1}{\max }cov(X_{h} a_{h}, Y_{h} b_{h}), $$

(1)

where X _h ,Y _h are residual matrices (obtained through a deflation step, as detailed in [18]). The PLS-DA algorithm is described in Additional file 1: Supplemental Material S1. The PLS-DA model assigns to each sample i a pair of H scores $(t_{h}^{i}, u_{h}^{i})$ which effectively represents the projection of that sample into the X- or Y- space spanned by those PLS components. As H<<P, the projection space is small, allowing for dimension reduction as well as insightful sample plot representation (e.g. graphical outputs in “Results” section). While PLS-DA ignores the data group structure inherent to each independent study, it can give satisfactory results when the between groups variance is smaller than the within group variance or when combined with extensive data subsampling to account for systematic variation across platforms [21].

mgPLS Multi-group PLS is an extension of the PLS framework we recently proposed to model grouped data [22, 23], which is relevant for our particular case where the groups represent independent studies. In mgPLS, the PLS-components of each group are constraint to be built based on the same loading vectors in X and Y. These global loading vectors thus allow the samples from each group or study to be projected in the same common space spanned by the PLS-components. We extended the original unsupervised approach to a supervised approach by using a dummy matrix Y as in PLS-DA to classify samples while modelling the group structure. For each dimension h=1,…,H mgPLS-DA seeks to maximize

$$ \underset{||a_{h}||_{2} = ||b_{h}||_{2} =1}{\max }\sum_{m=1}^{M} n_{m} cov\left(X^{(m)}_{h}a_{h}, Y^{(m)}_{h} b_{h}\right), $$

(2)

where a _h and b _h are the global loadings vectors common to all groups, $t_{h}^{(m)}=X_{h}^{(m)}a_{h}$ and $u_{h}^{(m)}=Y_{h}^{(m)}b_{h}$ are the group-specific (partial) PLS-components, and $X_{h}^{(m)}$ and $Y_{h}^{(m)}$ are the residual (deflated) matrices. The global loadings vectors (a _h ,b _h ) and global components (t _h =X _h a _h ,u _h =Y _h b _h ) enable to assess overall classification accuracy, while the group-specific loadings and components provide powerful graphical outputs for each study that is integrated in the analysis. Global and group-specific components and loadings are represented in Additional file 1: Figure S2. The next development we describe below is to include internal variable selection in mgPLS-DA for large dimensional data sets.

MINT

Our novel multivariate integrative method MINT simultaneously integrates independent studies and selects the most discriminant variables to classify samples and predict the class of new samples. MINT seeks for a common projection space for all studies that is defined on a small subset of discriminative variables and that display an analogous discrimination of the samples across studies. The identified variables share common information across all studies and therefore represent a reproducible signature that helps characterising biological systems. MINT further extends mgPLS-DA by including a ℓ ¹-penalisation on the global loading vector a _h to perform variable selection. For each dimension h=1,…,H the MINT algorithm seeks to maximize

$$ \underset{||a_{h}||_{2} = ||b_{s}||_{2} =1}{\max }\sum_{m=1}^{M} n_{m} cov(X_{h}^{(m)}a_{h}, Y_{h}^{(m)}b_{h}) + \lambda_{h}||a_{h}||_{1}, $$

(3)

where in addition to the notations from Eq. (2), λ _h is a non negative parameter that controls the amount of shrinkage on the global loading vectors a _h and thus the number of non zero weights. Similarly to Lasso [24] or sparse PLS-DA [18], the added ℓ ¹ penalisation in MINT improves interpretability of the PLS-components that are now defined only on a set of selected biomarkers from X (with non zero weight) that are identified in the linear combination $X_{h}^{(m)}a_{h}$. The ℓ ¹ penalisation in effectively solved in the MINT algorithm using soft-thresholding (see pseudo Algorithm 1).

In addition to the integrative classification framework, MINT was extended to an integrative regression framework (multiple multivariate regression, Additional file 1 Supplemental Material S2).

Class prediction and parameters tuning with MINT

MINT centers and scales each study from the training set, so that each variable has mean 0 and variance 1, similarly to any PLS methods. Therefore, a similar pre-processing needs to be applied on test sets. If a test sample belongs to a study that is part of the training set, then we apply the same scaling coefficients as from the training study. This is required so that MINT applied on a single study will provide the same results as PLS. If the test study is completely independent, then it is centered and scaled separately.

After scaling the test samples, the prediction framework of PLS is used to estimate the dummy matrix Y _test of an independent test set X _test [25], where each row in Y _test sums to 1, and each column represents a class of the outcome. A class membership is assigned (predicted) to each test sample by using the maximal distance, as described in [18]. It consists in assigning the class with maximal positive value in Y _test .

The main parameter to tune in MINT is the penalty λ _h for each PLS-component h, which is usually performed using Cross-Validation (CV). In practice, the parameter λ _h can be equally replaced by the number of variables to select on each component, which is our preferred user-friendly option. The assessment criterion in the CV can be based on the proportion of misclassified samples, proportion of false or true positives, or, as in our case, the balanced error rate (BER). BER is calculated as the averaged proportion of wrongly classified samples in each class and weights up small sample size classes. We consider BER to be a more objective performance measure than the overall misclassification error rate when dealing with unbalanced classes. MINT tuning is computationally efficient as it takes advantage of the group data structure in the integrative study. We used a “Leave-One-Group-Out Cross-Validation (LOGOCV)”, which consists in performing CV where group or study m is left out only once m=1,…,M. LOGOCV realistically reflects the true case scenario where prediction is performed on independent external studies based on a reproducible signature identified on the training set. Finally, the total number of components H in MINT is set to K−1, K= number of classes, similar to PLS-DA and ℓ ¹ penalised PLS-DA models [18].

Case studies

We demonstrate the ability of MINT to identify the true positive genes on the MAQC project, then highlight the strong properties of our method to combine independent data sets in order to identify reproducible and predictive gene signatures on two other biological studies.

The MicroArray quality control (MAQC) project. The extensive MAQC project focused on assessing microarray technologies reproducibility in a controlled environment [5]. Two reference samples, RNA samples Universal Human Reference (UHR) and Human Brain Reference (HBR) and two mixtures of the original samples were considered. Technical replicates were obtained from three different array platforms -Illumina, AffyHuGene and AffyPrime- for each of the four biological samples A (100% UHR), B (100% HBR), C (75% UHR, 25% HBR) and D (25% UHR and 75% HBR). Data were downloaded from Gene Expression Omnibus (GEO) - GSE56457. In this study, we focused on identifying biomarkers that discriminate A vs. B and C vs. D. The experimental design is referenced in Additional file 1: Table S1.

Stem cells. We integrated 15 transcriptomics microarray datasets to classify three types of human cells: human Fibroblasts (Fib), human Embryonic Stem Cells (hESC) and human induced Pluripotent Stem Cells (hiPSC). As there exists a biological hierarchy among these three cell types, two sub-classification problems are of interest in our analysis, which we will address simultaneously with MINT. On the one hand, differences between pluripotent (hiPSC and hESC) and non-pluripotent cells (Fib) are well-characterised and are expected to contribute to the main biological variation. Our first level of analysis will therefore benchmark MINT against the gold standard in the field. On the other hand, hiPSC are genetically reprogrammed to behave like hESC and both cell types are commonly assumed to be alike. However, differences have been reported in the literature [26–28], justifying the second and more challenging level of classification analysis between hiPSC and hESC. We used the cell type annotations of the 342 samples as provided by the authors of the 15 studies.

The stem cell dataset provides an excellent showcase study to benchmark MINT against existing statistical methods to solve a rather ambitious classification problem.

Each of the 15 studies was assigned to either a training or test set. Platforms uniquely represented were assigned to the training set and studies with only one sample in one class were assigned to the test set. Remaining studies were randomly assigned to training or test set. Eventually, the training set included eight datasets (210 samples) derived on five commercial platforms and the independent test set included the remaining seven datasets (132 samples) derived on three platforms (Table 1).

Table 1 Stem cells experimental design

The pre-processed files were downloaded from the http://www.stemformatics.org collaborative platform [29]. Each dataset was background corrected, log2 transformed, YuGene normalized and mapped from probes ID to Ensembl ID as previously described in [11], resulting in 13 313 unique Ensembl gene identifiers. In the case where datasets contained multiple probes for the same Ensembl ID gene, the highest expressed probe was chosen as the representative of that gene in that dataset. The choice of YuGene normalisation was motivated by the need to normalise each sample independently rather than as a part of a whole study (e.g. existing methods ComBat [7], quantile normalisation (RMA [30])), to effectively limit over-fitting during the CV evaluation process.

Breast cancer. We combined whole-genome gene-expression data from two cohorts from the Molecular Taxonomy of Breast Cancer International Consortium project (METABRIC, [31] and of two cohorts from the Cancer Genome Atlas (TCGA, [32]) to classify the intrinsic subtypes Basal, HER2, Luminal A and Luminal B, as defined by the PAM50 signature [20]. The METABRIC cohorts data were made available upon request, and were processed by [31]. TCGA cohorts are gene-expression data from RNA-seq and microarray platforms. RNA-seq data were normalised using Expectation Maximisation (RSEM) and percentile-ranked gene-level transcription estimates. The microarray data were processed as described in [32].

The training set consisted in three cohorts (TCGA RNA-seq and both METABRIC microarray studies), including the expression levels of 15 803 genes on 2 814 samples; the test set included the TCGA microarray cohort with 254 samples (Table 2). Two analyses were conducted, which either included or discarded the PAM50 genes from the data. The first analysis aimed at recovering the PAM50 genes used to classify the samples. The second analysis was performed on 15,755 genes and aimed at identifying an alternative signature to the PAM50.

Table 2 Experimental design of four breast cancer cohorts including 4 cancer subtypes: Basal, HER2, Luminal A (LumA) and Luminal B (LumB)

Performance comparison with sequential classification approaches

We compared MINT with sequential approaches that combine batch-effect removal approaches with classification methods. As a reference, classification methods were also used on their own on a naive concatenation of all studies. Batch-effect removal methods included Batch Mean-Centering (BMC, [9]), ComBat [7], linear models (LM) or linear mixed models (LMM), and classification methods included PLS-DA, sPLS-DA [18], mgPLS [22, 23] and Random forests (RF [12]). For LM and LMM, linear models were fitted on each gene and the residuals were extracted as a batch-corrected gene expression [33, 34]. The study effect was set as a fixed effect with LM or as a random effect with LMM. No sample outcome (e.g. cell-type) was included.

Prediction with ComBat normalised data were obtained as described in [19]. In this study, we did not include methods that require extra information -as control genes with RUV-2 [4]- and methods that are not widely available to the community as LMM-EH [10]. Classification methods were chosen so as to simultaneously discriminate all classes. With the exception of sPLS-DA, none of those methods perform internal variable selection. The multivariate methods PLS-DA, mgPLS and sPLS-DA were run on K−1 components, sPLS-DA was tuned using 5-fold CV on each component. All classification methods were combined with batch-removal method with the exception of mgPLS that already includes a study structure in the model.

MINT and PLS-DA-like approaches use a prediction threshold based on distances (see “Class prediction and parameters tuning with MINT ” section) that optimally determines class membership of test samples, and as such do not require receiver operating characteristic (ROC) curves and area under the curve (AUC) performance measures. In addition, those measures are limited to binary classification which do not apply for our stem cell and breast cancer multi-class studies. Instead we use Balanced classification Error Rate to objectively evaluate the classification and prediction performance of the methods for unbalanced sample size classes (“ MINT ” section). Classification accuracies for each class were also reported.

Results

Validation of the MINT approach to identify signatures agnostic to batch effect

The MAQC project processed technical replicates of four well-characterised biological samples A, B, C and D across three platforms. Thus, we assumed that genes that are differentially expressed (DEG) in every single platform are true positive. We primarily focused on identifying biomarkers that discriminate C vs. D, and report the results of A vs. B in the Additional file 1: Supplemental Material S3, Figure S3. Differential expression analysis of C vs. D was conducted on each of the three microarray platforms using ANOVA, showing an overlap of 1385 DEG (FDR <10⁻³ [35]), which we considered as true positive. This corresponded to 62.6% of all DEG for Illumina, 30.5% for AffyHuGene and 21.0% for AffyPrime (Additional file 1: Figure S4). We observed that conducting a differential analysis on the concatenated data from the three microarray platforms without accommodating for batch effects resulted in 691 DEG, of which only 56% (387) were true positive genes. This implies that the remaining 44% (304) of these genes were false positive, and hence were not DE in at least one study. The high percentage of false positive was explained by a Principal Component Analysis (PCA) sample plot that showed samples clustering by platforms (Additional file 1: Figure S4), which confirmed that the major source of variation in the combined data was attributed to platforms rather than cell types.

MINT selected a single gene, BCAS1, to discriminate the two biological classes C and D. BCAS1 was a true positive gene, as part of the common DEG, and was ranked 1 for Illumina, 158 for AffyPrime and 1182 for AffyHuGene. Since the biological samples C and D are very different, the selection of one single gene by MINT was not surprising. To further investigate the performance of MINT, we expanded the number of genes selected by MINT, by decreasing its sparsity parameter (see Methods), and compared the overlap between this larger MINT signature and the true positive genes. We observed an overlap of 100% for a MINT signature of size 100, and an overlap of 89% for a signature of size 1385, which is the number of common DEG identified previously. The high percentage of true positive selected by MINT demonstrates its ability to identify a signature agnostic to batch effect.

Limitations of common meta-analysis and integrative approaches

A meta-analysis of eight stem cell studies, each including three cell types (Table 1, stem cell training set), highlighted a small overlap of DEG lists obtained from the analysis of each separate study (FDR <10⁻⁵, ANOVA, Additional file 1: Table S2). Indeed, the Takahashi study with only 24 DEG limited the overlap between all eight studies to only 5 DEG. This represents a major limitation of merging pre-analysed gene lists as the concordance between DEG lists decreases when the number of studies increases.

One alternative to meta-analysis is to perform an integrative analysis by concatenating all eight studies. Similarly to the MAQC analysis, we first observed that the major source of variation in the combined data was attributed to study rather than cell type (Fig. 1 a). PLS-DA was applied to discriminate the samples according to their cell types, and it showed a strong study variation (Fig. 1 b), despite being a supervised analysis. Compared to unsupervised PCA (Fig. 1 a), the study effect was reduced for the fibroblast cells, but was still present for the similar cell types hESC and hiPSC. We reached similar conclusions when analysing the breast cancer data (Additional file 1: Supplemental Material S4, Figure S5).

Fig. 1

Stem cell study. a PCA on the concatenated data: a greater study variation than a cell type variation is observed. b PLSDA on the concatenated data clustered Fibroblasts only. c MINT sample plot shows that each cell type is well clustered, d MINT performance: BER and classification accuracy for each cell type and each study

MINT outperforms state-of-the-art methods

We compared the classification accuracy of MINT to sequential methods where batch removal methods were applied prior to classification methods. In both stem cell and breast cancer studies, MINT led to the best accuracy on the training set and the best reproducibility of the classification model on the test set (lowest Balanced Error Rate, BER, Fig. 2, Additional file 1: Figures S6 and S7). In addition, MINT consistently ranked first as the best performing method, followed by ComBat+sPLSDA with an average rank of 4.5 (Additional file 1: Figure S8).

Fig. 2

Classification accuracy for both training and test set for the stem cells and breast cancer studies (excluding PAM50 genes). The classification Balanced Error Rates (BER) are reported for all sixteen methods compared with MINT (in black)

On the stem cell data, we found that fibroblasts were the easiest to classify for all methods, including those that do not accommodate unwanted variation (PLS-DA, sPLS-DA and RF, Additional file 1: Figure S6). Classifying hiPSC vs. hESC proved more challenging for all methods, leading to a substantially lower classification accuracy than fibroblasts.

The analysis of the breast cancer data (excluding PAM50 genes) showed that methods that do not accommodate unwanted variation were able to rightly classify most of the samples from the training set, but failed at classifying any of the four subtypes on the external test set. As a consequence, all samples were predicted as LumB with PLS-DA and sPLS-DA, or Basal with RF (Additional file 1: Figure S7). Thus, RF gave a satisfactory performance on the training set (BER =18.5), but a poor performance on the test set (BER =75).

Additionally, we observed that the biomarker selection process substantially improved classification accuracy. On the stem cell data, LM+sPLSDA and MINT outperformed their non sparse counterparts LM+PLSDA and mgPLS (Fig. 2, BER of 9.8 and 7.1 vs. 20.8 and 11.9), respectively.

Finally, MINT was largely superior in terms of computational efficiency. The training step on the stem cell data which includes 210 samples and 13,313 was run in 1 s, compared to 8 s with the second best performing method ComBat+sPLS-DA (2013 MacNook Pro 2.6 Ghz, 16 Gb memory). The popular method ComBat took 7.1s to run, and sPLS-DA 0.9s. The training step on the breast cancer data that includes 2817 samples and 15,755 genes was run in 37s for MINT and 71.5s for ComBat(30.8s)+sPLS-DA(40.6s).

Study-specific outputs with MINT

One of the main challenges when combining independent studies is to assess the concordance between studies. During the integration procedure, MINT proposes not only individual performance accuracy assessment, but also insightful graphical outputs that are study-specific and can serve as Quality Control step to detect outlier studies. One particular example is the Takahashi study from the stem cell data, whose poor performance (Fig. 1 d) was further confirmed on the study-specific outputs (Additional file 1: Figure S9). Of note, this study was the only one generated through Agilent technology and its sample size only accounted for 4.2% of the training set.

The sample plots from each individual breast cancer data set showed the strong ability of MINT to discriminate the breast cancer subtypes while integrating data sets generated from disparate transcriptomics platforms, microarrays and RNA-sequencing (Fig. 3 a–c). Those data sets were all differently pre-processed, and yet MINT was able to model an overall agreement between all studies; MINT successfully built a space based on a handful of genes in which samples from each study are discriminated in a homogenous manner.

Fig. 3

MINT study-specific sample plots showing the projection of samples from a METABRIC Discovery, b METABRIC Validation and c TCGA-RNA-seq experiments, in the same subspace spanned by the first two MINT components. The same subspace is also used to plot the (d) overall (integrated) data. e Balanced Error Rate and classification accuracy for each study and breast cancer subtype from the MINT analysis

MINT gene signature identified promising biomarkers

MINT is a multivariate approach that builds successive components to discriminate all categories (classes) indicated in an outcome variable. On the stem cell data, MINT selected 2 and 15 genes on the first two components respectively (Additional file 1: Table S3). The first component clearly segregated the pluripotent cells (fibroblasts) vs. the two non-pluripotent cell types (hiPSC and hESC) (Fig. 1 c, d). Those non pluripotent cells were subsequently separated on component two with some expected overlap given the similarities between hiPSC and hESC. The two genes selected by MINT on component 1 were LIN28A and CAR which were both found relevant in the literature. Indeed, LIN28A was shown to be highly expressed in ESCs compared to Fibroblasts [36, 37] and CAR has been associated to pluripotency [38]. Finally, despite the high heterogeneity of hiPSC cells included in this study, MINT gave a high accuracy for hESC and hiPSC on independent test sets (93.9% and 77.9% respectively, Additional file 1: Figure S6), suggesting that the 15 genes selected by MINT on component 2 have a high potential to explain the differences between those cell types (Additional file 1: Table S3).

On the breast cancer study, we performed two analyses which either included or discarded the PAM50 genes that were used to define the four cancer subtypes Basal, HER2, Luminal A and Luminal B [20]. In the first analysis, we aimed to assess the ability of MINT to specifically identify the PAM50 key driver genes. MINT successfully recovered 37 of the 48 PAM50 genes present in the data (77%) on the first three components (7, 20 and 10 respectively). The overall signature included 30, 572 and 636 genes on each component (see Additional file 1: Table S4), i.e. 7.8% of the total number of genes in the data. The performance of MINT (BER of 17.8 on the training set and 11.6 on the test set) was superior than when performing a PLS-DA on the PAM50 genes only (BER of 20.8 on the training set and a very high 75 on the test set). This result shows that the genes selected by MINT offer a complementary characterisation to the PAM50 genes.

In the second analysis, we aimed to provide an alternative signature to the PAM50 genes by ommitting them from the analysis. MINT identified 11, 272 and 253 genes on the first three components respectively (Additional file 1: Table S5 and Figure S10). The genes selected on the first component gradually differentiated Basal, HER2 and Luminal A/B, while the second component genes further differentiated Luminal A from Luminal B (Fig. 3 d). The classification performance was similar in each study (Fig. 3 e), highlighting an excellent reproducibility of the biomarker signature across cohorts and platforms.

Among the 11 genes selected by MINT on the first component, GATA3 is a transcription factor that regulates luminal epithelial cell differentiation in the mammary glands [39, 40], it was found to be implicated in luminal types of breast cancer [41] and was recently investigated for its prognosis significance [42]. The MYB-protein plays an essential role in Haematopoiesis and has been associated to Carcinogenesis [43, 44]. Other genes present in our MINT gene signature include XPB1 [45], AGR3 [46], CCDC170 [47] and TFF3 [48] that were reported as being associated with breast cancer. The remaining genes have not been widely associated with breast cancer. For instance, TBC1D9 has been described as over expressed in cancer patients [49, 50]. DNALI1 was first identified for its role in breast cancer in [51] but there was no report of further investigation. Although AFF3 was never associated to breast cancer, it was recently proposed to play a pivotal role in adrenocortical carcinoma [52]. It is worth noting that these 11 genes were all included in the 30 genes previously selected when the PAM50 genes were included, and are therefore valuable candidates to complement the PAM50 gene signature as well as to further characterise breast cancer subtypes.

Discussion

There is a growing need in the biological and computational community for tools that can integrate data from different microarray platforms with the aim of classifying samples (integrative classification). Although several efficient methods have been proposed to address the unwanted systematic variation when integrating data [4, 7, 9–11], these are usually applied as a pre-processing step before performing classification. Such sequential approach may lead to overfitting and over-optimistic results due to the use of transductive modelling (such as prediction based on ComBat-normalised data [19]) and the use of a test set that is normalised or pre-processed with the training set. To address this crucial issue, we proposed a new Multivariate INTegrative method, MINT, that simultaneously corrects for batch effects, classifies samples and selects the most discriminant biomarkers across studies.

MINT seeks to identify a common projection space for all studies that is defined on a small subset of discriminative variables and that display an analogous discrimination of the samples across studies. Therefore, MINT provides sample plot and classification performance specific to each study (Fig. 3). Among the compared methods, MINT was found to be the fastest and most accurate method to integrate and classify data from different microarray and RNA-seq platforms.

Integrative approaches such as MINT are essential when combining multiple studies of complex data to limit spurious conclusions from any downstream analysis. Current methods showed a high proportion of false positives (44% on MAQC data) and exhibited very poor prediction accuracy (PLS-DA, sPLS-DA and RF, Fig. 2). For instance, RF was ranked second only to MINT on the breast cancer learning set, but it was ranked as the worst method on the test set. This reflects the absence of controlling for batch effects in these methods and supports the argument that assessing the presence of batch effects is a key preliminary step. Failure to do so, as shown in our study, can result in poor reproducibility of results in subsequent studies, and this would not be detected without an independent test set.

We assessed the ability of MINT to identify relevant gene signatures that are reproducible and platform-agnostic. MINT successfully integrated data from the MAQC project by selecting true positives genes that were also differentially expressed in each experiment. We also assessed MINT’s capabilities analysing stem cells and breast cancer data. In these studies, MINT displayed the highest classification accuracy in the training sets and the highest prediction accuracy in the testing sets, when compared to sixteen sequential procedures (Fig. 2). These results suggest that, in addition to being highly predictive, the discriminant variables identified by MINT are also of strong biological relevance.

In the stem cell data, MINT identified 2 genes LIN28A and CAR, to discriminate pluripotent cells (fibroblasts) against non-pluripotent cells (hiPSC and hESC). Pluripotency is well-documented in the literature and OCT4 is currently the main known marker for undifferentiated cells [53–56]. However, MINT did not selected OCT4 on the first component but instead, identified two markers, LIN28A and CAR, that were ranked higher than OCT4 in the DEG list obtained on the concatenated data (see Additional file 1: Figure S11, S12). While the results from MINT still supported OCT4 as a marker of pluripotency, our analysis suggests that LIN28A and CAR are stronger reproducible markers of differentiated cells, and could therefore be superior as substitutions or complements to OCT4. Experimental validation would be required to further assess the potential of LIN28A or CAR as efficient markers.

Several important issues require consideration when dealing with the general task of integrating data. First and foremost, sample classification is crucial and needs to be well defined. This required addressing in analyses with the stem cell and breast cancer studies generated from multiple research groups and different microarray and RNA-seq platforms. For instance, the breast cancer subtype classification relied on the PAM50 intrinsic classifier proposed by [20], which we admit is still controversial in the literature [31]. Similarly, the biological definition of hiPSC differs across research groups [26, 28], which results in poor reproducibility among experiments and makes the integration of stem cell studies challenging [21].

The expertise and exhaustive screening required to homogeneously annotate samples hinders data integration, and because it is a process upstream to the statistical analysis, data integration approaches, including MINT, can not address it.

A second issue in the general process of integrating datasets from different sources is data access and normalisation. As raw data are often not available, this results in integration of data sets that have each been normalised differently, as was the case with the breast cancer data in our study. Despite this limitation, MINT produced satisfactory results in that study. We were also able to overcome this issue in the stem cells data by using the stemformatics resource [29] where we had direct access to homogeneously pre-processed data (background correction, log2- and YuGene-transformed [11]). In general, variation in the normalisation processes of different data sets produces unwanted variation between studies and we recommend this should be avoided if possible.

A final important issue in data integration involves accounting for both between-study differences and platform effects. When samples clustered by study and the studies clustered by platform, then the experimental platform and not the study, is the biggest source of variation (e.g. 75% of the variance in the breast cancer data, Additional file 1: Figure S5). Indeed, there are inherent differences between commercial platforms that greatly magnify unwanted variability, as was discussed by [5] on the MAQC project. As platform information and study effects are nested, MINT and other data integration methods dismiss the platform information and focus on the study effect only. Indeed, each study is considered as included in a single platform. MINT successfully integrated microarray and RNA-seq data, which supports that such an approach will likely be sufficient in most scenarios.

When applying MINT, additional considerations need be taken into account. In order to reduce unwanted systematic variation, the method centers and scales each study as an initial step, similarly to BMC [9]. Therefore, only studies with a sample size >3 can be included, either in a training or test set. In addition, all outcome categories need to be represented in each study. Indeed, neither MINT nor any classification methods can perform satisfactorily in the extreme case where each study only contains a specific outcome category, as the outcome and the study effect can not be distinguished in this specific case.

Conclusion

We introduced MINT, a novel Multivariate INTegrative method, that is the first approach to integrate independent transcriptomics studies from different microarray and RNA-seq platforms by simultaneously, correcting for batch effects, classifying samples and identifying key discriminant variables. We first validated the ability of MINT to select true positives genes when integrating the MAQC data across different platforms. Then, MINT was compared to sixteen sequential approaches and was shown to be the fastest and most accurate method to discriminate and predict three human cell types (human Fibroblasts, human Embryonic Stem Cells and human induced Pluripotent Stem Cells) and four subtypes of breast cancer (Basal, HER2, Luminal A and Luminal B). The gene signatures identified by MINT contained existing and novel biomarkers that were strong candidates for improved characterisation the phenotype of interest. In conclusion, MINT enables reliable integration and analysis of independent genomic data sets, outperforms existing available sequential methods, and identifies reproducible genetic predictors across data sets. MINT is available through the mixMINT module in the mixOmics R-package.