Jackstraw Inference for AJIVE Data Integration

Following the standard in bioinformatics, the columns of these matrices represent data vectors or different experimental units and the rows represents features sometimes called variables or traits. In particular, all data matrices have the same number of columns (experimental units) but potentially different number of rows (features). This is the transpose of what is usual in classical statistics. Using this convention, AJIVE [7] assumes the following model for each of the mean centered data blocks:

Here the low rank matrix of joint scores $V_{J}$ is common for all the data blocks and is composed of $r_{J}$ orthonormal columns. The low rank matrices of individual scores $V_{I}^{m}$ are composed of $r_{I}^{m}$ orthonormal columns that are orthogonal to $V_{J}$. Notice that opposite to the usual PCA convention, to give the AJIVE analysis scale free properties, it is the scores that are normalized to have norm 1 and the loadings (columns of $L_{J}^{m}$ and $L_{I}^{m}$) reflect the variation quantified by singular values. The matrix $V_{J}$ is determined only up to a left multiplication by an orthogonal matrix. The matrices $V_{I}^{m}$ are chosen so that, just as in SVD, the columns of $L_{I}^{m}$ are orthogonal and decreasing in norm. The $E^{m}$ are full rank residual noise matrices.

The AJIVE modes of variation are obtained as outer products of the corresponding columns of the loading and score matrices, i.e., if $L_{i,J}^{m}$ and $V_{i,J}$, $i=1,\ldots,r_{J}$, are the $i$th column of $L_{j}^{m}$ and $V_{J}$ respectively, the $i$th joint mode of variation for block $m$ is the rank 1 matrix $L_{i,J}^{m}\times V_{i,J}^{\top}$. Analogously, the $j$th individual mode of variation is the rank 1 matrix $L_{j,I}^{m}\times V_{j,I}^{m\top}$, $j=1,\ldots,r_{I}^{m}$.

The AJIVE estimator of $V_{J}$ is called the Common Normalized Scores (CNS) matrix, and the estimators of $V_{I}^{m}$ are called the Block Specific Scores (BSS). Similar to PCA scores, the CNS and BSS indicate the relationships between the data vectors. Again an important difference from traditional PCA is that the columns of CNS and BSS are normalized to be orthonormal vectors, hence they are also direction vectors in $\mathbb{R}^{n}$. This is because AJIVE inference is based on studying uncertainty in the angle distribution of the scores. In particular, the joint space, i.e. subspace of $\mathbb{R}^{n}$ spanned by the columns of CNS, defines the estimated joint variation of all data blocks among the experimental subjects. Similarly the subspace of $\mathbb{R}^{n}$ spanned by columns of each of the BSS defines the estimated individual variation.

Estimated modes of variation are obtained from the estimated scores and loadings matrices as outer products of the corresponding columns. The joint and individual estimated modes of variation reflect how the data blocks work together and separately respectively. An important measure of the effectiveness of AJIVE in recovering these modes of variation is the angle between the estimated loadings vector and the true underlying direction.

Using a simple Gaussian simulation, we illustrate how AJIVE identifies modes of variation from multiple data blocks that is the basis of the inference developed here (Figure 1). This toy example is constructed as in (1) across two data blocks to visibly demonstrate the difference between joint and individual modes of variation. In real data, this is not as clear as in this toy example. Figure 1 shows a number of panels. Each is a heat map view of a data matrix whose entries are coded by colors. The color red is coding values that are greater than 1 and blue is coding smaller than -1. White codes zero and in-between values are lighter shades of red or blue. As discussed above, the columns of these matrices represent data vectors or different experimental units. The rows of the matrices represents features sometimes called variables or traits. This toy example has two data blocks: Datablock 1 (shown on the top) and Datablock 2 (bottom). Both have 120 features (rows) and 160 cases (columns). The input data matrices of the toy sample in the far-left panels are the sum of the other three matrices representing joint variation, individual variation, and noise respectively shown in the other panels. The second panels on the left are rank-1 joint matrices. The first 40 rows in the top (Datablock 1) have a similar pattern to the first 40 rows in the bottom, representing important joint features which should be flagged as significant by jackstraw and the remaining rows should be insignificant regarding the joint modes of variation. The third panels on the left are rank-1 individual matrices and the top and bottom panels are very different. They have a relatively weak signal characterized by a minimum value of -0.8 and a maximum of 0.8, as indicated by the paler red and blue colors. The individual jackstraw should label the bottom 80 rows as significant. The far-right panels show simulated independent Gaussian noise with mean 0 and variance 2.

The different joint/individual signal strengths in the toy example in Figure 1 are designed to show how AJIVE recovers signals that are poorly recovered by PCA. This example will also highlight the fact that this richer information improves the statistical power of AJIVE jackstraw tests relative to PCA.

As noted above the jackstraw method [4] for statistical inference on PCA loadings has inspired our approach to AJIVE inference. Technical details are given in Section 2. To demonstrate the improvement of jackstraw with AJIVE over PCA, we applied jackstraw to three different sets of data. 1) AJIVE extracted rank-1 Joint and Individual datasets (Figure 2, top row), 2) PCA modes of variation on Datablocks 1 and 2 separately (Figure 2, middle row) ordered by best correspondence to the row above, and 3) PCA on the concatenated datablocks 1 and 2 (Figure 2, bottom row). Jackstraw significance is indicated by horizontal black lines on the right of each heatmap matrix (using the same color code as Figure 1).

The first major take way from this analysis, was that stronger signal, indicated through deeper red and blue colors, are found in AJIVE compared to PCA (Figure 2, row 1 versus row 2). Concatenation of the datablocks did improve ability of PCA to find joint signal but was not able to resolve individual modes of variation similar to AJIVE (Figure 2, row 3 versus row 1). The second major take away is that the ability for jackstraw to identify the significant joint and individual components. This is observed by the higher number of joint (top 40 rows) and individual (bottom 80 rows) identified as significant by jackstraw in AJIVE compared to PCA.

Each of the $S$ resamples requires one round of the AJIVE algorithm. When the data matrices are large recomputing AJIVE $S$ times could be very computationally intensive. In many situation, the rows are highly correlated and individual features tend not to be particularly influential such as the data in Section 3. Thus permuting a small number $K$ out of a large number $d_{m}$ of rows will lead to minor changes in the AJIVE results. This is especially true in the case of a joint space estimated from multiple data blocks. In that case, one can save a large amount of computational effort by just using the original CNS rather than recomputing AJIVE $S$ times. Note this results in a classical permutation test.

We next demonstrate the effects of not recomputing AJIVE for each of the $S$ replications of jackstraw using TCGA breast cancer data, which will be analyzed in detail in Section 3. As noted above, we have 2 data blocks (data matrices): CNV that includes as features (rows) CNR, and GE that includes as features normalized log of gene expression.

After AJIVE inference, we identify 3 joint components and we proceed to apply jackstraw to find the significant genes/CNRs corresponding to each of the 3 joint components. As an investigation of the Monte Carlo variation, we run the 2 algorithms (full and approximate) 2 times using 2 different seeds, denoted as seed1 and seed2. For convenience of notation, let ‘App’ denote the approximate algorithm and ‘Full’ denote full algorithm. In Table 1, we will report the number of significant CNRs/genes. After a careful comparison, more than 95% of the significant CNRs/genes stay significant if we use different algorithms or random permutations. In each row, the subset corresponding to a smaller number of significant features is always included in the subset corresponding to the larger number. Thus, jackstraw analysis of TCGA breast cancer data is not sensitive to random permutations or algorithms.

While the algorithms gave similar results, the approximate algorithm took about 30 mins and the full algorithm takes about $30\times d_{m}=30\times 20249=607470$ mins for GE without parallel computation. With access to parallel computation, the full algorithm can be much faster. In practice, we recommend making a careful choice between the 2 algorithms based on whether any features are expected to be dominant.

Jackstraw typically involves a very large number of individual hypothesis tests. We recommend three diagnostic graphs (kernal density estimates, jackstraw p-values, and Kolmogorov-Smirnov tests) as useful methods for quickly summarizing the results of the numerous individual hypothesis tests such as those shown in Figure 3 of the toy example in Section 1.1. In the left panel, the black curve is a kernel density estimate visualization of the $F_{null}^{b}$ distribution (on the $log_{10}$ scale) calculated in Step 2 of the algorithm in Section 2. The red (significant) and blue (not significant) dots are the observed F statistics: $F_{j}$ calculated in Step 1. Significance is at the level of 0.05 using the Bonferroni correction. The middle panel shows the sorted jackstraw p-values, computed in Step 3, for all the genes. Out of 120 features, we find 40 significant features in the joint component. The right panel shows the results of a Kolmogorov-Smirnov test (K-S test)[10] for whether the first component overall has useful information. This tests the uniformity of the distribution of these p-values. In this panel, the colored curve is relatively far from the 45-degree line and the K-S p-value is $1.2\times e^{-11}$.

The following discussion demonstrates the kind of information that can be learned from the Common Normalized Scores (CNS). Figure 4 shows a scatter plot matrix view of the CNS, which indicates how the cancer patients relate to each other in terms of GE and CNV. This shows a strong visual separation of the subtypes, that is a consequence of subtypes playing an important role together in GE and CNV, as well as supervision using the third data block. The CNS ($V_{n\times 1}$ in (2)) is a 1038 by 3 matrix because of 3 joint components and 1038 cases/patients. The plots on the diagonal are one-dimensional views of the scores, columns of the CNS matrix. The subtypes are illustrated with subplot kernel density estimates (KDE). Scores are shown using jitter plots (i.e. random heights) to visually separate them in the one-dimensional plots. The overall KDEs are shown in black and the subtype KDEs are shown in subtype colors. The area under the curve of each colored KDE is proportional to its abundance so that the sum of the areas under the colored curves equals 1, the area under the black curve. The off-diagonal plots are scatter plots, i.e. two-dimensional projections of the data onto the planes generated by the corresponding CNS vectors. The first joint component strongly separates the Basal subtypes (red). Joint components 2 and 3 also seem to be strongly related to subtypes even though the separation is less than component 1. The second joint component separates the three other subtypes in the order of LumA (blue), LumB (cyan), and Her2 (magenta). This is sensible because these subtypes are known to be closer together and Her2 is known to be more distinct from LumA. The third joint component separates Her2 (magenta), LumA (blue), and LumB (cyan).

For a deeper understanding of this relationship between subtypes, Figure 5 shows the subtype loading for each joint mode of variation from Figure 4. These loadings indicate the impact of each subtype on the mode of variation. This gives another interpretation of the 3 joint components:

1. 1.

   contrast of Basal vs the rest; contains little Her2 information;

2. 2.

   contrast of (Her2 & LumB) vs LumA; contains little Basal information;

3. 3.

   contrast of LumB vs (Her2 & LumA); contains little Basal information.

While AJIVE analysis indicates distinct joint spaces, we implement the jackstraw method to focus on the driving set of significant features for each joint component. The column App-1 in Table 1 shows the number of jackstraw significant genes and CNRs for each joint component, using the approximate algorithm. For joint component 1, the strong differences between Basal-like subtype and non-Basals in Figures 4 and 5 appear in column App-1 in Table 1 as a large number of significant genes. This may be surprising to statisticians who have approached gene selection via sparsity methods, which has a goal of very few significant genes. However, it clearly underlines the biological complexity of this cancer. The differences within the non-Basal subtypes are weaker, which reflects itself in a smaller number of significant genes.

Our AJIVE version of jackstraw allows us to focus on a set of statistically significant features, giving insight into how GE and CNV work together. In cancer, one important mechanism is copy number variation that can influence expression by increasing or decreasing the numbers of copies per gene. Hence significant copy number variation is expected to lead to significant gene variation. We want to see if genes that are significant in joint space are joint because the copy number is directly regulating expression, however, it’s important to note that there can be many genes located in one CNR. We are interested in the cis effects where the copy number variation directly regulates the gene expression. As shown in the second row of Table 2, almost all significant CNRs contain some significant genes. This is consistent with GE being influenced by its local copy number. However, as shown in the third row of Table 2, not all significant genes in joint space have a dependency on significant CNRs. There are other mechanisms for GE regulation that are still important and show up in the joint space without having a direct dependency on the CNRs. These may be related to a trans effect (non-local regulation) or to presence of strongly correlated genes in the same pathway.

We consider how the jackstraw significance relates to the gene/CNR loadings in Figure 6.Each panel contains colored curves which show the sorted loadings for each gene/CNR. The x-axis represents the full set of the sorted features. The y-axis shows the corresponding loadings. Recall loading vectors have length one, so the sum of the squared entries must be one. In each panel, the curves show joint components 1 (black), 2 (green), and 3 (yellow). Each variable is plotted using jitter plots (i.e. random heights) for joint components 1,2,3. Jackstraw significant genes or CNRs are identified with a red dot and non-significant features with a blue dot. The GE loadings are very similar across the three curves in the left panel, yet the CNV loadings are substantially different in the right panel. In particular, CNV has very few CNRs with very negative loadings. Relatively speaking, the signal and gene expression is spread over more genes, so there are fewer loadings that are much different from zero. There are 20249 genes and only 806 CNRs, thus the jitter plots on the left are much denser than those on the right. In both panels, joint component 1 has more significant features than joint components 2 and 3. This is consistent with the fact that joint component 1 has the most shared variation and is driven mostly by Basal subtype information, which is known to be associated with large genomic changes [14]. Generally, features with large loadings tend to be more significant. Comparing these significant genes with prior work [15] shows the jackstraw analysis has picked up biologically important regions previously identified in these types of cancers. Detailed analysis is shown in Section 4.3 of Xi Yang’s dissertation [16].

Next we investigate feature selection using jackstraw. In particular, we study classification by comparing the results based on each of jackstraw significant features, jackstraw non-significant features and all features. Significance is assessed using the Direction-Projection-Permutation high dimensional hypothesis test (DiProPerm) [17], [18] to quantify subtype separations based on different feature sets. Test results are on the scale of Population Difference Criterion (PDC) as defined in [18]. PDC works on the scale of Gaussian Z-score, so larger than 2 gives statistical significance with larger values reflecting stronger significance. Error bars are included to reflect the simulated variation of the PDC using balanced permutations as explained in Section 3.2 in [18]. Each error bar in Figure 7 shows the confidence interval of the PDC of each test, where a large PDC indicates more significance.

In each panel of Figure 7, we quantify the amount of separation between subtypes and corresponding components indicated in the titles using each of the three feature subsets as the input GE data:

• 1

  All 20249 features (shown as ‘All’ in the left)

• 2

  Jackstraw significant features for each joint component (shown as ‘Sig.’ in the middle)

• 3

  Non-significant features for each joint component (shown as ‘Non.sig’ in the right)

The comparator groups in Figure 7 were taken from the relationships observed in Figure 5. Comp1 corresponds to Basal vs Rest, comp2 corresponds to (Her2&LumB) vs LumA, and comp 3 corresponds to (Her2&LumA) vs LumB. For all 3 panels, the tests using all features have very strong PDCs, indicating strong separations of the subtypes. When we focus on significant features only, the separations become even stronger as indicated by larger PDCs demonstrating the value of focusing on the jackstraw significance. However, the non-significant features have relatively small PDCs, but still retain some signal (PDCs greater than 2). The confidence intervals show that all of these differences are statistically significant relative to the natural permutation variation. In summary, the jackstraw significant feature sets provide the strongest subtype distinction.

As we have seen AJIVE is very effective at finding the joint structure between all three data blocks: GE, CNV, and subtype. Next we discuss the individual information left for each data block. The DiProPerm investigation of the individual GE did not find significant subtype information, PDC of 0.59 (note less than 2). The individual CNV did reveal a significant difference between LumA and LumB, i.e. PDC of 12, which is not joint with GE. This was surprising because class labels were determined by gene expression, which is expected to have information not contained in CNV. In particular, further investigation of the individual space show that there remains informative GE-independent CNRs that can differentiate LumA and LumB. Detailed analysis can be found in Section 4.3 of Xi Yang’s dissertation [16].