Abstract
Feng Chen1,*, Peng Peng1,*, Yi Zhou1, Zhen-Yu Yang1, Hai-Quan Zhang1, Xiang-Sheng Ao1, Da-Quan Zhou1 and Chun-Xiang Xiang2
1Department of Neurosurgery, Xiangyang Central Hospital, Hubei University of Arts and Science, Xiangyang, Hubei 441021, P. R. China
2Department of Pathology, Xiangyang Central Hospital, Hubei University of Arts and Science, Xiangyang, Hubei 441021, P. R. China
*These authors have contributed equally to this work
Correspondence to:
Da-Quan Zhou, email: [email protected]
Chun-Xiang Xiang, email: [email protected]
Keywords: glioblastoma, prognosis, biomarker, DNA microarray, LASSO
Received: February 04, 2017 Accepted: July 24, 2017 Published: August 10, 2017
ABSTRACT
Gliomas are the most common lethal brain tumours and remain great heterogeneity in terms of histopathology and clinical outcomes. Among them, glioblastomas are the most aggressive tumours that lead to a median of less than one-year survival in patients. Despite the little improvement of in diagnosis and treatments for last decades, there is an urgent need for prognostic markers to distinguish high- and low-risk patients before treatment. Here, we generated a list of genes associated with glioblastoma progressions and then performed a comprehensive statistical modelling strategy to derive a 10-gene (GLO10) score from genome wide expression profiles of a large glioblastoma cohort (n=844). Our study demonstrated that the GLO10 score could successfully distinguish high- and low-risk patients with glioblastomas regardless their traditional pathological factors. Validated in four independent cohorts, the utility of GLO10 score could provide clinicians a robust prognostic prediction tool to assess risk levels upfront treatments.
INTRODUCTION
Glioma is among the most common brain tumour type. With median survivals of less than a year for patients, glioblastoma is considered as a vital cause of cancer mortality in both adults and children despite aggressive surgery, chemotherapy and radiation [1, 2]. By the year of 2007, according to the World Health Organization (WHO), gliomas were assessed based on histopathological features and categorized into astrocytomas, oligodendrogliomas, mixed oligoastrocytomas, and ependymomas [3]. However, limitations of this grading system were addressed, such as high rate of inter- and intra-observer variability [4, 5]. Therefore, in 2016, molecular signature was introduced in diagnosing tumours due to the explosion of genomic information during the last decade [6].
DNA microarrays were introduced two decades ago. As this technology provides more comprehensive and objective information than traditional microscopic morphology, it has revolutionized cancer research. To date, there are molecular based diagnostic tests using DNA microarrays. Some have been incorporated in clinical practice guidelines, including Mamma Print and Oncotye Dx [7]. A number of DNA microarray based studies have identified prognostically distinct molecular subtypes of gliomas [8–12]. These approaches were based on unsupervised hierarchical or k-means clustering of genes. As reported in these studies, patients belonging to various subgroups (clusters) showed significantly different prognostic outcomes independent to some known clinical factors, such as age and grade. However, the clustering-based prognostic signatures remain challengeable in clinical practice as the results were dramatically affected by the selected genes. Previous studies remain inconsistency in the prognostic gene candidates, which increases the difficulties in making such methods practical applicable.
Here, we generated a 10-gene (GLO10) score for the prediction of glioblastoma overall survival from DNA microarray expression datasets. We trained the prognostic models in a large cohort (n=470) and demonstrated that the GLO10 score successfully distinguished high- and low-risk groups in other three independent cohorts comprising patients with glioblastomas (n=374). Our results unveiled novel glioma prognostic biomarkers that could be easily applied with great potential in producing robust results in clinical practices.
RESULTS
Identification of glioblastoma associated genes
We hypothesized that the prognostic biomarkers were associated with tumour genesis and progression. In this case, the glioblastoma associated genes were identified by comparing gene expression profiles of tumour tissues and normal brain samples. The data was retrieved from one previously published large scale study (n= 256), Repository of Molecular Brain Neoplasia Data (Rembrandt)[13]. After data preprocessing and differential expression analyses, a total of 723 probe sets corresponding to 552 genes were identified as significantly differentially expressed genes (DEGs) in glioblastomas. Among the DEGs, 137 were up-regulated and 415 were down-regulated (Figure 1). The gene functional analyses revealed that the glioblastoma associated genes were significantly enriched in the Gene Ontology Cellular Components including postsynapse (P=1.36×10-27), axon(P=8.14×10-23), ion channel complex(P=4.75×10-16), asymmetric synapse(P=1.90×10-15), neuron to neuron synapse(P=2.33×10-15), exocytic vesicle membrane(P=5.92×10-11), extracellular matrix component (P=3.34×10-06) and cytoplasmic region (P=5.43×10-03); Biological Processes including modulation of synaptic transmission(P=3.94×10-15), regulation of neurotransmitter levels(P=1.94×10-11), regulation of neuron projection development(P=3.24×10-10), regulation of cell-substrate adhesion(P=5.05×10-04), positive regulation of fibroblast migration(P=5.49×10-03); Molecular Functions including gated channel activity(P=1.143×10-13), GABA receptor activity(P=8.95×10-09), channel regulator activity(P=1.45×10-07), calmodulin binding(P=7.21×10-07), ligand-gated anion channel activity(P=1.02×10-05), ion channel regulator activity(P=2.50×10-05), ligand-gated cation channel activity(P=1.27×10-03) and structural constituent of myelin sheath(P=1.54×10-03)(Figure 2). The pathway enrichment analysis (Figure 3) revealed that the identified glioblastoma associated genes were enriched in various pathways that involved in glioblastoma progression, for example GABA receptor activation [14], Ion channel transport [15, 16] and Extracellular matrix organization [17–19].
Figure 1: Volcano plot of differentially expressed genes in glioblastomas. The red dots indicate those were considered as significant up- and down- regulated genes.
Figure 2: Gene set enrichment analysis of glioblastoma associated genes using Gene Ontology (Cellular Component, Biological Process and Molecular Function).
Figure 3: Gene set enrichment analysis of glioblastoma associated genes using Reactome Pathway database.
Univariate prognostic analysis by Cox proportional hazards models
To investigate the prognostic potentiality, we integrated a larger multi-institutional data set of 470 glioblastoma patients from the TCGA cohort (Discovery set). Previously identified glioblastoma associated genes were fitted into Cox promotional hazards models. Based on the overall survival association significance assessed by log rank test, these probe sets were ranked and filtered at the cutoff of P < 0.001. A total of 28 probe sets were selected for statistical modelling.
LASSSO statistical modelling to derive GLO10 score
A statistical regression algorithm based on the least absolute shrinkage and selection operator (LASSO) was applied to relate the 28 probe sets to the patient survival in the discovery set. We have optimized a 10-gene (Table 2) signature (GLO10) as the weighted sum of expression levels of the 10 genes for each patient. High GLO10 scores were significantly associated with shorter overall survival times. After the parameter optimization, we have categorized the patients into high-, or low-risk groups based on the GLO10 score threshold of 3.58. The high- and low-risk groups showed significant differences in terms of overall survival (Figure 4, Table 1).
Figure 4: Evaluation of the association of GLO10 score with glioblastoma survival. Left: Kaplan-Mert plot of high- and low-risk groups, P value was calculated using the log-rank test. Middle: Comparison of gender ratios in high- and low-groups, P value was calculated using Pearson’s chi-squared test. Right: Comparison of the distributions of age in high- and low-risk groups, P value was calculated using Wilcoxon rank-sum test.
Table 1: Clinical characteristics of discovery and validation sets
Discovery Set: TCGA | ||||
|---|---|---|---|---|
Characteristics | Overall | High risk group | Low risk group | P-value |
(n=470) | (n=193) | (n=277) | ||
Median overall survival[months](range) | 10.37 (0.10-127.55) | 10.03 (0.10-127.55) | 10.58 (0.10-108.81) | 0.020§ |
Female [n (%)] | 181 (38.51%) | 77 (39.90%) | 104 (37.54%) | 0.68‡ |
Median age[years](range) | 59.00 (5.00-89.00) | 58.00 (10.00-88.00) | 59 (5.00-89.00) | 0.20† |
Validation Set 1: GSE13041 | ||||
Characteristics | Overall | High risk group | Low risk group | P-value |
(n=160) | (n=24) | (n=136) | ||
Median overall survival[months](range) | 12.84 (0.23-110.22) | 9.04 (1.74-72.02) | 13.58 (0.23-110.22) | 0.044§ |
Female [n (%)] | 63 (39.38%) | 7 (29.17%) | 56 (41.18%) | 0.38‡ |
Median age[years](range) | 51.50 (18.00-86.00) | 57.50 (40.00-85.00) | 49.50 (18.00-86.00) | 0.049† |
Validation Set 2: GSE16011 | ||||
Characteristics | Overall | High risk group | Low risk group | P-value |
(n=155) | (n=89) | (n=66) | ||
Median overall survival[months](range) | 8.76 (0.24-150.72) | 7.08 (0.24-66.72) | 14.16 (0.48-150.72) | <0.0001§ |
Female [n (%)] | 50 (32.26%) | 27 (30.34%) | 23 (34.85%) | 0.67‡ |
Median age[years](range) | 55.00 (14.00-80.00) | 58.00 (14.00-79.00) | 46.50 (15.00-80.00) | 1.45×10-6† |
Validation Set 3: GSE83294 | ||||
Characteristics | Overall | High risk group | Low risk group | P-value |
(n=59) | (n=5) | (n=54) | ||
Median overall survival[months](range) | 7.80 (0.23-41.00) | 5.03 (1.74-7.80) | 9.29 (0.23-41.00) | 0.014§ |
Female [n (%)] | 32 (54.24%) | 1 (20.00%) | 31 (57.41%) | 0.26‡ |
Median age[years](range) | 47.00 (1.00-82.00) | 56.00 (39.00-66.00) | 44.50 (1.00-82.00) | 0.35† |
‡ P value calculated using the Pearson’s chi-squared test.
† P value calculated using the Wilcoxon rank-sum test.
§ P value calculated using the log-rank test.
To assess the robustness of the GLO10 score, we then evaluated its performance in other three independent glioblastoma cohorts with GEO accessions GSE13041 (validation set 1, n=160), GSE16011 (validation set 2, n=155) and GSE83294 (validation set 3, n=59). Patients with GLO10 scores higher than 3.58 were assigned into the high risk group, while those with GLO10 scores lower than the threshold were considered as low risk. As shown in Figure 4, there are significant differences between high-risk and low-risk groups (validation set 1: P = 0.044, validation set 2: P < 0.0001, validation set 3: P=0.014). Surprisingly, although high-risk patients tended to have higher age, the gender ratio showed insignificant different patterns (Figure 4, Table 1). This indicated that the GLO10 score retained prognostic power independent of other traditional factors, such as age and gender. The dominance of GLO10 score in the prognosis prediction revealed the importance of the 10 genes (Table 2) in glioblastoma diagnosis and treatment researches.
Table 2: List of gens in the GLO10 score
Probeset ID§ | Symbol | Gene ID* | Chromosomal location | Description |
|---|---|---|---|---|
221898_at | PDPN | 10630 | 1p36.21 | podoplanin |
202133_at | WWTR1 | 25937 | 3q25.1 | WW domain containing transcription regulator 1 |
203706_s_at | FZD7 | 8324 | 2q33.1 | frizzled class receptor 7 |
201792_at | AEBP1 | 165 | 7p13 | AE binding protein 1 |
221766_s_at | FAM46A | 55603 | 6q14.1 | family with sequence similarity 46 member A |
202718_at | IGFBP2 | 3485 | 2q35 | insulin like growth factor binding protein 2 |
203729_at | EMP3 | 2014 | 19q13.33 | epithelial membrane protein 3 |
212063_at | CD44 | 960 | 11p13 | CD44 molecule (Indian blood group) |
203504_s_at | ABCA1 | 19 | 9q31.1 | ATP binding cassette subfamily A member 1 |
201761_at | MTHFD2 | 10797 | 2p13.1 | methylenetetrahydrofolate dehydrogenase (NADP+ dependent) 2, methenyltetrahydrofolate cyclohydrolase |
§ Affymetrix Human Genome Array probe set identifier.
* NCBI Gene Database https://www.ncbi.nlm.nih.gov/gene.
DISCUSSION
For the prognostic biomarker discovery in glioblastomas, there remains discordance among scientific publications given various sample requirements, data complexity, evolving technologies and lack of golden standard practice guidelines [7, 20]. However, due to high complexity of genomics in glioma patients, considerable challenges still present in implementing such strategies. To address these challenges, we applied statistical modelling approaches and leveraged large scale of gene expression profiles from multi-institutional cohorts. The robustness of the GLO10 score was demonstrated by successfully implemented in more than three independent cohorts.
Functional analyses revealed that genes consisting the GLO10 score were involved in the genesis and progression of glioblastomas (Figures 2 and 3). Among them, the prognostic associations with glioblastoma outcomes have already been reported. For example, Insulin-like growth factor-binding protein 2 (IGFBP2) has been considered as a glioma oncogene [21]. Increasing expression of IGFBP2 could associated to poor glioma prognosis as it may play major role in glioma tumour progression [22–27]. It is also reported that in IDH-mutant glioma, IGFBP2 was inhibited so that patients’ survival could be improved [28]. The expressions of FZD7, along with other two genes, SFRP1 and SFRP4, were identified to be associated with poor prognosis in glioma patients [29]. The up-regulated expression of FZD7 could promote glioma cell proliferation [30]. However, to our understanding, GLO10 the first gene signature that utilize the expression data and sum into a unified score that ease its applications in clinical practices.
For the future work, we will collect more datasets and re-train the model to improve the robustness and predictive power. We would also assess its prognostic value by comparing other biomarkers such as IDH1 and 1p19q status. In summary, our analysis of data from different independent cohorts demonstrates that the utility of GLO10 score as a tool for glioblastoma prognosis prediction. The incorporation of the GLO10 into the prognosis prediction for newly diagnosed glioma patients will facilitate the development of biomarker and drug target discovery.
MATERIALS AND METHODS
Datasets collection and preprocess
In this study, we collected datasets from various cohorts including The Cancer Genome Atlas (TCGA)[31], Repository for Molecular Brain Neoplasia Data (Rembrandt)[13] and other published studies. Gene expression profiling datasets and clinical information of five independent cohorts were retrieved NCBI GEO Database with accession numbers: GSE68848 (Rembrandt cohort), GSE83130 (TCGA cohort), GSE13041 [11], GSE16011 [32] and GSE83294 [8]. Given the complexity of data sets, we restricted the samples to grade three and four glioblastomas in this study.
For the Rembrandt cohort, a total to 28 normal and 228 glioblastoma samples were considered. For the glioblastoma prognostic gene signature training, the numbers of samples considered for this study were 470, 160, 155 and 59 for the TCGA cohort (discovery set), GSE13041 (validation set 1), GSE16011 (validation set 2) and GSE83294 (validation set 3). While searching for the data sets, Affymetrix HG-U133A and HG-U133 Plus 2.0 microarrays were considered. Since they are sharing 22,277 probe sets, the performances of the derived gene signatures could be easily evaluated in different studies without gene ID mapping.
The raw fluorescence intensity profiles (*.CEL) were preprocessed, background corrected and normalized with RMA algorithm [33] using Bioconductor package, affy [34], in the R environment.
Differential gene expression analysis
We applied the differentially expressed gene analysis by using the limma (linear models for microarray data, [35]) algorithm. All the gene expression data was converted into the log base-2 scale before comparison using the Welch’s t-test with Benjamini & Hochberg correction [36]. The functional analysis was based on the gene set enrichment analysis in Gene Ontology and Reactome Pathway database [37].
Statistical modelling
For the development of prognostic gene signature, we used the datasets from the discovery set (TCGA cohort, n=470). A linear regression modelling based on the LASSO algorithm was implemented by the glmnet [38] R package. Briefly, the glmnet package fits a generalized linear model via penalized maximum likelihood and extract the prognostic gene signature by10-fold cross validation approaches. A subset of 10 genes was selected as their weighted combined gene expression data was significantly correlated to the overall survival outcomes of glioblastoma patients in the discovery set.
CONFLICTS OF INTEREST
The authors do not have any conflicts of interest.
REFERENCES
1. Furnari FB, Fenton T, Bachoo RM, Mukasa A, Stommel JM, Stegh A, Hahn WC, Ligon KL, Louis DN, Brennan C, Chin L, DePinho RA, Cavenee WK. Malignant astrocytic glioma: genetics, biology, and paths to treatment. Genes Dev. 2007; 21:2683-2710.
2. Yan W, Zhang W, Jiang T. Oncogene addiction in gliomas: implications for molecular targeted therapy. J Exp Clin Cancer Res. 2011; 30:58.
3. Louis DN, Ohgaki H, Wiestler OD, Cavenee WK, Burger PC, Jouvet A, Scheithauer BW, Kleihues P. The 2007 WHO classification of tumours of the central nervous system. Acta Neuropathol. 2007; 114:97-109.
4. Wang Y, Jiang T. Understanding high grade glioma: molecular mechanism, therapy and comprehensive management. Cancer Lett. 2013; 331:139-146.
5. van den Bent MJ. Interobserver variation of the histopathological diagnosis in clinical trials on glioma: a clinician’s perspective. Acta Neuropathol. 2010; 120:297-304.
6. Perry A, Wesseling P. Histologic classification of gliomas. Handb Clin Neurol. 2016; 134:71-95.
7. Weigelt B, Baehner FL, Reis-Filho JS. The contribution of gene expression profiling to breast cancer classification, prognostication and prediction: a retrospective of the last decade. J Pathol. 2010; 220:263-280.
8. Freije WA, Castro-Vargas FE, Fang Z, Horvath S, Cloughesy T, Liau LM, Mischel PS, Nelson SF. Gene expression profiling of gliomas strongly predicts survival. Cancer Res. 2004; 64:6503-6510.
9. Nutt CL, Mani DR, Betensky RA, Tamayo P, Cairncross JG, Ladd C, Pohl U, Hartmann C, McLaughlin ME, Batchelor TT, Black PM, von Deimling A, Pomeroy SL, et al. Gene expression-based classification of malignant gliomas correlates better with survival than histological classification. Cancer Res. 2003; 63:1602-1607.
10. Phillips HS, Kharbanda S, Chen R, Forrest WF, Soriano RH, Wu TD, Misra A, Nigro JM, Colman H, Soroceanu L, Williams PM, Modrusan Z, Feuerstein BG, Aldape K. Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis. Cancer Cell. 2006; 9:157-173.
11. Lee Y, Scheck AC, Cloughesy TF, Lai A, Dong J, Farooqi HK, Liau LM, Horvath S, Mischel PS, Nelson SF. Gene expression analysis of glioblastomas identifies the major molecular basis for the prognostic benefit of younger age. BMC Med Genomics. 2008; 1:52.
12. Liang Y, Diehn M, Watson N, Bollen AW, Aldape KD, Nicholas MK, Lamborn KR, Berger MS, Botstein D, Brown PO, Israel MA. Gene expression profiling reveals molecularly and clinically distinct subtypes of glioblastomamultiforme. Proc Natl Acad Sci U S A. 2005; 102:5814-5819.
13. Madhavan S, Zenklusen JC, Kotliarov Y, Sahni H, Fine HA, Buetow K. Rembrandt: helping personalized medicine become a reality through integrative translational research. Mol Cancer Res. 2009; 7:157-167.
14. Blanchart A, Fernando R, Haring M, Assaife-Lopes N, Romanov RA, Andang M, Harkany T, Ernfors P. Endogenous GABAA receptor activity suppresses glioma growth. Oncogene. 2017; 36:777-786.
15. McFerrin MB, Sontheimer H. A role for ion channels in glioma cell invasion. Neuron Glia Biol. 2006; 2:39-49.
16. Joshi AD, Parsons DW, Velculescu VE, Riggins GJ. Sodium ion channel mutations in glioblastoma patients correlate with shorter survival. Mol Cancer. 2011; 10:17.
17. Netti PA, Berk DA, Swartz MA, Grodzinsky AJ, Jain RK. Role of extracellular matrix assembly in interstitial transport in solid tumors. Cancer Res. 2000; 60:2497-2503.
18. Hynes RO. The extracellular matrix: not just pretty fibrils. Science. 2009; 326:1216-1219.
19. Ulrich TA, de Juan Pardo EM, Kumar S. The mechanical rigidity of the extracellular matrix regulates the structure, motility, and proliferation of glioma cells. Cancer Res. 2009; 69:4167-4174.
20. Subramanian J, Simon R. What should physicians look for in evaluating prognostic gene-expression signatures? Nat Rev Clin Oncol. 2010; 7:327-334.
21. Holmes KM, Annala M, Chua CY, Dunlap SM, Liu Y, Hugen N, Moore LM, Cogdell D, Hu L, Nykter M, Hess K, Fuller GN, Zhang W. Insulin-like growth factor-binding protein 2-driven glioma progression is prevented by blocking a clinically significant integrin, integrin-linked kinase, and NF-kappaB network. Proc Natl Acad Sci U S A. 2012; 109:3475-3480.
22. Wang H, Wang H, Zhang W, Fuller GN. Overexpression of IGFBP5, but not IGFBP3, correlates with the histologic grade of human diffuse glioma: a tissue microarray and immunohistochemical study. Technol Cancer Res Treat. 2006; 5:195-199.
23. McDonald KL, O’Sullivan MG, Parkinson JF, Shaw JM, Payne CA, Brewer JM, Young L, Reader DJ, Wheeler HT, Cook RJ, Biggs MT, Little NS, Teo C, et al. IQGAP1 and IGFBP2: valuable biomarkers for determining prognosis in glioma patients. J Neuropathol Exp Neurol. 2007; 66:405-417.
24. Zhou YH, Hess KR, Raj VR, Yu L, Liu L, Yung AW, Linskey ME. Establishment of prognostic models for astrocytic and oligodendroglial brain tumors with standardized quantification of marker gene expression and clinical variables. Biomark Insights. 2010; 5:153-168.
25. Zheng S, Houseman EA, Morrison Z, Wrensch MR, Patoka JS, Ramos C, Haas-Kogan DA, McBride S, Marsit CJ, Christensen BC, Nelson HH, Stokoe D, Wiemels JL, et al. DNA hypermethylation profiles associated with glioma subtypes and EZH2 and IGFBP2 mRNA expression. Neuro Oncol. 2011; 13:280-289.
26. Ahani N, Karimi Arzenani M, Shirkoohi R, Rokouei M, Alipour Eskandani M, Nikravesh A. Expression of insulin-like growth factor binding protein-2 (IGFBP-2) gene in negative and positive human cytomegalovirus glioblastomamultiforme tissues. Med Oncol. 2014; 31:812.
27. Martino-Echarri E, Fernandez-Rodriguez R, Bech-Serra JJ, Plaza-Calonge Mdel C, Vidal N, Casal C, Colome N, Seoane J, Canals F, Rodriguez-Manzaneque JC. Relevance of IGFBP2 proteolysis in glioma and contribution of the extracellular protease ADAMTS1. Oncotarget. 2014; 5:4295-4304. https://doi.org/10.18632/oncotarget.2009.
28. Huang LE, Cohen AL, Colman H, Jensen RL, Fults DW, Couldwell WT. IGFBP2 expression predicts IDH-mutant glioma patient survival. Oncotarget. 2017; 8:191-202. https://doi.org/10.18632/oncotarget.13329.
29. Kierulf-Vieira KS, Sandberg CJ, Grieg Z, Gunther CC, Langmoen IA, Vik-Mo EO. Wnt inhibition is dysregulated in gliomas and its re-establishment inhibits proliferation and tumor sphere formation. Exp Cell Res. 2016; 340:53-61.
30. Qiu X, Jiao J, Li Y, Tian T. Overexpression of FZD7 promotes glioma cell proliferation by upregulating TAZ. Oncotarget. 2016; 7:85987-85999. https://doi.org/10.18632/oncotarget.13292.
31. Cancer Genome Atlas Research Network. Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature. 2008; 455:1061-1068.
32. Gravendeel LA, Kouwenhoven MC, Gevaert O, de Rooi JJ, Stubbs AP, Duijm JE, Daemen A, Bleeker FE, Bralten LB, Kloosterhof NK, De Moor B, Eilers PH, van der Spek PJ, et al. Intrinsic gene expression profiles of gliomas are a better predictor of survival than histology. Cancer Res. 2009; 69:9065-9072.
33. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003; 4:249-264.
34. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004; 20:307-315.
35. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47.
36. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995:289-300.
37. Fabregat A, Sidiropoulos K, Viteri G, Forner O, Marin-Garcia P, Arnau V, D’Eustachio P, Stein L, Hermjakob H. Reactome pathway analysis: a high-performance in-memory approach. BMC Bioinformatics. 2017; 18:142.
38. Friedman JH, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. 2010; 33:22.
