Integrative Analysis of Prognosis Data on Multiple Cancer Subtypes using Penalization

In cancer research, profiling studies have been extensively conducted, searching for genes/SNPs associated with prognosis. Cancer is a heterogeneous disease. Examining similarity and difference in the genetic basis of multiple subtypes of the same cancer can lead to better understanding of their connections and distinctions. Classic meta-analysis approaches analyze each subtype separately and then compare analysis results across subtypes. Integrative analysis approaches, in contrast, analyze the raw data on multiple subtypes simultaneously and can outperform meta-analysis. In this study, prognosis data on multiple subtypes of the same cancer are analyzed. An AFT (accelerated failure time) model is adopted to describe survival. The genetic basis of multiple subtypes is described using the heterogeneity model, which allows a gene/SNP to be associated with the prognosis of some subtypes but not the others. A compound penalization approach is developed to conduct gene-level analysis and identify genes that contain important SNPs associated with prognosis. The proposed approach has an intuitive formulation and can be realized using an iterative algorithm. Asymptotic properties are rigorously established. Simulation shows that the proposed approach has satisfactory performance and outperforms meta-analysis using penalization. An NHL (non-Hodgkin lymphoma) prognosis study with SNP measurements is analyzed. Genes associated with the three major subtypes, namely DLBCL, FL, and CLL/SLL, are identified. The proposed approach identifies genes different from alternative analysis and has reasonable prediction performance.


Introduction
Profiling studies have been extensively conducted in cancer research, searching for SNPs (single nucleotide polymorphisms) and genes that are associated with prognosis. Cancer is a heterogeneous disease. Different subtypes of the same cancer usually have different prognosis patterns and different associated genes/SNPs. Consider NHL (non-Hodgkin lymphoma), which is a heterogeneous group of malignancies ranging from very indolent forms to aggressive ones. As discussed in Zhang et al. (2011), different subtypes of NHL are largely different. For example, DLBCL, the largest subtype, is aggressive, whereas FL, the second largest subtype, is indolent. Chromosomal translocations such as t (3,22) are specific to DLBCL, whereas others such as t (14,18)  With data on multiple subtypes, the goal is to identify genes associated with prognosis. For marker identification, we adopt penalization, which has been extensively applied to the analysis of cancer prognosis data with high-dimensional genetic measurements. Single-dataset penalization methods, such as Lasso, SCAD, bridge, MCP and their group counterparts, cannot be directly applied to the analysis of multiple datasets. With multiple datasets, the homogeneity model assumes that, if a functional unit (gene or SNP) is identified, it is concluded as associated with prognosis in all datasets (Liu et al. 2012). An alternative to the homogeneity model is the heterogeneity model, under which a gene or SNP can be associated with prognosis in some datasets but not the others. Under the heterogeneity model, research on penalization methods has been limited (Liu et al. 2012). Compared with the existing studies which analyze gene expression data, the present one has additional complexity. One gene may consist of multiple SNPs, and it is important to allow different effects for SNPs within the same gene. In addition, theoretical properties of the methods developed in Liu et al. (2012) and others have not been established. To the best of our knowledge, the only available method that is tailored to the type of data analyzed in this study is Ma et al. (2012), which adopts thresholding for marker selection. The thresholding method does not have a well-defined objective function. Thus its properties can be very difficult to establish. In addition, it may have more tuning parameters than the penalization method. Compared with the existing studies, another advancement of this study is the analysis of a prognosis data on NHL, which may provide insights into the genetic basis of this deadly disease.
The integrative analysis of data on multiple subtypes of cancer can be challenging. With some cancers, the subtype information may be only partial or even wrong. In addition, the definitions 3 of subtypes are still evolving. For NHL subtypes, we refer to Zhang et al. (2011) and references therein for relevant discussions. When there are a large number of subtypes, the set of subtypes chosen for analysis needs to be jointly determined by the scientific question of interest, quality of data, sample size, evidence from epidemiologic studies and other factors. We acknowledge the importance and difficulty of these issues. In this study, we focus on the development of a new analysis approach and refer to other publications for relevant discussions. The subscript "o" is used to discriminate the original (versus weighted) covariates. For simplicity of notation, assume that different subtypes measure the same set of covariates. In practice, if a covariate is not measured for a specific subtype, its corresponding regression coefficient will be set as zero. In penalization, rescaling is used to accommodate partially matched covariate sets.
For subject i of subtype m, the AFT (accelerated failure time) model assumes that where β m 0 is the intercept, β m ⊆ R d is the regression coefficient, and ǫ m i is the error term. As Compared with alternatives such as the Cox model, the AFT model has a much simpler objective function, as demonstrated in the next section, and hence significantly lower computational cost. Such a property is particularly desirable for high-dimensional data. In addition, it directly describes event times, and its regression coefficients may have more lucid interpretations than those 4 in alternative models. As there is a lack of model diagnostics tools for high-dimensional data, alternative models will not be discussed.

Weighted least squares estimation
LetF m be the Kaplan-Meier estimator of the distribution function F m of T m . Following Stute (1996),F m can be written asF where the ω m i 's are the jumps in the Kaplan-Meier estimator and can be expressed as . Using the weighted centered values, the intercept is zero.
The weighted least squares objective function can be written as This simple form makes computation affordable even with high-dimensional data.
Assume independence between data for the M subtypes. Consider the overall objective function Here we normalize L m by n m so that the analysis is not dominated by large subtypes. When larger subtypes are of more interest, the unnormalized objective function may be considered. 5  To more explicitly describe the data and model settings, heterogeneity model, and our analysis strategy, consider a hypothetical cancer study with three subtypes and eight SNPs representing four genes (Table 1). Gene 1 is associated with the prognosis of all three subtypes; Gene 2 is associated with the first two subtypes but not the third one; Gene 3 is associated with the third subtype only; And gene 4 is not associated with any subtype. In Table 1, we show the regression coefficient matrix whose main characteristics reflect the essence of integrative analysis under the heterogeneity model. Unimportant genes/SNPs not associated with prognosis have no effects and so zero regression coefficients. With penalization approaches including the proposed one, marker identification amounts to identifying the sparsity structure of models. For an important gene/SNP (for example SNP 1 1), its strengths of association with multiple subtypes, which are measured with regression coefficients, may be different for different subtypes. With SNP data, analysis can be conducted at multiple levels, particularly including SNP-level and gene-level. In this study, we conduct gene-level analysis, which complements SNP-level analysis. As the goal is to identify important genes that contain prognosis-associated SNPs, within an important gene, no further selection is conducted. Thus, SNPs within the same gene have the "all in or all out" property.

Heterogeneity model
Such a strategy has been adopted in Ma et al. (2012) and is different from that for SNP-level analysis approaches. 6

Penalized Marker Selection
Penalization is adopted for marker selection. Under the data and model settings described in the previous section, existing penalization approaches are not directly applicable. A new penalization approach is described in Section 3.1. An effective computational algorithm is proposed in Section 3.2. Tuning parameter selection is discussed in Section 3.3. Some practical concerns are discussed in Section 3.4. Asymptotic properties and proofs are established in Appendix.

Penalty function
Assume that the d SNPs belong to J genes. To accommodate the scenario with partially matched gene sets, suppose that M j subtypes (studies) measure gene j. Without loss of generality, assume that gene j is measured in the first M j subtypes. Denote d jk as the number of SNPs in the jth gene and kth subtype with coefficient vector The subscript k is kept to accommodate partially matched SNP sets for the same genes. Then β j = (β ′ j1 , . . . , β ′ jM j ) ′ is the regression coefficient for all SNPs in the jth gene across all subtypes.

A nonzero component ofβ indicates an association between the corresponding gene (SNP) and
subtype's prognosis. Consider the penalty function where λ n > 0 is a data-dependent tuning parameter, c j ∝ M 1−γ j is a constant accommodating partially matched gene sets, || · || is the L 2 norm, and 0 < γ < 1 is the fixed bridge parameter.

7
The above penalty has been designed to tailor the special model characteristics as described in Table 1. In our analysis, genes are the basic functional units. The penalty is the sum of J individual terms, with one for each gene. For a specific gene, two levels of selection need to be conducted. The first is to determine whether it is associated with any subtype at all. This step of selection is achieved using a bridge penalty. For a gene associated with at least one subtype, the second level of selection is to determine which subtype(s) it is associated with. This step of selection is achieved using a Lasso-type penalty. The composition of the bridge-type penalty and the Lasso-type penalty can achieve the desired two-level selection. One gene may contain multiple SNPs. The effect of gene j for subtype k is represented by the vector β jk . Here the penalty is Lasso, which is computationally more expensive, and cannot accommodate the "SNP-within-gene" structure.

Computational algorithm
For subtype m, denote Y m as the vector composed of √ n √ n m Y m ω s, and X m as the matrix composed 3. Repeat Step 2 until convergence.
In numerical study, we use the L 2 norm of the difference between two consecutive estimates less than 0.001 as the convergence criterion. The proposed algorithm always converges, since at each step, the nonnegative objective function decreases. It is noted that as the group bridge type penalty is not convex, the algorithm may converge to a local minimizer depending on the initial value β (0) .
Using the proposed initial value works well in our numerical study.

Tuning parameter selection
The proposed penalty involves two tuning parameters γ and λ n . In the study of bridge type penalties ), the value of γ is usually fixed. Theoretically speaking, different values of γ, as long as in the interval (0, 1), lead to similar asymptotic results. In practice as γ → 1, the bridge type penalty goes to the Lasso type penalty; On the other hand, as γ → 0, it behaves similarly to AIC/BIC type penalties. In our numerical study, we experiment with a few γ values, particularly including 0.5, 0.7 and 0.9. The effect of λ n is similar to that with other penalties. As λ n → ∞, fewer genes/SNPs are identified.
As the function L(β) has a least squares form, we propose using BIC for tuning parameter selection. Particularly, with a fixed γ, the optimal λ n minimizes Here we use the notationβ(λ n ) to emphasize the dependence ofβ on λ n . Motivated by Yuan and Lin (2006), an approximation of the degree of freedom is adopted as Hereβ LS jk is obtained by fitting an AFT model (with least squares estimation) using the jth gene and kth subtype only.

Practical considerations
With practical data, minor allele frequencies in some loci can be low. This may cause an instability problem in the Cholesky decomposition when some eigenvalues of the correlation matrices are too small. In the proposed penalized selection, within-gene-SNP level selection is not of interest.
To reduce the dimensionality within genes and to tackle the colinearity problem, when there is evidence of a lack of stability, we first conduct principal component analysis (PCA) within genes.
Specifically, we choose the number of PCs such that at least 90% of the total variation is explained.
Then the PCs, as opposed to the original SNP measurements, are used in downstream analysis.
Our empirical study suggests that this simple step may ensure that the smallest eigenvalues of the covariance matrices are not too small and that the Cholesky decomposition is stable. Two scenarios under the heterogeneity model are considered. Under the first scenario, all three subtypes share three common susceptibility genes, and each subtype has one subtype-specific susceptibility gene. The "unmatching rate" of susceptibility genes is thus 25%. Under the second scenario, the three subtypes share two common susceptibility genes, and each subtype has two subtype-specific susceptibility genes. The unmatching rate of susceptibility genes is 50%. As a special case of the heterogeneity model, the homogeneity model is also considered, under which all three subtypes have the same susceptibility genes.

Simulation Study
The logarithms of event times are generated from the AFT models with intercept equal to 0.5 and normally distributed random errors. The logarithms of censoring times are generated as uniformly distributed and independent of the event times. The censoring distribution parameters are adjusted so that overall censoring rate is about 30%.
Beyond the proposed approach, simulated data are also analyzed using a meta-analysis approach. Here each subtype is analyzed using the group Lasso (GLasso) approach, where a "group" corresponds to one gene with multiple SNPs. Then the identified gene lists are combined across subtypes. With both approaches, the tuning parameters are chosen using BIC. With the proposed approach, we experiment with γ =0.5, 0.7, and 0.9. Summary statistics on gene identification accuracy based on 100 replicates are shown in Table 2-3 and 5-8 (Appendix).
Simulation suggests that performance of the proposed approach depends on the correlation structure, values of nonzero regression coefficients, and γ value. As correlation gets stronger, in general, more true positives and more false positives are identified. We fail to observe a clear pattern for the dependence (of the performance of proposed approach) on nonzero regression coefficients.
As γ gets larger, also more true positives and more false positives are identified. This observation is reasonable, considering that when γ → 1, bridge penalization becomes close to Lasso, and that Lasso-type penalization tends to over-select. Under almost all simulated scenarios, the proposed approach identifies more true positives than GLasso. For example in Table 2  Analysis results using the proposed approach are shown in Table 4 and 9-11 (Appendix). In particular, Table 4 contains the L 2 norms of the identified genes, whereas Table 9 Table 4 shows that only gene IL10 for DLBCL has a low occurrence index (0.21). All other observed occurrence indexes are high, suggesting relatively satisfactory stability. Prediction performance is also evaluated using a random sampling approach. In particular, genes are identified and models are constructed using 3/4 of randomly sampled subjects. Then prediction is made for the rest 1/4 subjects. Based on the predicted X m′βm , subjects are separated into two risk groups.
The logrank statistic is computed to compare the survival risk of the two groups. This process is repeated 100 times, and the mean logrank statistic is computed as 7.1 (p-value 0.0077), suggesting satisfactory prediction.
For comparison, we also analyze each subtype separately using GLasso (results shown in Ta shows that the identified genes also have satisfactory stability. Prediction evaluation generates a logrank statistic of 0.2 (p-value 0.65), which is considerably smaller than that using the proposed approach.

Discussion
In this study, with prognosis data on multiple subtypes of the same cancer, we develop a penalization approach which can conduct integrative analysis, identify important genes that contain SNPs associated with multiple subtypes, and allow for subtype-specific susceptibility genes. The proposed approach can be realized using an effective iterative algorithm. Under mild conditions, it has the much desired consistency properties. Simulation shows that the proposed approach outperforms penalization-based meta-analysis, with more true positives and fewer false positives. In the analysis of NHL prognosis data, it identifies multiple genes shared by two or three subtypes as well as subtype-specific genes. The shared genes have important biological implications. The proposed approach also leads to significantly better prediction performance.
To avoid confusion, in our description we focus on the scenario with multiple subtypes of the same cancer and the "SNP-within-gene" structure. The proposed approach is directly applicable to the analysis of multiple types of cancers and "gene-within-cluster (pathway)" and other structures.
In addition, with minor modifications, analysis of prognosis data under other models and analysis of diagnosis data can be conducted. The proposed penalty is built on bridge-type penalties. We conjecture that it is possible to build on other penalties such as MCP. Our limited investigation shows that under the present setup, the proposed penalty may have the lowest computational cost.
In data analysis, our preliminary search shows that the common genes shared by multiple subtypes have important implications. However, because of the following limitations, the analysis results should be interpreted with caution. First, the sample size is still limited. Second, the NHL study takes a candidate gene approach. It is possible that important genes have been missed in the profiling stage. Third, the proposed evaluation is cross-validation based. Although it can compare different approaches on the same ground, it does not use completely independent data. More, independent studies are needed to fully comprehend the data analysis results.    Table 3: Simulation under the heterogeneity model: unmatching rate=50% and nonzero regression coefficients under case 1. In each cell, the first row is number of true positives (standard deviation), and the second row is model size (standard deviation).

Statistical properties and proofs
In this section, we use the same notations as in Section 2 and 3.1. The PCA step described in Section 3.4 is optional and mainly for practical consideration. If PCA is actually conducted, the d jk values may be smaller.
Let B 1 and B 2 be the index sets of genes with nonzero and zero norms of regression coefficients, respectively. Let β 0 be the true parameter vector of β, q m be the size of the set {j : β 0jm = 0}, and q = M m=1 q m (that is, the total number of associations between genes and prognosis across 1 all subtypes). Without loss of generality, suppose that β 0jk = 0 for all j = 1, . . . , J 1 and some k = 1, . . . , M j . With X as the full design matrix (defined in Section 3.2), denote X B1 as the design matrix corresponding to B 1 . Define Σ n = n −1 X ′ X, Σ 1n = n −1 X ′ B1 X B1 . Further for any index set A, let X A denote the corresponding design matrix, and Σ A = n −1 X ′ A X A . Let We make the following assumptions: ( That is, there exist constants 0 < c * < c * < ∞, such that for q * = (3 + 4C)q and C = c * /c * , with probability converging to 1, c * ≤ ν ′ Σ A ν ν 2 ≤ c * , ∀A with |A| = q * and ν ∈ R q * . In addition, the SRC condition holds for the design matrix of each subtype separately with rank q m * for subtype m. (A4) and η n = O(1).
The above assumptions are in parallel with those in Huang and Ma (2010). Further complications are introduced to accommodate the multi-datatsets setting, heterogeneity across subtypes, and "SNP-within-gene" structure. Main properties of the penalized estimate can be summarized as follows. The above result establishes that the number of identified genes is a finite multiply of the true number of associated genes, which is assumed to be finite in (A1). In addition, if log(d)/n → 0, the estimate is L 2 estimation consistent. The selection and estimation results can be further strengthened as follows.
2. Suppose that {B 1 , β 0B 1 } are fixed unknown. In the AFT models, the subgaussian condi- Then, in distribution, with u jk corresponding to the component of gene j and subtype k.
Proof of Theorem 7.1. Part (i) follows from the proof of Theorem 1 of Zhang and Huang (2008).
One difference is that here a subgaussian assumption is made, which is weaker than the normality assumption in Zhang and Huang (2008). Since subgaussian random variables have the same tail behaviors as normal random variables, the argument of Zhang and Huang (2008) goes through. In addition, the SRC condition is imposed at a group level, as opposed to individual covariate level, to accommodate the "SNP-within-gene" structure.
Following from (16) and (17), we have Let ζ B be the projection of ζ to the span of X B , i.e., Therefore, by the Cauchy-Schwarz inequality, Combining (18) and (19), By the SRC condition (A3), It follows that We have Therefore, The result follows from (21) and (22).
Proof of Theorem 7.2. (i) We prove the first part of Theorem 7.2.

Defineβ byβ
The Karush-Kuhn-Tucker condition for (9) implies that where X jl is the sub-matrix corresponding to the jth gene of the lth subtype.
By the definition ofβ, we have Thus, by SRC condition (A3), with n −1 X(β −β) 2 ≤ c * β −β 2 , we have which implies, by Theorem 7.1 (ii), We still need to find a lower bound of J Since c j ≥ 1 by assumption, If β B 2 > 0, the combination of (24) and (25) yields Since λ n (logd/n) γ/2−1 → ∞ by assumption, this implies that (ii) We now prove the second part. Let h n = n −1/2 define with 0 being the zero vector of dimension |B 2 |. By (i), the following holds with large probability: The function V 1n (u), u ∈ R d 1 , can be written as For the first term, we have For the second term,   In each cell, the first row is number of true positives (standard deviation), and the second row is model size (standard deviation).