nipalsMCIA: Flexible Multi-Block Dimensionality Reduction in R via Non-linear Iterative Partial Least Squares

Motivation: With the increased reliance on multi-omics data for bulk and single cell analyses, the availability of robust approaches to perform unsupervised analysis for clustering, visualization, and feature selection is imperative. Joint dimensionality reduction methods can be applied to multi-omics datasets to derive a global sample embedding analogous to single-omic techniques such as Principal Components Analysis (PCA). Multiple co-inertia analysis (MCIA) is a method for joint dimensionality reduction that maximizes the covariance between block- and global-level embeddings. Current implementations for MCIA are not optimized for large datasets such such as those arising from single cell studies, and lack capabilities with respect to embedding new data. Results: We introduce nipalsMCIA, an MCIA implementation that solves the objective function using an extension to Non-linear Iterative Partial Least Squares (NIPALS), and shows significant speed-up over earlier implementations that rely on eigendecompositions for single cell multi-omics data. It also removes the dependence on an eigendecomposition for calculating the variance explained, and allows users to perform out-of-sample embedding for new data. nipalsMCIA provides users with a variety of pre-processing and parameter options, as well as ease of functionality for down-stream analysis of single-omic and global-embedding factors. Availability: nipalsMCIA is available as a BioConductor package at https://bioconductor.org/packages/release/bioc/html/nipalsMCIA.html, and includes detailed documentation and application vignettes. Supplementary Materials are available online.


Introduction
Multiple co-inertia analysis (MCIA) is a member of the family of joint dimensionality reduction (jDR) methods that extend unsupervised dimension reduction techniques such as Principal Components Analysis (PCA) and Non-negative Matrix Factorization (NMF) to datasets with multiple data blocks (alternatively called views) [1,2].Such methods, also known as multi-block or multi-view analysis algorithms, are becoming increasingly important in the field of bioinformatics, where data is often collected simultaneously using multiple omics technologies such as transcriptomics, proteomics, epigenomics, metabolomics, etc. [3].
Here, we present a new implementation in R/Bioconductor of MCIA, nipalsMCIA, that uses an extension with proof of monotonic convergence of Non-linear Iterative Partial Least Squares (NIPALS) to solve the MCIA optimization problem [4].This implementation shows significant speed-up over existing Singular Value Decomposition (SVD)-based approaches for MCIA [5,6] on large datasets.Furthermore, nipalsMCIA offers users several options for pre-processing and deflation to customize algorithm performance, methodology to perform out-of-sample global embedding, and analysis and visualization capabilities for efficient results interpretation.We show application of nipalsMCIA to both bulk and single cell multi-omics data.The overall workflow that includes the optimization steps and analyses for nipalsMCIA is outlined in Figure 1.

Notation and preliminaries
Scalars, vectors, and matrices are represented in lowercase script (), lowercase script with a vector symbol ( ⃗ ) and bold uppercase script (), respectively.The th column vector of a matrix  is denoted ⃗  () .
Since we are evaluating several datasets (termed blocks) simultaneously, the sample-by-feature data matrix for the th block is labeled as   .We denote the column-wise concatenation of  data blocks as the 'global' data matrix  = [ 1 |...|  ].

Loadings and Scores
MCIA extends the concept from PCA of deriving principal components (which we term scores) and loadings (which we also term loadings) to the multi-block setting.The loadings are a set of optimal axes in feature space, while the scores are the projection coefficients of the samples onto these axes.Unlike PCA, MCIA generates two types of scores and loadings, one set for all the data (global scores/loadings), and the other for the individual omics (block scores/loadings).The number of scores/loadings generated is equal to the dimension of the MCIA embedding of the data, which we will denote as .
Originally, the optimization criteria for MCIA were presented using the concept of statistical triplets [7,5].The criteria can equivalently be represented as a parameterized member of the Regularized Canonical Correlation Analysis (RGCCA) family of multi-variate dimension reduction methods [2,8], which is consistent with the optimization criteria that is solved by an extension of the NIPALS algorithm [4].We review these criteria below.
Scores and loadings are computed by nipalsMCIA to satisfy the objective function and orthogonality constraints where ⃗  () = ( )  is a vector of block contributions to the  th order global score, with constraint || ⃗  () || 2 = 1 for all orders  = 1, ...,  as in [4], and   is the Kronecker delta function.Equation ( 2) is solved separately for each order () up to the dimension of the embedding, .The block scores { ⃗ , ⃗ , ..., ⃗  ()  } represent a -dimensional embedding of the samples in the orthonormal set of block loadings vectors for block .This contrasts with Consensus PCA (CPCA), which solves for the same objective function as MCIA, but with an orthogonality constraint on the global scores instead of the block loadings [9].In nipalsMCIA, users can choose to use either method.

NIPALS strategy for computing MCIA
Several methods exist for computing MCIA, including direct computation from the principal components of the covariance matrix (see [2]).The implementation in nipalsMCIA uses an extension of Nonlinear Iterative Partial Least Squares method (NIPALS) [4].NI-PALS was first introduced as an iterative (power) method to estimate principal components [10,11], and later extended to the multi-block setting [12].A modification of the multi-block algorithm was proven to have monotonic convergence [4].Since the NIPALS procedure is iterative, it does not require a full eigendecomposition.Moreover, it easily enables a choice of deflation methods.In nipalsMCIA, the stable multi-block extension to NIPALS [4] is implemented with deflation options for both MCIA and CPCA.Additionally, variance explained by each component is also calculated without reference to an eigendecomposition calculation.

Usage and functionality
Since MCIA is designed to handle multiple omics data blocks, preprocessing options are available both at within-and whole-block levels.The latter is recommended to account for potential disparities

Analysis & Interpretation
The nipals_multiblock function is used to run MCIA in nipalsMCIA.The function outputs an object of the NipalsResult class, which includes the global scores and loadings, block scores and loadings, the global score eigenvalues, and the block score contributions vector for all orders up to the maximum specified via the num_PCs argument.The global scores represent the projection of the multi-block data in the reduced space, and can be plotted with or without corresponding block scores (Figure 1C, ii).The contribution of each block to the global score can be easily visualized (Figure 1C, iii), along with high-scoring features (Figure 1C, iv).
Vignettes providing full analysis pipelines using nipalsMCIA for bulk and single cell data are available with the package.The example bulk data is a subset of the National Cancer Institute 60 tumor-cell line screen (NCI60 data) [13,8].It includes RNA-Seq, miRNA, and protein data from 21 cell lines that correspond to three cancer subtypes (brain, leukemia, and melanoma).The single cell data is sourced from 10x Genomics and includes both gene expression and cell surface antibody data [14].The single cell analysis vignette includes instruction on how to obtain, process, and prepare the dataset for nipalsMCIA, along with a demonstration of the capability of nipalsMCIA for effectively clustering known cell types in a computationally efficient manner.

Out-of-sample embedding
The loadings vectors generated by MCIA on a dataset  represent linear combinations of the original features of .Therefore, after computing MCIA on a training dataset, one can use the associated loadings vectors to predict global embeddings for a test dataset of new observations of the same features.nipalsMCIA provides the predict_gs function for this task.
This can be valuable for testing the quality of the embedding, as well as embedding new data without rerunning the decomposition.We provide a vignette in the package showing how this can be done using the NCI60 data set, using 70% of the data to train the model, and then deriving global scores for the remaining 30%.

Computation time comparison for MCIA algorithms
We used three datasets to compare the performance of nipalsMCIA with two other implementations of MCIA: MOGSA [6], and Omicade [5].The three datasets are composed of the NCI60 data, the 10x single cell data filtered for the top 2000 most variable genes, and the same single cell data without filtering.Data pre-processing was standardized across all algorithms and a decomposition for 10 factors was performed across all datasets and implementations.All experiments were performed in R 4.3.0 on a MacBook with 3.2GHz and 16GB RAM.The dimensions of the datasets and performance are shown in Table 1.We observe that while MOGSA has slightly faster performance than nipalsMCIA and Omicade on the smaller NCI60 dataset, nipalsMCIA is an order of magnitude faster for both the filtered and full single cell datasets, even when using the 'fast SVD' option in MOGSA.The speedup offered by nipalsMCIA thus opens up capabilities for practical deployment of nipalsMCIA on a larger variety of datasets, including high-dimensional single cell data.

Discussion
The accessibility of next-generation sequencing and other highthroughput biological assays are resulting in an increase of multiblock (or multi-modal) datasets [15,16,17,18].Analysis of these data are facilitated by the application of joint dimensionality reduction methods such as MCIA.nipalsMCIA is a comprehensive R package that implements MCIA in a highly efficient manner using the NIPALS algorithm.The package features various pre-processing and analysis options, is much faster for large input datsets compared with existing packages, supports the projection for out-of-sample scores, and offers visualization options for scores and top-magnitude loadings at each order.

)Figure 1 .
Figure 1.Workflow overview for nipalsMCIA performed on the three-block NCI60 data from the main text.a) A breakdown of the NIPALS algorithm for performing MCIA.Data blocks are normalized before scores and loadings are computed to satisfy the objective function.Higher-order results are then computed after the data has been deflated with the current scores or loadings.b) Scree plot for the proportion of variance explained by each order of global score/loading.c) Scheme for interpreting the global loadings and scores.(i) Global scores are calculated from the global data matrix and global loadings.(ii) Global scores represent low-dimensional embeddings of the data used to cluster samples via hierarchical clustering.Colors represent the three different cancer types associated with each sample (iii) Block contributions vectors plotted to visualize the weight of each block on each order of global score.(iv) The first global loadings vector is plotted to identify the top features for the first global score.
Mattessich et al. nipalsMCIA: Flexible Multi-Block Dimensionality Reduction in R via Nonlinear Iterative Partial Least Squares in block size.

Table 1 .
Computation time (in seconds) comparison for different MCIA implementations and datasets.