A multiple genomic data fused SF2 prediction model, signature identification, and gene regulatory network inference for personalized radiotherapy.

Radiotherapy is one of the most important cancer treatments, but its response varies greatly among individual patients. Therefore, the prediction of radiosensitivity, identification of potential signature genes, and inference of their regulatory networks are important for clinical and oncological reasons. Here, we proposed a novel multiple genomic fused partial least squares deep regression method to simultaneously analyze multi-genomic data. Using 60 National Cancer Institute cell lines as examples, we aimed to identify signature genes by optimizing the radiosensitivity prediction model and uncovering regulatory relationships. A total of 113 signature genes were selected from more than 20,000 genes. The root mean square error of the model was only 0.0025, which was much lower than previously published results, suggesting that our method can predict radiosensitivity with the highest accuracy. Additionally, our regulatory network analysis identified 24 highly important 'hub' genes. The data analysis workflow we propose provides a unified and computational framework to harness the full potential of large-scale integrated cancer genomic data for integrative signature discovery. Furthermore, the regression model, signature genes, and their regulatory network should provide a reliable quantitative reference for optimizing personalized treatment options, and may aid our understanding of cancer progress mechanisms.

Introduction

Radiotherapy is a major cancer treatment, but the radiosensitivity of different tumors or even the same type of tumor in different patients varies widely.[1] Therefore, predicting the radiosensitivity of patients before radiation therapy, identifying underlying molecular signatures, and constructing their regulatory network have high clinical and oncological importance.

Radiosensitivity can be measured as the fraction of cells surviving a single 2 Gy dose of radiation (SF2), with high values indicating radiation resistance. While other methods are available to measure cellular radiation sensitivity in cell lines, SF2 is considered to be the gold standard and is supported by strong clinical evidence.[2]

Variations in mRNA expression values (GE), copy numbers, and other genomic patterns are thought to be the main underlying factors for different radiation responses. Accumulating large amounts of these data provides an effective but challenging way to predict the radiosensitivity of tumor cells.

Torres-Roca et al predicted the radiosensitivity of 35 human cell lines in a NCI-60 panel using a linear classifier of expression values of tens of genes selected by the significance analysis of microarrays (SAM) method in 2005.[3] They developed a radiosensitivity index (RSI) as a biomarker of cellular radiosensitivity in 48 NCI-60 cancer cell lines in 2009. SF2 was the central criterion for both feature selection and final model training for RSI development. Ten of the selected ‘hub’ genes were then used to construct a linear prediction model of SF2.[4] Additionally, Tewari et al investigated the feasibility of integrating an in vitro chemo-radiosensitivity assay with a gene microarray system,[5] identifying 54 genes correlated with radiosensitivity using an integrated nearest neighbor model with Pearson correlation coefficient. Moreover, Amundson et al performed large-scale comparisons of gene expression variations in response to different doses (2, 5, and 8 Gy) of γ-ray radiation,[6] and identified 22 genes that could discriminate the SF2 values of 63 cell lines (including NCI-60 and three other cell lines) between low and high groups.

Besides GE data, copy number variation (CNV) and methylation (ME) data are also related to radiosensitivity. Work by Moelans et al indicated that allelic loss and amplification at the 8p11-12 breakpoint region are associated with poor radiotherapy responses,[7] while Zhu et al reported a pivotal role for DNA ME in tumor radiosensitivity.[8]

Unfortunately, none of the individual types of genomic data thoroughly capture the complexity of the cancer genome or precisely pinpoint the cancer-driving mechanism.[9] Additionally, it has become increasingly clear that the integrative analysis of multi-omic data types can aid the search for potential “drivers” by uncovering genomic features dysregulated by multiple mechanisms.[10] More importantly, true oncogenic mechanisms are more visible when combining evidence across patterns of alterations in DNA CNV, ME, GE and mutational profiles.[11,12] A well-known example is the HER2 oncogene that can be activated through DNA amplification and mRNA over-expression.[9] Therefore, the development of tumor molecular analysis using multiple genomic data may lead to a more comprehensive prediction of molecular signatures.

No widely accepted threshold exists between radiation-sensitive or -resistant phenotypes and SF2 values from 0 to 1. Therefore, instead of roughly dividing cell lines into different groups by subjective cutoffs of SF2 values, it may be more useful to consider a regression issue for continuous variables of SF2. To this end, we propose a novel integrated multiple genomic data regression method for SF2 prediction, focusing on identifying signature genes for functional and genetic network analysis, rather than “hotspot” or “hot-loci” from GE, CNV, and ME data.[13,14]

Studying the gene regulatory network (GRN) structure provides important insights into the mechanisms of complex diseases.[15,16] Several studies have shown that gene expression is influenced not only by the expression of other genes but also by CNV or other biological variations.[17] Therefore, it is also necessary to infer GRN using multi-genomic data. Correspondingly, the aims of this study were two-fold: 1) to identify signature genes strongly associated with radiosensitivity from fused multiple genomic data to further corroborate and expand the evidence of radiosensitivity-associated signature genes in the prediction of radioresponses; and 2) to uncover regulatory relationships among identified signatures using fused multiple genomic data, employing least absolute shrinkage and selection operator (LASSO) regression based on coordinate descent algorithms to construct GRN. Figure 1 shows the study outline.

Figure 1.

The overall structure of this paper. Based on the fact that the current standard approaches rely on separate mono-genomics data analyses followed by manual integration, multiple genomic data fused regression approach (MGPLS) is proposed to identifying signature genes. MGPLS method can analysis all data types simultaneously using a single integrated regression model as well as eliminating noise effects. VIP: variable importance on projection; CV: cross-validation; UVE: uninformative variable elimination.

Materials and Methods

Datasets

GE, CNV, and ME data of NCI-60 cell lines were downloaded from (https://discover.nci.nih.gov/cellminer/loadDownload.do). GE data were collected via five platforms (Affymetrix: HG-U95, HG-U133, HG-U133 Plus 2.0, HG Exon 1.0 ST; and Agilent: Whole Human Genome Oligo array).[18] CNV data were collected via four different platforms (Agilent Human Genome CGH Microarray 44A, Nimblegen HG19 CGH 385K WG Tiling v2.0, Affymetrix GeneChip Human Mapping 500k Array Set, and Illumina Human1Mv1_C Beadchip).[19] ME data were collected using the Infinium HumanMethylation450 BeadChip Kit platform. The measured SF2 values of corresponding cell lines were collected from the study of Eschrich et al[4] and are listed in Table 1.

Table 1.

Measured and predicted SF2 values for NCI-60 cell lines

Cell line	Measured SF2	Predicted SF2(MGPLS)	Predicted SF2(SVM)	Error(MGPLS)	Error(SVM)
CNS:U251	0.57	0.568	0.57	0.002	0
OV:OVCAR-4	0.29	0.293	0.291	-0.003	-0.001
LE:CCRF-CEM	0.185	0.180	0.186	0.005	-0.001
CNS:SNB-19	0.43	0.425	0.43	0.005	0
RE:SN12C	0.62	0.611	0.618	0.009	0.002
BR:T-47D	0.52	0.531	0.52	-0.011	0
RE:ACHN	0.72	0.707	0.722	0.013	-0.002
LE:HL-60(TB)	0.315	0.335	0.319	-0.020	-0.004
ME:MALME-3M	0.8	0.779	0.797	0.021	0.003
ME:SK-MEL-5	0.72	0.697	0.72	0.023	0
OV:OVCAR-8	0.6	0.572	0.597	0.028	0.003
ME:UACC-257	0.48	0.510	0.48	-0.030	0
ME:SK-MEL-28	0.74	0.709	0.737	0.031	0.003
LE:RPMI-8226	0.1	0.069	0.1	0.031	0
PR:DU-145	0.52	0.488	0.52	0.032	0
BR:HS 578T	0.79	0.757	0.79	0.033	0
CO:HCT-15	0.4	0.435	0.4	-0.035	0
CO:HCT-116	0.38	0.418	0.38	-0.038	0
PR:PC-3	0.484	0.445	0.486	0.039	-0.002
LC:EKVX	0.7	0.660	0.694	0.040	0.006
OV:OVCAR-5	0.408	0.452	0.409	-0.044	-0.001
RE:TK-10	0.52	0.475	0.522	0.045	-0.002
LE:K-562	0.05	0.100	0.054	-0.050	-0.004
OV:NCI/ADR-RES	0.57	0.520	0.572	0.050	-0.002
CNS:SNB-75	0.55	0.602	0.55	-0.052	0
ME:M14	0.42	0.477	0.42	-0.057	0
ME:UACC-62	0.52	0.461	0.519	0.059	0.001
OV:OVCAR-3	0.55	0.491	0.548	0.059	0.002
LC:NCI-H322M	0.65	0.587	0.65	0.063	0
RE:UO-31	0.62	0.686	0.619	-0.066	0.001
CO:COLO 205	0.69	0.762	0.687	-0.072	0.003
OV:IGROV1	0.39	0.463	0.39	-0.073	0
LE:SR	0.07	0.143	0.072	-0.073	-0.002
CNS:SF-539	0.82	0.746	0.817	0.074	0.003
BR:MCF7	0.576	0.500	0.574	0.076	0.002
ME:SK-MEL-2	0.66	0.737	0.66	-0.077	0
RE:RXF 393	0.67	0.754	0.669	-0.084	0.001
LC:NCI-H522	0.43	0.344	0.431	0.086	-0.001
ME:LOX IMVI	0.68	0.588	0.68	0.092	0
LC:HOP-92	0.43	0.522	0.43	-0.092	0
CNS:SF-268	0.45	0.543	0.45	-0.093	0
ME:MDA-MB-435	0.1795	0.273	0.183	-0.094	-0.0035
BR:BT-549	0.632	0.537	0.635	0.095	-0.003
ME:MDA-N	0.45	0.352	0.449	0.098	0.001
CNS:SF-295	0.73	0.631	0.73	0.099	0
LE:MOLT-4	0.05	0.149	0.052	-0.099	-0.002
RE:786-0	0.66	0.551	0.659	0.109	0.001
LC:HOP-62	0.164	0.277	0.166	-0.113	-0.002
CO:KM12	0.42	0.535	0.418	-0.115	0.002
LC:A549/ATCC	0.61	0.730	0.61	-0.120	0
RE:A498	0.61	0.734	0.62	-0.124	-0.001
CO:HCC-2998	0.44	0.572	0.439	-0.132	0.001
OV:SK-OV-3	0.9	0.767	0.894	0.133	0.006
RE:CAKI-1	0.37	0.517	0.37	-0.147	0
CO:SW-620	0.62	0.473	0.622	0.147	-0.002
LC:NCI-H226	0.63	0.786	0.626	-0.156	0.004
LC:NCI-H460	0.84	0.671	0.835	0.169	0.005
BR:MDA-MB-231	0.82	0.613	0.82	0.207	0
CO:HT29	0.79	0.567	0.785	0.223	0.005
LC:NCI-H23	0.086	0.315	0.0925	-0.229	-0.0065

* Cell lines sorted by multiple genomic data fused partial least square deep regression (MGPLS).

Data preprocessing

Before training the dataset, we performed the following preprocessing on the downloaded data:

Considering that there are only 3–4 samples per cancer type in the NCI-60 cell line panel, if the value of a gene is missing in one sample then it is absent from >25% of the samples of that cancer type. Therefore, these genes were removed from the analysis.

Common genes in all GE, CNV, and ME datasets were identified.

Because GEs and CNVs of genes that are strongly related to each other are more likely to contain less noise,[19] they were selected using Pearson correlation coefficients (cutoff, 0.5).

GE, CNV, and ME data were respectively standardized with the Z-Score method for further regression analysis.

A multiple genomic data fused partial least square deep regression (MGPLS) method

Partial least squares (PLS) is a widely used algorithm for modeling relationships between sets of observed variables using latent variables. It comprises regression and classification tasks as well as dimension reducing and modeling.[20] Instead of identifying hyperplanes of minimum variance between the response and independent variables, it finds a linear regression model by projecting the predicted variables (regression or classification labels) and the observed variables (fused GE, CNV, or ME values of genes in our case) to a new lower space.[21] This is highly suited to the analysis of high-dimension, low-sample-size data in bioinformatics.

To integrate multiple genome data for improved regression performance, we proposed an MGPLS method:

Given the predictor matrix and the response matrix ,

where p = g + c + m, g is the number of variables in GE data, c is the number of variables in CNV data, and m is the number of variables in ME data.

MGPLS may be applied to cases where the aim is to describe the linear relationship between and ,

based on the basic latent variable decomposition:

where is the regression coefficient matrix, is the residual matrix, is the orthogonal latent variable (LV) matrix, and are loading matrices, and are residual matrices, and k is the number of LVs. According to the regular PLS, is a linear transformation of ,

where is a matrix of weights.

From Eq.(1) and Eq.(5), it is obvious that

where is the ith LV (i.e., ith column of ), and w is ith column of . Hence, by the projection of the MGPLS algorithm, the p-dimensional -space, consisting of GE, CNV, and ME, is integrated and compressed into the k-dimensional LV-space (k<