Supervised Regularized Canonical Correlation Analysis: integrating histologic and proteomic measurements for predicting biochemical recurrence following prostate surgery

Abstract

Background

Multimodal data, especially imaging and non-imaging data, is being routinely acquired in the context of disease diagnostics; however, computational challenges have limited the ability to quantitatively integrate imaging and non-imaging data channels with different dimensionalities and scales. To the best of our knowledge relatively few attempts have been made to quantitatively fuse such data to construct classifiers and none have attempted to quantitatively combine histology (imaging) and proteomic (non-imaging) measurements for making diagnostic and prognostic predictions. The objective of this work is to create a common subspace to simultaneously accommodate both the imaging and non-imaging data (and hence data corresponding to different scales and dimensionalities), called a metaspace. This metaspace can be used to build a meta-classifier that produces better classification results than a classifier that is based on a single modality alone. Canonical Correlation Analysis (CCA) and Regularized CCA (RCCA) are statistical techniques that extract correlations between two modes of data to construct a homogeneous, uniform representation of heterogeneous data channels. In this paper, we present a novel modification to CCA and RCCA, Supervised Regularized Canonical Correlation Analysis (SRCCA), that (1) enables the quantitative integration of data from multiple modalities using a feature selection scheme, (2) is regularized, and (3) is computationally cheap. We leverage this SRCCA framework towards the fusion of proteomic and histologic image signatures for identifying prostate cancer patients at the risk of 5 year biochemical recurrence following radical prostatectomy.

Results

A cohort of 19 grade, stage matched prostate cancer patients, all of whom had radical prostatectomy, including 10 of whom had biochemical recurrence within 5 years of surgery and 9 of whom did not, were considered in this study. The aim was to construct a lower fused dimensional metaspace comprising both the histological and proteomic measurements obtained from the site of the dominant nodule on the surgical specimen. In conjunction with SRCCA, a random forest classifier was able to identify prostate cancer patients, who developed biochemical recurrence within 5 years, with a maximum classification accuracy of 93%.

Conclusions

The classifier performance in the SRCCA space was found to be statistically significantly higher compared to the fused data representations obtained, not only from CCA and RCCA, but also two other statistical techniques called Principal Component Analysis and Partial Least Squares Regression. These results suggest that SRCCA is a computationally efficient and a highly accurate scheme for representing multimodal (histologic and proteomic) data in a metaspace and that it could be used to construct fused biomarkers for predicting disease recurrence and prognosis.

Background

With the plentitude of multi-scale, multi-modal, disease pertinent data being routinely acquired for diseases such as breast and prostate cancer, there is an emerging need for powerful data fusion (DF) methods to integrate the multiple orthogonal data streams for the purpose of building diagnostic and prognostic meta-classifiers for disease characterization [1]. Combining data derived from multiple sources has the potential to significantly increase classification performance relative to performance trained on any one modality alone [2]. A major limitation in constructing integrated meta-classifiers that can leverage imaging (histology, MRI) and non-imaging (proteomics, genomics) data streams is having to deal with data representations spread across different scales and dimensionalities [3].

For instance, consider two different data streams F _A (x) and F _B (x) describing the same object x. If F _A (x) and F _B (x) correspond to the same scale or resolution and also have the same dimensionality, then one can envision, concatenating the two data vectors into a single unified vector [F _A (x), F _B (x)] which could then be used to train a classifier. However when F _A (x) and F _B (x) correspond to different scales, resolutions, and dimensionalities, it is not immediately obvious as to how one would go about combining the different types of measurements to build integrated classifiers to make predictions about the class label of x. For instance, directly aggregating data from very different sources without accounting for differences in the number of features and relative scaling, can not only lead to the curse of dimensionality (too many features and not enough corresponding samples [4]), but can lead to classifier bias towards the modality with more attributes. A possible solution is to first project the data streams into a space where the scale and dimensionality differences are removed; a meta-space allowing for a homogeneous, fused, multi-modal data representation.

DF methods try to overcome these obstacles by creating such a metaspace, on which a proper meta-classifier can be constructed. Methods leveraging embedding techniques have been proposed to try and fuse such heterogeneous data for the purpose of classification and prediction [2, 3, 5–7]. However, all of these DF techniques have their own weaknesses in creating an appropriate representation space that can simultaneously accommodate multiple imaging and non-imaging modalities. Generalized Embedding Concatenation [5] is a DF scheme that relies on dimensionality reduction (DR) methods to first eliminate the differences in scales and dimensionalities between the modalities before fusing them. However, these DR methods face the risk of extracting noisy features which degrade the metaspace [8]. Other variants of the embedding fusion idea, including Consensus embedding [6] and Boosted embedding [3] have yielded promising results, but come at a high computational cost. Consensus embedding attempts to combine multiple low dimensional data projections via a majority voting scheme while the Boosted embedding scheme leverages the Adaboost classifier [9] to combine multiple weak embeddings. In the case of weighted multi-kernel embedding using graph embedding [7] and support vector machine classifiers [2], insufficient training data can lead to overfitting and inaccurate weights to the various kernels, which can lower the performance of the meta-classifier [10].

CCA is a statistical DF technique that extracts linear correlations, by using cross-covariance matrices, between 2 data sources, X and Y. It capitalizes on the knowledge that the different modalities represent different sets of descriptors for characterizing the same object. For this reason, the mutual information that is most correlated between the two modalities will provide the most meaningful transformation into a metaspace. In recent years, CCA has been used to fuse heterogeneous data such as pixel values of images and the text attached between these images [11], assets and liabilities in banks [12], and audio and face images of speakers [13].

Regularized CCA (RCCA) is an improved version of CCA which in the presence of insufficient training data prevents overfitting by using a ridge regression optimization scheme [14]. Denote p and q as the number of features in X and Y, and n as the sample size. When n < < p or n < < q, the features in X and Y tend to be highly collinear. This leads to ill-conditioned matrices C _xx and C _yy , which denote the covariance matrix of X with itself and Y with itself, such that their inverses are no longer reliable resulting in an invalid computation of CCA and an unreliable metaspace [15]. The condition placed on the data to guarantee that C _xx and C _yy will be invertible is n ≥ p + q + 1 [16]. However, that condition is usually not met in the bioinformatics domain, where samples (n) are usually limited, and modern technology has enabled very high dimensional data streams to be routinely acquired resulting in very high dimensional feature sets (p and q). This creates a need for regularization, which works by adding small positive quantities to the diagonals of C _xx and C _yy to guarantee their invertibility [17]. RCCA has been used to study expressions of genes measured in liver cells and compare them with concentrations of hepatic fatty acids in mice [18]. However, the regularization process required by RCCA is computationally very expensive. Both CCA and RCCA also fail to take complete advantage of class label information, when available [19].

In this paper, we present a novel efficient Supervised Regularized Canonical Correlation Analysis (SRCCA) DF algorithm that is able to incorporate a supervised feature selection scheme to perform regularization. Mainly, it makes better use of labeled information that in turn allows for significantly better stratification of the data in the metaspace. While SRCCA is more expensive than the overfitting-prone CCA, it provides the needed regularization while also being computationally cheaper than RCCA. SRCCA first produces an embedding of the most correlated data in both modalities via a low dimensional metaspace. This representation is then used in conjunction with a classifier (K-Nearest Neighbor [20] and Random Forest [21] are used in this study) to create a highly accurate meta-classifier.

Along with CCA and RCCA, SRCCA is compared with 2 other low dimensional data representation techniques: Principal Component Analysis (PCA) and Partial Least Squares Regression (PLSR). PCA [22] is a linear DR method that reduces high dimensional data to dominant orthogonal eigenvectors that try to represent the maximal amount of variance in the data. PLSR [23] is a DR method that uses one modality as a set of predictors to try to predict the other modality. Tiwari et al. [24] employed PCA in conjunction with a wavelet based representation of different MRI protocols to build a fused classifier to detect prostate cancer in vivo. PLSR has been used with heterogeneous multivariate signaling data collected from HT-29 human colon carcinoma cells stimulated to undergo programmed cell death to uncover aspects of biological cue-signal-response systems [25].

In this work, we apply SRCCA to the problem of predicting biochemical recurrence in prostate cancer (CaP) patients, following radical prostatectomy, by fusing histologic imaging and proteomic signatures. Biochemical recurrence is commonly defined as a detectable elevation of Prostate Specific Antigen (PSA), a key biomarker for CaP [26–28]. However, the nonspecificity of PSA leads to over-treatment of CaP, resulting in many unnecessary treatments, which are both stressful and costly [29–33]. Even the most widely used prognostic markers such as pathologist assigned Gleason grade [34], which attempts to capture the morphometric and architectural appearance of CaP on histopathology, has been found to be a less than perfect predictor of biochemical recurrence [35]. Additionally, Gleason grade has been found to be subject to inter-, and intra-observer variability [36–38]. While some researchers have proposed quantitative, computerized image analysis approaches [1, 39, 40] for modeling and predicting Gleason grade (a number that goes from 1 to 5 based on morphologic appearance of CaP on histopathology), it is still not clear that an accurate, reproducible grade predictor from histology will also be accurate in predicting biochemical recurrence and long term patient outcome [41].

Recent studies have shown that proteomic markers can be used to predict aggressive CaP [42, 43]. Techniques such as mass spectrometry hold promise in their ability to identify protein expression profiles that might be able to distinguish more aggressive from less aggressive CaP and identify candidates for biochemical recurrence [44–46]. However, more and more, it is becoming apparent that a single prognostic marker may not possess sufficient discriminability to predict patient outcome which suggests that the solution might lie in an integrated fusion of multiple markers [47]. This then begs the question as to what approaches need to be leveraged to quantitatively fuse imaging and non-imaging measurements to build an integrated prognostic marker for CaP recurrence. The overarching goal of this study is to leverage SRCCA to construct a fused quantitative histologic, proteomic marker, and a subsequent meta-classifier, for predicting 5 year biochemical recurrence in CaP patients following surgery.

Our main contributions in this paper are:

• A novel data fusion algorithm, SRCCA, that builds an accurate metaspace representation that can simultaneously represent and accommodate two heterogeneous imaging and non-imaging modalities.

• Leveraging SRCCA to build a meta-classifier to predict risk of 5 year biochemical recurrence in prostate cancer patients following radical prostatectomy by integrating histological image and proteomic features.

The organization of the rest of the paper is as follows: In the methods section, we first review the 4 statistical methods, PCA, PLSR, CCA and RCCA. Next, we introduce our novel algorithm, Supervised Regularized Canonical Correlation Analysis (SRCCA). We then discuss the DF algorithm for metaspace creation and the computational complexities for CCA, RCCA and SRCCA. In the Experimental Design section, we briefly discuss the prostate cancer dataset considered in this study and the subsequent proteomic and histologic feature extraction schemes before moving on to the experiments performed on the dataset where we try to determine the ability of PCA, PLSR, CCA, RCCA and SRCCA to identify patients at risk for biochemical recurrence following surgery. The results are discussed in the subsequent section and the concluding remarks are presented at the end of the paper.

Methods

Review of PCA and PLSR

Principal Component Analysis (PCA) and Partial Least Squares Regression (PLSR) are common statistical methods used to analyze multi-modal data and they are briefly discussed in the following sections. However, further information, explaining how these two methods can be viewed as special cases of the generalized eigenproblem, can be found in [48].

Principal Component Analysis (PCA)

PCA [22] constructs a low dimensional subspace of the data by finding a series of linear orthogonal bases called principal components. Each component seeks to explain the maximal amount of variance in the dataset. Denote two multidimensional variables, X ∈ ℝ^{n × p}and Y ∈ ℝ^n×q, where p and q are the number of features in X and Y and n the number of overall samples. PCA is usually performed on the data matrix, Z ∈ ℝ^n×(p+q), obtained by concatenating the individual modalities such that: Z = [X Y] [24]. $\bar{Z}\in {ℝ}^{n\times \left(p+q\right)}$ is then obtained by subtracting the means of all features for a certain sample from its original feature value in Z so that the resultant $\bar{Z}$ has rows with a 0 mean. $\bar{Z}$ is further broken using singular value decomposition into [22]:

$$\bar{Z}=UE{V}^T$$

(1)

where E ∈ ℝ^n×nis a diagonal matrix containing the eigenvalues of the eigenvectors which are stored in U ∈ ℝ^p×p, and V^T∈ ℝ^m×n. The eigenvalues stored in E explain how much variance of the original $\bar{Z}$ is stored in the corresponding eigenvector, or principal component. Using these eigenvalues as a rank, the top d embedding components can be chosen to best represent the original data in a lower dimensional subspace.

Partial Least Squares Regression(PLSR)

PLSR [49] is a statistical technique that generalizes PCA and multiple regression. The general underlying model behind PLSR is [23]:

$$X=T{P}^T+E$$

(2)

$$Y=T{C}^T+F$$

(3)

where T ∈ ℝ^n×lis a score matrix, P ∈ ℝ^p×land C ∈ ℝ^q×lare loading matrices for X and Y, and E ∈ ℝ^n×pand F ∈ ℝ^n×pare the error terms. PLSR is an iterative process and works by continually approximating, and improving the approximation of the matrices T, P and C [50].

Review of CCA and RCCA

Canonical Correlation Analysis (CCA)

CCA [51] is a way of using cross-covariance matrices to obtain a linear relationship between the two multidimensional variables, X ∈ ℝ^n×pand Y ∈ ℝ^n×q. CCA obtains two directional vectors w _x ∈ ℝ^p×1and w _y ∈ ℝ^q×1such that Xw _x and Yw _y will be maximally correlated. It is defined as the optimization problem [11]:

$$\rho =\underset{w_x,w_y}{max}\frac{w_x^T{C}_{xy}{w}_y}{\sqrt{w_x^T{C}_{xx}{w}_x{w}_y^T{C}_{yy}{w}_y}}$$

(4)

where C _xy ∈ ℝ^p×qis the covariance matrix of the matrices X and Y, C _xx ∈ ℝ^p×pis the covariance matrix of the matrix X with itself and C _yy ∈ ℝ^{q × q}is the covariance matrix of the matrix Y with itself. The solution to CCA reduces to the solution of the following two generalized eigenvalue problems [52]:

$$C_{xy}{C}_{yy}^{-1}{C}_{yx}=\lambda C_{xx}{w}_x$$

(5)

$$C_{yx}{C}_{xx}^{-1}{C}_{xy}=\lambda C_{yy}{w}_y$$

(6)

where λ is the generalized eigenvalue representing the canonical correlation, and w _x and w _y are the corresponding generalized eigenvectors. CCA can further produce exactly min{p, q) orthogonal embedding components (sets of w _x X and w _y Y) which can be sorted in order of decreasing correlation, λ.

Regularized Canonical Correlation Analysis (RCCA)

RCCA [53, 54] corrects for noise in X and Y by first assuming that X and Y are contaminated with noise, N _x ∈ ℝ^n×pand N _Y ∈ ℝ^n×q. We assume that these noise vectors in the p and q columns of N _X and N _Y , respectively, are gaussian, independent and identically distributed. For this reason, all combinations of the covariances of the p columns of N _X and q columns of N _Y will be 0 except the covariance of a particular column vector with itself. This variance of each column of N _X and N _Y is labeled λ _x and λ _y and these labels are called the regularization parameters. The matrix C _xy will not be affected but the matrices C _xx and C _yy become C _xx + λ _x I _x and C _yy + λ _x I _x . The solution to RCCA now becomes the solution to these generalized eigenvalue problems [52]:

$$C_{xy}{\left(C_{yy}+{\lambda }_y{I}_y\right)}^{-1}{C}_{yx}=\lambda \left(C_{xx}+{\lambda }_x{I}_x\right)w_x$$

(7)

$$C_{yx}{\left(C_{xx}+{\lambda }_x{I}_x\right)}^{-1}{C}_{xy}=\lambda \left(C_{yy}+{\lambda }_y{I}_y\right)w_y$$

(8)

The regularization parameters next have to be chosen. For i ∈ {1, 2, . . . , n}, let $w_x^i$ and $w_y^i$ denote the weights calculated from RCCA when samples X _i and Y _i are removed. λ _x and λ _y are varied in a certain range θ₁ ≤ λ _x , λ _y ≤ θ₂ and chosen via a grid search [55] optimization of the following cost function [18]:

$$\underset{{\lambda }_x,{\lambda }_y}{max}\left[corr\left({\left\{X_i{w}_x^i\right\}}_{i=1}^n,{\left\{Y_i{w}_y^i\right\}}_{i=1}^n\right)\right]$$

(9)

where corr (·, ·) refers to the Pearson's correlation coefficient [56]. The above cost function essentially measures the change in the produced $w_x^i$ and $w_y^i$ when a sample i is omitted and seeks the optimal λ _x and λ _y where this change is minimized. λ _x and λ _y are chosen using the embedding component with the highest λ and then adjusted for the remaining dimensions [18].

Extending RCCA to SRCCA

Supervised Regularized Canonical Correlation Analysis (SRCCA) chooses λ _x and λ _y using a supervised feature selection method (t-test, Wilcoxon Rank Sum Test and Wilks Lambda Test are used in this study). Denote ${\mathbb{W}}_{\mathsf{\text{1}}}$ and ${\mathbb{W}}_2$ as class 1 and class 2 and μ₁ and μ₂, ${\sigma }_1^2$ and ${\sigma }_2^2$, n₁ and n₂ as the means, variances, and sample sizes of ${\mathbb{W}}_{\mathsf{\text{1}}}$ and ${\mathbb{W}}_2$. The data in the metaspace, Xw _x or Yw _y , can be split using its labels into the n₁ samples that belong to ${\mathbb{W}}_{\mathsf{\text{1}}}$ and the n₂ samples that belong to class ${\mathbb{W}}_2$, where n₁ + n₂ = n. These two partitions can then be used to calculate the discrimination level between the samples of the two classes in the metaspace representation. In this study, we implement RCCA with the t-test (SRCCA _TT ), the Wilcoxon Rank Sum Test (SRCCA _WRST ) and the Wilks Lambda Test (SRCCA _WLT ) to try to choose more appropriate regularization parameters, λ _x and λ _y , that can more successfully stratify the samples in the metaspace compared to the parameters chosen by RCCA. Similar to RCCA, for SRCCA, λ _x and λ _y are chosen using the embedding component with the most discriminatory score as chosen by the feature selection schemes below and then adjusted for the remaining dimensions.

SRCCA_TT

The t-test [57] is a parametric test that assumes the distributions of the two samples are normal and tests whether these distributions have the same means. The t-score, which measures the number of standard deviations the two means of n₁ samples of ${\mathbb{W}}_{\mathsf{\text{1}}}$ and n₂ samples of ${\mathbb{W}}_2$ are away from each other, is maximized using a grid search algorithm as:

$$\underset{{\lambda }_x,{\lambda }_y}{max}\frac{\left|\right|{\mu }_1-{\mu }_2\left|\right|}{\sqrt{\frac{{\sigma }_1^2}{n_1}+\frac{{\sigma }_2^2}{n_2}}}.$$

(10)

SRCCA_WRST

Wilcoxon Rank Sum Test [58] sorts both the samples in order from lowest value to highest value. It then uses their respective ranks within the population to calculate the discriminatory score:

$$\underset{{\lambda }_x,{\lambda }_y}{max}\left\{\left(\sum _{i=1}^{n_2}{b}_i-\frac{n_2\left(n_2+1\right)}{2}\right),\left(n_1{n}_2-\sum _{i=1}^{n_2}{b}_i+\frac{n_2\left(n_2+1\right)}{2}\right)\right\},$$

(11)

where b _i represents the rank of the sample $i\in {\mathbb{W}}_2$ with respect to the rest of the samples.

SRCCA_{W LT}

In an ideal metaspace representation, samples from each class will be grouped together while the samples from different classes will be grouped separately. The WLT [59] capitalizes on this knowledge and calculates the ratio of within class variance of both samples to the total variance of both samples combined. Wilks Lambda (Λ) is minimized using a grid search algorithm as:

$$\underset{{\lambda }_x,{\lambda }_y}{min}\frac{n_1{\sigma }_1^2+n_2{\sigma }_2^2}{n{\sigma }^2}.$$

(12)

Data Fusion in the context of CCA, RCCA and SRCCA

DF is performed as described in Foster et al. [60]. When the Xw _x and Yw _y are maximally correlated, each modality represents similar information, and thus either Xw _x or Yw _y can be used to represent the original two modalities in the metaspace. Moreover, X and Y are both descriptors of the same object and thus, the most relevant information is the data that exists and is correlated in both modalities. Thus, a high correlation of Xw _x and Yw _y is indicative that meaningful data, measuring the object of interest, is being added to the metaspace.

In order of decreasing λ, the top d embedding components, up to φ = min{p, q} can be chosen to represent the two modalities in a metaspace. However, the lower embedding components will have a lower λ, and thus a lower correlation between Xw _x and Yw _y which might imply that non-relevant data is being added to the metaspace. To avoid this issue, a threshold, λ₀, can be selected such that only embedding components with λ ≥ λ₀ will be included in the metaspace.

Computational Complexity

Given φ = min{p, q}, CCA has a computational complexity of φ! (based on the source code in [61]). The regularization algorithm requires a grid search process for each ordered pair (λ _x , λ _y ). Assume v potential λ _x and λ _y sampled evenly between θ₁ and θ₂. RCCA requires a training/testing cross-validation strategy, at each ordered pair (λ _x , λ _y ), to find the optimal λ _x and λ _y . It will require CCA to be performed an order of n times at each of the v intervals leading to a complexity of vnφ!. SRCCA only requires a CCA factorization once at each of the v intervals leading to a complexity of vφ!.

The computational complexities for each of the CCA schemes are summarized in Table 1. Table 1 indicates that SRCCA is an order of n times faster compared to RCCA. However, SRCCA is also more complex compared to CCA and will have a longer execution time.

Table 1 The computational complexities of all 3 DF methods used in this study

Experimental Design

Data Description

A total of 19 prostate cancer patients at the Hospital at the University of Pennsylvania were considered for this study. All patient identifiers are stripped from the data at the time of acquisition. The data was deemed to be exempt for review by the internal review board at Rutgers University and the protocol was approved by the University of Pennsylvania internal review board. Hence, the data was deemed eligible for use in this study. All of these patients had been found to have prostate cancer on needle core biopsy and subsequently underwent radical prostatectomy. 10 of these patients had biochemical recurrence within 5 years following surgery (BR) and the other 9 did not (NO BR). The 19 patient studies were randomly chosen from a larger cohort of 110 patient studies at the University of Pennsylvania all of whom had been stage and grade matched (Gleason score of 6 or 7) and had undergone gland resection. Of these 110 cases, 55 had experienced biochemical recurrence within 5 years while the other 55 had not. The cost of the mass spectrometry to acquire the proteomic data limited this study to only 19 patient samples. Following gland resection, the gland was sectioned into a series of histological slices with a meat cutter. For each of the 19 patient studies, a representative histology section on which the dominant tumor nodule was observable was identified. Mass Spectrometry was performed at this site to yield a protein expression vector. The representative histologic sections were then digitized at 40 × magnification using a whole slide digital scanner.

In the next two sections, we briefly describe the construction of the proteomic and histologic feature spaces. Subsequently we describe the strategy for combination of quantitative image descriptors from the tumor site on the histological prostatectomy specimen and the corresponding proteomic measurements obtained from the same tumor site, via mass spectrometry. The resultant meta-classifier, constructed in the fused meta-space, is then used to distinguish the patients at 5 year risk of biochemical recurrence following radical prostatectomy from those who are not.

Proteomic Feature Selection

Prostate slides were deparaffinized, and rehydrated essentially as described in [62]. Tumor areas previously defined on a serial H&E section were collected by needle dissection, and formalin cross-links were removed by heating at 99°C. The FASP (Filter-Aided Sample Preparation) method [63] was then used for buffer exchange and tryptic digest. After peptide purification on C-18 StageTips [64] samples were analyzed using nanoflow C-18 reverse phase liquid chromatography/tandem mass spectrometry (nLC-MS/MS) on an LTQ Orbitrap mass spectrometer. A top-5 data-dependent methodology was used for MS/MS acquisition, and data files were processed using the Rosetta Elucidator proteomics package, which is a label-free quantitation package that uses extracted ion chromatograms to calculate protein abundance rather than peptide counts. A high dimensional feature vector was obtained, denoted ϕ^P∈ ℝ^{19 × 953}, characterizing each patient's protein expression profile following surgery. This data underwent quantile normalization, log(2) transformation, and mean and variance normalization on a per-protein basis.

Quantitative Histologic Feature Extraction

In prostate whole-mount histology, denoted ϕ^H∈ ℝ^{19 × 151} (Figure 1 (a), (f)), the objects of interest are the glands (shown in Figure 1 (b), (g)), whose shape and arrangement are highly correlated with cancer progression [1, 39, 65, 66]. We briefly describe this process below. Prior to extracting image features, we employ an automatic region-growing gland segmentation algorithm presented by Monaco et al. [67]. The boundaries of the interior gland lumen and the centroids of each gland, allow for extraction of 1) morphological and 2) architectural features from histology as described briefly below. More extensive details on these methods are in our other publications [5, 39, 68].

Figure 1

Multi-modal patient data (top row: relapsed case, bottom row: non-relapsed case). (a), (f) Original prostate histology section showing region of interest, (b), (g) Magnified ROI showing gland segmentation boundaries, (c), (h) Voronoi Diagram (d), (i) Delaunay Triangulation depicting gland architecture, (e), (j) Plot of the proteomic profile obtained from the dominant tumor nodule regions (white box in (a), (f) respectively) via mass spectrometry.

Glandular Morphology

The set of 100 morphological features [1], (denoted ϕ^M∈ ℝ^{19 × 100}), of attributes, consists of the average, median, standard deviation, and min/max ratio for features such as gland area, maximum area, area ratio, and estimated boundary length (See Table 2).

Table 2 Description of 25 Proteomic Features, 100 Morphological, and 51 Architectural

Architectural Feature Extraction

51 architectural image features, which have been shown to be predictors of cancer [69], (denoted ϕ^A∈ ℝ^{19 × 51}), were extracted in order to quantify the arrangement of glands present in the section (See Table 2). Voronoi diagrams, Delaunay Triangulation and Minimum Spanning Trees were constructed on the digital histologic image using the gland centroids as vertices, the gland centroids having previously been identified via the scheme in [68].

Fusing Proteomic, Histologic Features for Predicting Biochemical Recurrence in CaP Patients Post-Surgery

Experiment 1 - Comparing SRCCA with CCA and RCCA

We performed CCA, RCCA, and SRCCA on selected multimodal combinations, ϕ^Pand ϕ^J, where J ∈ {M, A, H}. ϕ^Pwas reduced to 25 features as ranked by the t-test, with a p-value cutoff of p = .05, using a leave-one-out validation strategy. For CCA, ϕ^Pand ϕ^Jwere used as the two multidimensional variables, X and Y, as mentioned above in Section 2. For RCCA and SRCCA, ϕ^Pand ϕ^Jwere used in a manner similar to CCA except they are tested with regularization parameters λ _x and λ _y evenly spaced from θ₁ = .001 to θ₂ = .2 with v = 200.

The top d = 3 embedding components (which were experimentally found to meet the criteria of λ₀ = .99 for all SRCCA on all 3 multimodal combinations) were produced from CCA, RCCA, SRCCA _TT , SRCCA _WRST , and SRCCA _WLT . The classification accuracies were determined with the classifiers K-Nearest Neighbor, denoted via ϕ^KNN[20], with K = 1, and Random Forest, denoted via ϕ^RF[21], with 50 Trees. Both these classifiers were used because of their high computational speed. Accuracies were determined using leave-one-out validation, which was implemented because of the small sample size. In this process, 18 samples were used for the initial feature pruning, determining the optimal regularization parameter and training the classifier while the remaining sample was used as the testing set for evaluating the classifier. This procedure was repeated till all the samples were used in the testing set.

Experiment 2 - Comparing SRCCA with PCA and PLSR

In addition to the steps performed in Experiment 1, metaspaces were also produced with PCA and PLSR. ϕ^Pand ϕ^Jwere concatenated and PCA was then performed on this new data matrix. For PLSR, a regression of ϕ^Jon ϕ^Pwas performed.

Similarly, using the top d = 3 embedding components produced from PCA, PLSR, SRCCA _TT , SRCCA _WRST , and SRCCA _WLT , the classification accuracies of ϕ^KNN, with K = 1, and ϕ^RF, with 50 Trees, were determined using leave-one-out validation.

Experiment 3 - Comparing classifier accuracy for PCA, PLSR and CCA variants using metaspace representations

Using the 10 different values for d ∈ {1, 2, ..10}, and the 3 fusion schemes considered (ϕ^P, ϕ^M), (ϕ^P, ϕ^A), and (ϕ^P, ϕ^H), 30 different embeddings were obtained for PCA, PLSR, CCA, RCCA, SRCCA _TT , SRCCA _WRST , and SRCCA _WLT . The maximum and median of these 30 different measurements for each classifier were calculated.

In addition, we denote as α₁(i), the classification accuracy obtained by the DF scheme i, where i ∈ {PCA, PLSR, CCA, RCCA} and α₂(j) as the accuracy obtained by the DF scheme i, where j ∈ {SRCCA _TT , SRCCA _WRST , SRCCA _WLT }. A two paired student t-test was employed to identify whether there were statistically significant improvements in the 3 SRCCA variants by comparing the classification accuracies with the null hypothesis:

$$H_o:{\alpha }_1\left(i\right)={\alpha }_2\left(j\right)$$

(13)

for all i ∈ {PCA, PLSR, CCA, RCCA} and for all j ∈ {SRCCA _TT , SRCCA _WRST , SRCCA _WLT }.

Experiment 4 - Computational consideration for RCCA and SRCCA

We measured the 3 individual single run completion times for RCCA and SRCCA to fuse (ϕ^P, ϕ^M), (ϕ^P, ϕ^A), and (ϕ^P, ϕ^H), with the null hypothesis:

$$H_o:\mathsf{\text{completion}}\mathsf{\text{time}}\mathsf{\text{of}}\mathsf{\text{RCCA}}\mathsf{\text{ = }}\mathsf{\text{completion}}\mathsf{\text{time}}\mathsf{\text{of}}\mathsf{\text{SRCCA}}$$

(14)

These experiments were performed on a quadcore computer with a clock speed of 1.8 GHz, and the programs were written on MATLAB(R) platform.

Results and Discussion

Experiment 1

Across both classifiers for d = 3, the 3 SRCCA variants, SRCCA _TT , SRCCA _WRST , SRCCA _WLT , had a combined median classification accuracy of 80% compared to 60% for CCA and 42% for RCCA. SRCCA also performed better in all 36 of 36 direct comparisons with CCA and RCCA (see Tables 3 and 4). The higher classification accuracy results indicate that SRCCA produces a metaspace, where the samples are more stratified, compared to CCA and RCCA. This also seems to indicate that the supervised scheme of choosing regularization parameters, by the 3 SRCCA variants, is a more appropriate scheme for classification purposes compared to the ridge regression scheme used by RCCA.

Table 3 Experiment 1: Classification Accuracy with K-Nearest Neighbor

Table 4 Experiment 1: Classification Accuracy with Random Forest

These results, which seem to suggest that SRCCA outperforms the other two CCA based approaches for this dataset, CCA and RCCA, are observable in the embedding plots of Figure 2, which show the metaspace produced by CCA, RCCA, SRCCA _TT , SRCCA _WRST and SRCCA _WLT with d = 2 components. It may be seen that because CCA lacks regularization, the corresponding covariance matrices are singular and lack inverses. For this reason, in Figure 2 the embedding components are not orthogonal but are highly correlated to each other and yield the same information. RCCA overcomes this regularization problem but still does not produce the same level of discrimination between patient classes compared to the 3 variations of SRCCA. Note that SRCCA _TT , SRCCA _WRST and SRCCA _WLT chose similar regularization parameters, λ _x and λ _y , and have similar embedding plots.

Figure 2

2-dimensional representation of ( ϕ^P, ϕ^A). 2-dimensional representation of (ϕ^P, ϕ^A) using (a) CCA, (b) RCCA, (c) SRCCA _TT , (d) SRCCA _WRST and (e) SRCCA _WLT where the X and Y axes are the two most significant embedding components produced by the 3 different algorithms. CCA (a) suffers from lack of regularization, RCCA (b) is regularized but does not produce the best metaspace while the three variations of SRCCA (c)(d)(e) result in the best embedding components in terms of classification accuracy distinguished via best fit ellipses with one outlier.

Experiment 2

We see that SRCCA _TT , SRCCA _WRST , SRCCA _WLT are able to outperform PCA and PLSR in all 36 of 36 direct comparisons (see Tables 5 and 6). Even though, across both classifiers for d = 3, PCA and PLSR have median classification accuracies of 64% and 61%, which is higher than the accuracies for CCA and RCCA, it is still much lower than the 80% for SRCCA _TT , SRCCA _WRST , SRCCA _WLT . These results also seem to indicate that SRCCA _TT , SRCCA _WRST , SRCCA _WLT could also create a more appropriate metaspace than, not only CCA and RCCA, but also PCA and PLSR.

Table 5 Experiment 2: Classification Accuracy with K-Nearest Neighbor

Table 6 Experiment 2: Classification Accuracy with Random Forest

Experiment 3

In Tables 7 and 8 we see that the maximum and median ϕ^KNNand ϕ^RFof the 3 SRCCA variants for fusion of (ϕ^I, ϕ^J) were much higher than the corresponding values of PCA, PLSR, CCA or RCCA. We also see that SRCCA _WLT attains a maximum classifier accuracy of 93.16% (see Table 7). In Tables 9 and 10, the 3 SRCCA variants are statistically significantly better than PCA, PLSR, CCA or RCCA even at the p = .001 level using either classifiers, ϕ^KNNor ϕ^RF. We further see that SRCCA _WLT tends to marginally outperform SRCCA _TT and SRCCA _WRST . However given the small sample size it is difficult to draw any definitive conclusions about which of SRCCA _TT , SRCCA _WRST , or SRCCA _WLT might be the better SRCCA variant.

Table 7 Experiment 3: Maximum ϕ^KNNand ϕ^RFof DF schemes across d ∈ {1, 2, ..10}

Table 8 Experiment 3: Median ϕ^KNNand ϕ^RFof DF schemes across d ∈ {1, 2, ..10}

Table 9 Experiment 3: Statistical Significance (p-value) of SRCCA for ϕ^KNN

Table 10 Experiment 3: Statistical Significance (p-value) of SRCCA for ϕ^RF

In Figures 3 and 4, we see the classification accuracies of the 7 DF methods, PCA, PLSR, CCA, RCCA, SRCCA _TT , SRCCA _WRST , or SRCCA _WLT over a range of d ∈ {1, 2, ..10} embedding components for the fusion (ϕ^P, ϕ^H). Importantly, we see that the SRCCA _TT , SRCCA _WRST , and SRCCA _WLT all outperform PCA, PLSR, CCA and RCCA for a majority of the embedding dimensions, across both the ϕ^KNNand ϕ^RFclassifiers.

Figure 3

Classification accuracies of ( ϕ^P, ϕ^H) across dimensions d ∈ {1, 2, ..10} using the classifier ϕ^KNN. Accuracies were obtained for fusing (ϕ^P, ϕ^H), with PCA, PLSR, CCA, RCCA, SRCCA _TT , SRCCA _WRST , and SRCCA _WLT across d ∈ {1, 2, ...10} components, using ϕ^KNN, with K = 1, and leave-one-out validation to identify patients at the risk of biochemical recurrence from those who are not.

Figure 4

Classification accuracies of ( ϕ^P, ϕ^H) across dimensions d ∈ {1, 2, ..10} using the classifier ϕ^RF. Accuracies were obtained for fusing (ϕ^P, ϕ^H), with PCA, PLSR, CCA, RCCA, SRCCA _TT , SRCCA _WRST , and SRCCA _WLT across d ∈ {1, 2, ...10} components, using ϕ^KNN, with K = 1, and leave-one-out validation to identify patients at the risk of biochemical recurrence from those who are not.

Experiment 4

Figure 5 reveals that the completion time of SRCCA is significantly lower than the completion time of RCCA. Even though the differences in these times are visibly different, a p-value of 1.9 × 10^-3 even with just 3 samples, indicates that SRCCA appears to be statistically significantly faster compared to RCCA.

Figure 5

Computational run times for SRCCA and RCCA for fusing ( ϕ^P, ϕ^M), ( ϕ^P, ϕ^A), and ( ϕ^P, ϕ^H). SRCCA significantly outperforms RCCA across all fusion experiments. SRCCA significantly outperforms RCCA across all fusion experiments.

Note that the canonical factorization stage is the most time consuming part of the of the algorithm. The Feature Selection stage computation, in comparison, is not as time consuming. SRCCA _TT , SRCCA _WRST , and SRCCA _WLT (whose results are reported in Figure 5) all have similar execution times.

Conclusions

In this paper, we presented a novel data fusion (DF) algorithm called Supervised Regularized Canonical Correlation Analysis (SRCCA) that, unlike CCA and RCCA, is (1) able to fuse with a feature selection (FS) scheme, (2) regularized, and (3) computationally cheap. We demonstrate how SRCCA can be used for quantitative integration and representation of multi-scale, multi-modal imaging and non-imaging data. In this work we leveraged SRCCA for the purpose of constructing a fused quantitative histologic-proteomic classifier for predicting which prostate cancer patients are at risk for 5 year biochemical recurrence following surgery. We have demonstrated that SRCCA is able to (1) produce a metaspace, where the samples are more stratified than the metaspace produced by CCA or RCCA, (2) better identify patients at the risk of biochemical recurrence compared to Principal Component Analysis (PCA), Partial Least Squares Regression (PLSR), CCA or RCCA, (3) perform regularization, all the while being statistically significantly faster compared to RCCA.

While the fused prognostic classifier for predicting biochemical recurrence in this work appears to be promising, we also acknowledge the limitations of this work: (1) As previously mentioned, the cost of mass spectrometry limited this study to only 19 datasets. By using a minimum sample size derivation model [70, 71], we were able to determine that our fused SRCCA classifier would yield an accuracy of 93%, more than 95% of the time if our dataset were expanded to 56 studies. We intend to evaluate our classifier on such a cohort in the future. (2) Ideally, a randomized cross validation strategy should have been employed for the training and evaluation of the classifier. Unfortunately, this was also limited by the size of the cohort. While both parametric and non-parametric feature selection strategies were employed in this work, the availability of a larger dataset for classification in conjunction with SRCCA would allow for employment of parametric selection strategies, assuming that the underlying distribution can be estimated. For small sample datasets, a non-parametric feature selection strategy might be more approrpriate. In future work, we also plan to apply SRCCA in the context of data fusion for other imaging and non-imaging datasets in the context of other problem domains and applications.