Supervised group Lasso with applications to microarray data analysis

Colon data

In this dataset, expression levels of 40 tumor and 22 normal colon tissues for 6500 human genes are measured using the Affymetrix gene chips. A selection of 2000 genes with the highest minimal intensity across the samples has been made by [2], and these data are publicly available at [21]. The colon data have been analyzed in several previous studies using other statistical approaches, see for example [3, 22, 23].

Nodal data

This dataset was first presented by [24, 25]. It includes expression values of 7129 genes of 49 breast tumor samples. The expression data were obtained using the Affymetrix gene chip technology and are available at [26]. The response describes the lymph node (LN) status, which is an indicator for the metastatic spread of the tumor, a very important risk factor for the disease outcome. Among the 49 samples, 25 are positive (LN+) and 24 are negative (LN-). We threshold the raw data with a floor of 100 and a ceiling of 16000. Genes with max(expression)/min(expression) < 10 and/or max(expression) – min(expression) < 1000 are also excluded [1]. 3332 (46.7%) genes pass the first step screening. A base 2 logarithmic transformation is then applied. The Nodal data have also been studied by [22].

Although there is no limitation on the number of genes that can be used in the proposed approach, we first identify 500 genes for each dataset based on marginal significance to gain further stability as in [23]. Compute the sample standard errors of the d biomarkers se₍₁₎,..., se_(d)and denote their median as med.se. Compute the adjusted standard errors as 0.5(se₍₁₎ + med.se),..., 0.5(se_(d)+ med.se). Then the genes are ranked based on the t-statistics computed with the adjusted standard errors. The 500 genes with the largest absolute values of the adjusted t-statistics are used for classification. The adjusted t-statistic is similar to a simple shrinkage method discussed in [27].

For the Colon and Nodal data, clusters are constructed using the K-means approach and the Gap statistic is used to select the optimal number of clusters. We show in Figure 1 the Gap statistic as a function of the number of clusters. 9 clusters are constructed for the Colon data and 20 clusters are constructed for the Nodal data. Details of the clustering information are available upon request. With the generated clusters, we apply the proposed SGLasso approach. Tuning parameters are chosen using 3-fold cross validation. Summary model features are shown in Table 1. For the Colon data, 22 genes are present in the final model, representing 8 clusters. For the Nodal data, 66 genes are selected, representing 17 clusters. We list the identified genes in Tables 2 and 3.

For the Colon data, gene Has.1039 has also been identified to be associated with Colon cancer in [28, 29]. Gene Hsa.42949 is estrogen sulfotransferase. Research show that certain compounds, such as soy, have protective effect for colon cancer. The protective role of these compounds could be due to an ability to inhibit competitively the activation of promutagenic estrogen metabolites into carcinogens by estrogen sulfotransferases. The official symbol of gene Hsa.1454 is CSNK1E. Studies have revealed a negative regulatory function of CK1 in the Wnt signaling pathway, where CK1 acts as a negative regulator of the LEF-1/beta-catenin transcription complex, thereby protecting cells from development of cancer. Gene Hsa.8214 has official symbol DCGR6. It has been shown to be associated with mammary cancer and tumor cell proliferation in general. Gene Hsa.462 (official symbol SERPINC1) has been shown to be related to cancer cell proliferation. Gene Hsa.627 is also identified as a Colon cancer biomarker in [28]. Gene Hsa.696 has official symbol BTN1A1. RT-PCR analysis has revealed strongest expression of BTNL3 in small intestine, colon, testis, and leukocytes. Tristetraproline (gene Hsa.1682) has been reported to negatively regulate tumor necrosis factor alpha (TNF-alpha) production by binding the AU-rich element within the 3' noncoding sequences of TNF-alpha mRNA. Gene Hsa.3016 is S100P protein. 100P is expressed in human cancers, including breast, colon, prostate, and lung. In colon cancer cell lines, its expression level was correlated with resistance to chemotherapy.

For the Nodal data, gene U27185_at has official Symbol RARRES1. Also known as TIG1, the expression of this gene is upregulated by tazarotene as well as by retinoic acid receptors. Silencing of TIG1 promoter by hypermethylation is common in human cancers and may contribute to the loss of retinoic acid responsiveness in some neoplastic cells. The role of the matrilins (gene U69263) in tumorigenesis has not been studied. However, a related family of proteins (fibulins) has been implicated in cancer. Increased fibulin expression is seen in breast cancer, lung adenocarcinoma, colon cancer, and other solid tumors, suggesting that these proteins might play a role in tumor formation or progression. Gene U07223 is a member of the chimerin family and encodes a protein with a phorbol-ester/DAG-type zinc finger, a Rho-GAP domain and an SH2 domain. Decreased expression of this gene is associated with high-grade gliomas and breast tumors, and increased expression of this gene is associated with lymphomas. Findings suggest that CEACAM1 (gene X16354) participates in immune regulation in physiological conditions and in pathological conditions, such as inflammation, autoimmune disease, and cancer. Gene D87071_at is a confirmed breast cancer biomarker. Ras (gene J00277) and c-Myc play important roles in the up-regulation of nucleophosmin/B23 during proliferation of cells associated with a high degree of malignancy, thus outlining a signaling cascade involving these factors in the cancer cells. GNA11 (gene M69013) is involved in signaling of gonadotropin-releasing hormone receptor, which negatively regulates cell growth.

Down-regulation is suggested to be involved in human breast cancers. Gene D59532 encodes a member of the C-type lectin/C-type lectin-like domain (CTL/CTLD) superfamily. Members of this family share a common protein fold and have diverse functions, such as cell adhesion, cell-cell signaling, glycoprotein turnover, and roles in inflammation and immune response. Gene X53587 is human mRNA for integrin beta 4. Colonization of the lungs by human breast cancer cells is correlated with cell surface expression of the alpha(6)beta(4) integrin and adhesion to human CLCA2 (hCLCA2), Tumor cell adhesion to endothelial hCLCA2 is mediated by the beta(4) integrin, establishing for the first time a cell-cell adhesion property for this integrin that involves an entirely new adhesion partner. This adhesion is augmented by an increased surface expression of the alpha(6) beta(4) integrin in breast cancer cells selected in vivo for enhanced lung colonization but abolished by the specific cleavage of the beta(4) integrin with matrilysin. Wnt-5a (gene L20861) has been shown to influence the metastatic behavior of human breast cancer cells, and the loss of Wnt-5a expression is associated with metastatic disease. NFAT1, a transcription factor connected with breast cancer metastasis, is activated by Wnt-5a through a Ca2+ signaling pathway in human breast epithelial cells. Endogenous RXR beta (gene M84820) contributes to ERE binding activity in nuclear extracts of the human breast cancer cell line MCF-7. Detailed microscopic analysis of the morphology of MCF7 breast cancer cells lacking CtBPs (gene U37408) reveals an increase in the number of cells containing abnormal micronucleated cells and dividing cells with lagging chromosomes, indicative of aberrant mitotic chromosomal segregation. Methyl-CpG-binding domain protein-2 (gene X99687) mediates transcriptional repression associated with hypermethylated GSTP1 CpG islands in MCF-7 breast cancer cells. The nmt55/p54nrb protein (gene U02493) is post-transcriptionally regulated in human breast tumors leading to reduced expression in ER- tumors and the expression of an amino terminal altered isoform in a subset of ER+ tumors. Expression of elastin (gene U77846) in breast carcinoma cells has been demonstrated by immunohistochemistry and in situ hybridization. Cytochrome P450 1B1 (CYP1B1, gene X07618) is active in the metabolism of estrogens to reactive catechols and of different procarcinogens. The CYP1B1 gene polymorphisms do not influence breast cancer risk overall but may modify the risk after long-term menopausal hormone use. Genetic deficiency of DPYD enzyme (gene U09178) results in an error in pyrimidine metabolism associated with thymine-uraciluria and an increased risk of toxicity in cancer patients receiving 5-flourouracil chemotherapy.

For the alternative approaches, we also compute the LOO cross validation errors. Tuning parameters when presented are also chosen via 3-fold cross validation. Comparison results are shown in Table 1.

We can see from Table 1 that the SGLasso is capable of feature selection at both the cluster level and the within cluster gene level. The number of genes selected is much less than its counterpart from the simple clustering approach and GLasso. The Lasso is also capable of selecting a small number of genes. Especially we note that the number of genes selected by Lasso is smaller than by SGLasso. For simple data sets such as the Colon data, the Lasso prediction error is the same as the SGLasso. However for data sets that are more difficult to classify (Nodal), the SGLasso prediction error is much smaller. Both data sets have also been analyzed by other approaches. For the Colon data, ROC based approach has prediction error 0.14 [23]; LogitBoost has classification errors 0.145, 0.194 and 0.161 [22]; and classification tree has classification error 0.145 [22]. We note that since different sets of genes are used in those studies, Table 1 only provides rough comparisons. For the Nodal data, in [22], LogitBoost yields prediction error 0.184, 0.265 and 0.184, while classification tree has prediction error 0.204 and 1-nearest neighbor has prediction error 0.367.

Follicular lymphoma data

Follicular lymphoma is the second most common form of non-Hodgkin's lymphoma, accounting for about 22 percent of all cases. A study was conducted to determine whether the survival probability of patients with follicular lymphoma can be predicted by the gene-expression profiles of the tumors at diagnosis [5]. Fresh-frozen tumor-biopsy specimens and clinical data from 191 untreated patients who had received a diagnosis of follicular lymphoma between 1974 and 2001 were obtained. The median age at diagnosis was 51 years (range 23 to 81), and the median follow up time was 6.6 years (range less than 1.0 to 28.2). The median follow up time among patients alive at last follow up was 8.1 years. Eight records with missing survival information are excluded from the downstream analysis. Detailed experimental protocol can be found in [5].

Affymetrix U133A and U133B microarray genechips were used to measure gene expression levels from RNA samples. A log2 transformation was applied to the Affymetrix measurements. We first filter the 44928 gene measurements with the following criteria: (1) the max expression value of each gene across 191 samples must be greater than 9.186 (the median of the maximums of all probes). (2) the max-min should be greater than 3.874 (the median of the max-min of all probes). (3) Compute correlation coefficients of the uncensored survival times with gene expressions. Select the genes whose correlation with survival time is greater than 0.2. There are 729 genes that pass this screening process. We normalize genes across samples to have mean 0 and variance 1.

Mantel cell lymphoma data

[4] reported a study using microarray expression analysis of mantle cell lymphoma (MCL). The primary goal of this study was to discover genes that have good predictive power of patient's survival risk. Among 101 untreated patients with no history of previous lymphoma included in this study, 92 were classified as having MCL, based on established morphologic and immunophenotypic criteria. Survival times of 64 patients were available and other 28 patients were censored. The median survival time was 2.8 years (range 0.02 to 14.05 years). Lymphochip DNA microarrays [15] were used to quantify mRNA expression in the lymphoma samples from the 92 patients. The gene expression data that contain expression values of 8810 cDNA elements are available at [30].

We pre-process the data as follows to exclude noises and gain further stability: (1) Compute the variances of all gene expressions; (2) Compute correlation coefficients of the uncensored survival times with gene expressions; and (3) Select the genes with variances larger than the first quartile and with correlation coefficients larger than 0.25. 834 out of 8810 genes (16.5%) pass the above initial screening. We standardize these genes to have zero mean and unit variance.

For the Follicular data, gene 23098_x_a is associated with tumor protein p53. In transfected cells and KSHV-infected B lymphoma cells, KSHV-encoded latency-associated nuclear antigen (LANA) expression stimulates degradation of tumor suppressors von Hippel-Lindau and p53. In a recent case study with a Japanese girl who had EEC3 and developed diffuse large B-cell type non-Hodgkin lymphoma, researchers identified heterozygosity for a 1079A-G transition in exon 8 of the TP73L gene, resulting in a germline asp312-to-gly (D312G) mutation. Gene 223333_s_a is a member of the angiopoietin/angiopoietin-like gene family and encodes a glycosylated, secreted protein with a fibrinogen C-terminal domain. The encoded protein may play a role in several cancers and it also has been shown to prevent the metastatic process by inhibiting vascular activity as well as tumor cell motility and invasiveness. Gene 224357_s_a encodes a member of the membrane-spanning 4A gene family. Members of this nascent protein family are characterized by common structural features and similar intron/exon splice boundaries and display unique expression patterns among hematopoietic cells and nonlymphoid tissues. Chemokines (genes 204470_at, 205114_s_a) are a group of small (approximately 8 to 14 kD), mostly basic, structurally related molecules that regulate cell trafficking of various types of leukocytes through interactions with a subset of 7-transmembrane, G protein-coupled receptors. Chemokines also play fundamental roles in the development, homeostasis, and function of the immune system, and they have effects on cells of the central nervous system as well as on endothelial cells involved in angiogenesis or angiostasis. Serpin Al (gene 211429_s_a) has an invasion-promoting effect in anaplastic large cell lymphoma. Gene 216950_s_a has official symbol FCGR1A. Findings showed that both Fcgamma RIA and FcgammaRIIA mediated enhanced dengue virus immune complex infectivity but that FcgammaRIIA appeared to do so far more effectively. TRPM4-mediated (gene 219360_s_a) depolarization modulates Ca2+ oscillations, with downstream effects on cytokine production in T lymphocytes. The protein encoded by gene 210973_s_a is a member of the fibroblast growth factor receptor (FGFR) family, where amino acid sequence is highly conserved between members and throughout evolution. Chromosomal aberrations involving this gene are associated with stem cell myeloproliferative disorder and stem cell leukemia lymphoma syndrome. FGFR-1 is expressed in early hematopoietic/endothelial precursor cells, as well as in a subpool of endothelial cells in tumor vessels. Gene 227697_at encodes a member of the STAT-induced STAT inhibitor (SSI), also known as suppressor of cytokine signaling (SOCS), family. Over expression of suppressor of cytokine signaling 3 is associated with anaplastic large cell lymphoma.

Genes identified in the MCL study also have sound biological basis. When the positive cells are treated with phosphatidylinositol-specific phospholipase C (gene Hs.522568), a significant decrease in both stain intensity and percentage of positive cells is demonstrated by immunofluorescence. The protein encoded by gene Hs. 120949 is the fourth major glycoprotein of the platelet surface and serves as a receptor for thrombospondin in platelets and various cell lines. Since thrombospondins are widely distributed proteins involved in a variety of adhesive processes, this protein may have important functions as a cell adhesion molecule. Mutations in this gene cause platelet glycoprotein deficiency. The protein encoded by gene Hs.84113 belongs to the dual specificity protein phosphatase family. It was identified as a cyclin-dependent kinase inhibitor, and has been shown to interact with, and dephosphorylate CDK2 kinase, thus prevent the activation of CDK2 kinase. This gene was reported to be deleted, mutated, or overexpressed in several kinds of cancers. Studies show TOP2A (gene Hs.156346) is a proliferation marker, indicator of drug sensitivity, and prognostic factor in mantle cell lymphoma. This gene encodes a DNA topoisomerase, an enzyme that controls and alters the topologic states of DNA during transcription. The gene encoding this enzyme functions as the target for several anticancer agents and a variety of mutations in this gene have been associated with the development of drug resistance. Reduced activity of this enzyme may also play a role in ataxia-telangiectasia. Gene Hs.497741 encodes a protein that associates with the centromere-kinetochore complex. Autoantibodies against this protein have been found in patients with cancer or graft versus host disease. DNA Pol theta (gene Hs.241517) has a specialized function in lymphocytes and in tumor progression. Gene Hs.298990 encodes tumor suppressor proteins. The protein encoded by gene Hs. 105956 catalyzes the transfer of galactose to lactosylceramide to form globotriaosylceramide, which has been identified as the P(k) antigen of the P blood group system. PLG (gene Hs.368912) has the potential to simultaneously regulate calcium signaling pathways and regulate pHi via an association with NHE3 linked to DPP IV, necessary for tumor cell proliferation and invasiveness.

In the leave-one-out (LOO) based evaluation, denote $\widehat{\beta }$(-i) as the LOO estimate of β based on the reduced data set with the i ^th subject removed. We then compute the predicted risk score $\widehat{\beta }$(-i)${Z^{\prime }}_i$ for the i ^th subject. Since the Cox model is a special form of the transformation model, the uncensored survival time depends on β Z' via a generalized linear model. So the prediction evaluation can be based on comparing the survival functions of groups composed of different range of β Z'. A simple approach is to first dichotomize the predicted risk scores at the median to create two risk groups with equal sizes. We then compare the survival functions of the two generated risk groups. A significant difference (measured by the logrank test statistic with degree of freedom 1) indicates satisfactory prediction performance.

We show the prediction comparison in Table 1. We can again see that the models obtained under SGLasso are much smaller than those from the simple cluster approach and GLasso. However the SGLasso models are larger than their Lasso counterparts. For both data sets, the proposed SGLasso has the largest logrank statistics, indicating the best prediction performance. For the Follicular data, only Lasso and SGLasso have logrank statistics with p-value less than 0.05. The GLasso and simple approaches cannot properly predict survival based expressions. For the MCL data, all approaches have logrank statistics with p-value less than 0.05, with the largest logrank statistic from the SGLasso.

Algorithm I: L₁boosting Lasso

The β ⁱ at convergence is the Lasso estimate. We conclude convergence if the absolute value of ${\phi }^{k^{*}}$ (β) computed in step 3 is less than a pre-defined criteria, and/or if R _n (β) is larger than a pre-defined threshold. Alternative algorithm can be LARS based. Since it is not the focus of this study, we omit discussions of other computational algorithms.

For the GLasso, a LARS based approach is proposed in [19]. With high dimensional cases, a computationally more affordable approach is proposed in [39]. This approach shares the same spirit as the L₁ boosting Lasso and they are both special cases of the gradient-based constraint maximization discussed in [40]. This boosting based algorithm can be summarized into the following iterations.

Algorithm II: boosting group Lasso