The present invention relates to a method for identifying radiation induced genes and long non-coding RNAs (lncRNAs). In particular, some of the radiation induced genes and long non-coding RNAs (lncRNAs) are applied as markers for predicting radiotherapy response of glioblastoma in patients with glioblastoma.
Radiotherapy has become a popular and standard approach for treating cancer patients because it greatly improves patient survival. However, some of the patients receiving radiotherapy suffer from adverse effects and do not obtain survival benefits. This may be attributed to the fact that most radiation treatment plans are designed based on cancer type, without consideration of each individual's radiosensitivity.
Conventional radiotherapy is now combined with images processed by computed tomography to provide a non-invasive and highly tumor-specific treatment plan. Currently, most radiation treatment plans are designed based purely on the cancer type. Challenges arise, however, when individual genetic differences in radiosensitivity lead to different treatment responses among patients receiving the same radiation dose. For example, International Journal of Radiation Oncology* Biology* Physics. 2005; 61(1):196-202 disclosed genetic variants in ATM are associated with differential hypersensitivity to radiation exposure. These results suggest that radiation sensitivity should be taken into consideration when designing the treatment plan.
Based on the aforementioned description, a method for identifying markers for identifying radiation induced markers and then predicting the radiosensitivity of tumor in patients with tumor by the radiation induced markers is required in the near future.
In one aspect, the present invention provides a method for identifying radiation induced genes and long non-coding RNAs (lncRNAs). The invention method comprises the steps of (1). Provide expression values of genes and long non-coding RNAs; (2). Execute weighted gene correlation network analysis (WGCNA) by a computer system to calculate Pearson correlation coefficients of pairs of the genes and long non-coding RNAs based on the expression values of the genes and long non-coding RNAs; and (3). Perform a screening step by the computer system to identify radiation induced genes and long non-coding RNAs based on the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs with a value more than 0.75.
The aforementioned radiation induced genes and long non-coding RNAs are identified from the pairs of the genes and long non-coding RNAs based on the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs with a value more than 0.75. That is the Pearson correlation coefficients of the pairs of the radiation induced genes and long non-coding RNAs are more than 0.75.
In general, the expression values of genes and long non-coding RNAs are measured in vitro from immortalized B cells (GSE26835) after irradiation by a method comprises microarray.
For developing the method for identifying radiation induced genes and long non-coding RNAs, some tests and algorithm, such as the analysis of variance (ANOVA) test, Tukey's honest significant difference (HSD) test and weighted gene correlation network analysis are performed by a computer system for analyzing the expression values of genes and long non-coding RNAs, and further identify the radiation induced genes and long non-coding RNAs based on the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs.
Firstly, the analysis of variance (ANOVA) test and Tukey's honest significant difference (HSD) test were performed in microarray dataset of immortalized B cells (GSE26835) to identify differentially expressed genes and lncRNAs (P<0.0001). In total, 1,086 samples were classified into 3 groups, which were 0, 2, and 6 h after irradiation (As shown in
Secondly, a network-based co-expression algorithm, weighted gene correlation network analysis (WGCNA), was performed to identify modules composed of correlated genes and lncRNAs. Comparison of the samples without irradiation (0 h) to the identified set of differentially expressed probes (N=1,244) via WGCNA revealed 2 and 10 modules at 2 h and 6 h post-irradiation, respectively. Among the 12 identified modules, highly correlated interactions between genes and lncRNAs were selected based on their Pearson correlation coefficients (r>0.75). Consequently, a total of 43 interaction pairs triggered by irradiation were identified (As listed in Table 1), which contains 34 genes and 7 lncRNAs. The 34 genes and 7 lncRNAs are determined to the claimed radiation induced genes and lncRNAs.
In another aspect, the present invention provides a method for predicting radiotherapy response of glioblastoma in patients with glioblastoma. The method comprises the steps of: (1). determine expression values of markers in a test sample from the patients with glioblastoma to obtain a test dataset of expression values of the markers, wherein the markers are selected from the group consisting of DDX3Y, EIF1AY, RPS4Y1, USP9Y, RHOC, EPS8L2, ACTA2, MR1, TRAPPC6A, ANXA4, ETHE1, PYCARD, JUP///KRT17, ETFDH, ALDH6A1, RTN1, TP53TG1, LOC100653017 and LOC100506948; (2). compare the test dataset of expression values of the markers to a radiosensitive sample and to a radioresistant sample; (3). classify whether a test dataset of expression values of the markers is significantly within the radiosensitive sample or within the radioresistant sample to assess whether the test sample is radiosensitive (RS) or radioresistant (RR); (4). determine the patients with glioblastoma as in a radiosensitive (RS) group if the test sample is radiosensitive, or determine the patients with glioblastoma as in a radioresistant (RR) group if the test sample is radioresistant; and (5). predict radiotherapy response of glioblastoma in the patients with glioblastoma based on whether the radiosensitive (RS) group or radioresistant (RR) group is treated with radiotherapy.
For developing the aforementioned method, the NCI-60 cell lines, along with the fraction of surviving cells after a radiation dose of 2 Gy (SF2), were utilized as the training set. For the claimed radiation induced genes and long non-coding RNAs, a genetic algorithm (GA) was used to select the best combination of genes and lncRNAs as the markers for developing a model by the support vector machine (SVM) algorithm (As shown in
The SVM algorithm is applied to calculate a decision value for determining whether a sample generated from the expression of the 20 markers as listed in Table 2 in NCI-60 cell lines is radiosensitive (RS) or radioresistant (RR). When the decision value is negative, a sample is defined to a radiosensitive sample. In contrast, when the decision value is positive, a sample is defined to a radioresistant sample.
Two microarray datasets from patients with glioblastoma multiforme (GBM) were retrieved, one from The Cancer Genome Altas (TCGA) and another from the Gene Expression Omnibus (GEO) (accession number GSE16011). The datasets were analyzed to evaluate the prediction performance of the developed SVM model. The model was applied to classify the patients into radiosensitive (RS) and radioresistant (RR) groups.
To sum up, the present invention provides the method for identifying radiation induced genes and lncRNAs. In particular, some of the radiation induced genes and lncRNAs are applied as markers for predicting the radiosensitivity of glioblastoma in the patients with glioblastoma. Also the method for predicting radiotherapy response of glioblastoma in patients with glioblastoma is disclosed in the present invention.
In a first embodiment, the present invention discloses a method for identifying radiation induced genes and long non-coding RNAs (lncRNAs). The invention method comprises the steps of (1). Provide expression values of genes and long non-coding RNAs; (2). Execute weighted gene correlation network analysis (WGCNA) by a computer system to calculate Pearson correlation coefficients of pairs of the genes and long non-coding RNAs based on the expression values of the genes and long non-coding RNAs; and (3). Perform a screening step by the computer system to identify the radiation induced genes and long non-coding RNAs based on the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs with a value more than 0.75.
Accordingly, the invention method for identifying the radiation induced genes and long non-coding RNAs is mainly based on the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs more than 0.75. That is the radiation induced genes and long non-coding RNAs have the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs more than 0.75.
In one example of the first embodiment, the expression values of genes and long non-coding RNAs are measured in vitro from immortalized B cells (GSE26835) after irradiation by a method comprises microarray.
In one example of the first embodiment, the radiation induced genes are selected from the group consisting of DDX3Y, EIF1AY, RPS4Y1, USP9Y, KDM5D, TP53I3, RHOC, ASTN2, GAMT, EPS8L2, ACTA2, E2F5, MR1, UTY, TRAPPC6A, FHL2, ANXA4, GLS2, CEACAM1, ETHE1, TSPAN31, PYCARD, CDK2, JUP///KRT17, ATP6V1D, PROCR, ETFDH, ALDH6A1 and RTN1.
In one example of the first embodiment, the radiation induced long non-coding RNAs are selected from the group consisting of TTTY15, TP53TG1, LOC100653079, LOC100653017 and LOC100506948.
For identifying the claimed radiation induced genes and long non-coding RNAs, some tests and algorithm, such as the analysis of variance (ANOVA) test, Tukey's honest significant difference (HSD) test and weighted gene correlation network analysis (WGCNA) are executed by a computer system for analyzing the expression values of genes and long non-coding RNAs, and further identify the radiation induced genes and long non-coding RNAs based on the Pearson correlation coefficients of the pairs of the genes and long non-coding RNAs.
Firstly, the analysis of variance (ANOVA) test and Tukey's honest significant difference (HSD) test were performed in microarray dataset of immortalized B cells (GSE26835) to identify differentially expressed genes and lncRNAs (P<0.0001). In total, 1,086 samples were classified into 3 groups, which were 0, 2, and 6 h after irradiation (As shown in
Secondly, a network-based co-expression algorithm, weighted gene correlation network analysis (WGCNA), was performed to identify modules composed of correlated genes and lncRNAs. Comparison of the samples without irradiation (0 h) to the identified set of differentially expressed probes (N=1,244) via WGCNA revealed 2 and 10 modules at 2 h and 6 h post-irradiation, respectively. Among the 12 identified modules, highly correlated interactions between genes and lncRNAs were selected based on their Pearson correlation coefficients (r>0.75). Consequently, a total of 43 interaction pairs triggered by irradiation were identified as listed in Table 1, which contains 34 genes and 7 lncRNAs. The 34 genes and 7 lncRNAs are defined as the claimed radiation induced genes and lncRNAs.
In one preferred example of the first embodiment, the radiation induced genes selected from the group consisting of DDX3Y, EIF1AY, RPS4Y1, USP9Y, RHOC, EPS8L2, ACTA2, MR1, TRAPPC6A, ANXA4, ETHE1, PYCARD, JUP///KRT17, ETFDH, ALDH6A1 and RTN1 are applied as markers for patients with glioblastoma.
In one preferred example of the first embodiment, the radiation induced long non-coding RNAs selected from the group consisting of TP53TG1, LOC100653017 and LOC100506948 are applied as markers for patients with glioblastoma.
In a second embodiment of the present invention, a method for predicting radiotherapy response of glioblastoma in patients with glioblastoma is provided. The method comprises the steps of: (1). determine expression values of markers in a test sample from the patients with glioblastoma to obtain a test dataset of expression values of the markers, wherein the markers are selected from the group consisting of DDX3Y, EIF1AY, RPS4Y1, USP9Y, RHOC, EPS8L2, ACTA2, MR1, TRAPPC6A, ANXA4, ETHE1, PYCARD, JUP///KRT17, ETFDH, ALDH6A1, RTN1, TP53TG1, LOC100653017 and LOC100506948; (2). compare the test dataset of expression values of the markers to a radiosensitive sample and to a radioresistant sample; (3). classify whether a test dataset of expression values of the markers is significantly within the radiosensitive sample or within the radioresistant sample to assess whether the test sample is radiosensitive or radioresistant; (4). determine the patients with glioblastoma as in a radiosensitive (RS) group if the test sample is radiosensitive, or determine the patients with glioblastoma as in a radioresistant (RR) group if the test sample is radioresistant; and (5). predict radiotherapy response of glioblastoma in the patients with glioblastoma based on whether the radiosensitive (RS) group or radioresistant (RR) group is treated with radiotherapy.
For developing the aforementioned method, the NCI-60 cell lines, along with the fraction of surviving cells after a radiation dose of 2 Gy (SF2), were utilized as the training set. For the claimed radiation induced genes and long non-coding RNAs, a genetic algorithm (GA) was used to select the best combination of genes and lncRNAs as the markers for developing a model by the support vector machine (SVM) algorithm (As shown in
The SVM algorithm is applied to calculate a decision value for determining whether a sample generated from the expression of the 20 markers as listed in Table 2 in NCI-60 cell line is radiosensitive (RS) or radioresistant (RR).
The SVM feature weight (w) is firstly calculated by the following function using a computer system with SVM algorithm. The coefficients (coefs) of the 20 markers are listed in TABLE 3
SVM feature weight (w)=t(model$coefs) %*% model SV
The decision value is then obtained by the following function and equation using a computer system with SVM algorithm.
Data.scaled=scale(testing_data,model$x.scale[[1]],mode.$x.scale[[2]])
Decision value=t(w %*% t(as.matrix(data.scaled)))−model$rho
When the decision value is negative (decision value <0), the sample is defined to a radiosensitive sample. In contrast, when the decision value is positive (decision value >0), the sample is defined to a radioresistant sample.
Two microarray datasets from patients with glioblastoma multiforme (GBM) were retrieved, one from The Cancer Genome Altas (TCGA) and another from the Gene Expression Omnibus (GEO) (accession number GSE16011). The datasets were analyzed to evaluate the prediction performance of the developed SVM model. The model was applied to classify the patients into radiosensitive (RS) and radioresistant (RR) groups.
In one example of the second embodiment, the method for predicting radiotherapy response of glioblastoma in patients with glioblastoma is applied to a prognosis analysis of the patients with glioblastoma treated with radiotherapy.
In one example of the second embodiment, the expression values of the markers are measured in vitro by a method comprises microarray and real-time polymerase chain reaction (RT-PCR).
In one example of the second embodiment, the test sample comprises lymphocytes sample and peripheral blood sample.
In one example of the second embodiment, the radiosensitive sample and radioresistant sample are determined by a decision value calculated from the expression values of the markers in NCI-60 cell lines by support vector machine, when the decision value is negative, the radiosensitive sample is defined; and when the decision value is positive, the radioresistant sample is defined.
In one example of the second embodiment, the step of classifying whether the test dataset of expression values of the markers is significantly within the radiosensitive sample or within the radioresistant sample comprises use of a predictive algorithm. Preferably, the predictive algorithm is a support vector machine.
In one example of the second embodiment, the steps of (2), (3), (4) and (5) are performed by a computer system.
The following descriptions are to explain both of the method for identifying the claimed radiation induced genes and long non-coding RNAs (lncRNAs) and the method for predicting radiotherapy response of glioblastoma in patients with glioblastoma in more details.
Microarray Datasets
Four datasets analyzed in this study were retrieved from the Gene Expression Omnibus (GEO), CellMiner, and The Cancer Genome Altas (TCGA). The dataset with the largest sample size, GSE26835, was utilized as the identification set to select differentially expressed genes and lncRNAs triggered by irradiation. Subsequently, a prediction model for radiosensitivity was developed based on the expression profiles of the NCI-60 cell lines. Two GBM datasets were used as the validation sets to evaluate the performance of the developed prediction model. All preprocessing procedures and normalization algorithms, including robust multiarray averaging (RMA) and quantile normalization, were performed in R. Although the microarray platforms of these 4 datasets were originally designed to examine gene expression profiles, the studies have shown that a re-analysis of the probe sequences can select the probes targeting lncRNAs. In this study, all those probe sets targeting lncRNAs can be mapped to the Ensembl database.
Identification of Differentially Expressed Genes and lncRNAs Triggered by Irradiation
A protocol to identify the radiation-induced genes and lncRNAs and develop a prediction model is illustrated in
Identification of Gene-lncRNA Interaction Pairs
To take the biological functions of lncRNAs into consideration, a network-based co-expression algorithm, WGCNA, was used. The main concept of the WGCNA method is to cluster the genes and lncRNAs with similar expression patterns into a single module based on the hypothesis that genes and lncRNAs involved in the same functional pathway will have highly correlated expression values. All parameters were set as their default values, except that “deepSplit” was used to explore more possible regulations between genes and lncRNAs. Lastly, the Pearson correlation coefficients were calculated for the genes and lncRNAs classified into the same module in order to select the highly correlated gene-lncRNA interaction pairs (r>0.75), and then the claimed radiation induced genes and lncRNAs are identified.
Development of a Prediction Model Using a Genetic Algorithm
The gene expression profiles and radiation parameters of the NCI-60 cell lines retrieved from CellMiner were analyzed to develop a prediction model. To mimic the situation of radiotherapy in real patients, the survival fraction under 2Gy (SF2) was utilized, and the NCI-60 cell lines were dichotomized into radiosensitive (RS) and radioresistant (RR) groups based on the threshold of 0.4. A genetic algorithm (GA) was designed to select genes and lncRNAs with the highest prediction accuracy in the NCI-60 cell lines (As shown in
Validation of Prediction Model in GBM Patients
The developed prediction model was applied to 2 microarray datasets of GBM patients from GSE16011 and TCGA. Considering the high mortality rate of GBM, we focused on the overall survival rate within 2 years. Patients' data were excluded from the analyses if no information was provided about their status of radiotherapy and survival. The stratification resulted in 324 RT+ and 57 RT− patients in the GBM dataset from TCGA and 193 RT+ and 70 RT− patients in GSE16011. The developed prediction model was utilized to classify those patients into RS or RR groups, and log-rank tests were used to examine the differences in overall survival between them. To further evaluate whether the prediction model was associated with radiotherapy only, the RT+ and RT− patients were compared accordingly. Kaplan-Meier survival curves (as shown in
A multivariable Cox hazard regression model was performed to compare the prediction performance of the model with other clinical parameters, including age, grade, and Karnofsky Performance Score (KPS). As shown in Table 4, the prediction model remains an independent predictor after being adjusted for other confounding factors (P=0.00698 for TCGA and P=0.0277 for GSE16011). Therefore, these results demonstrate that the prediction model can effectively identify patients who will be responsive to radiotherapy and thus experience significantly better survival outcomes.
As so far, the present invention is the first study to incorporate the expression levels of lncRNAs into a prediction model for radiosensitivity. Furthermore, a permutation test was performed in this study to assess whether it is superior to add lncRNAs into a prediction model for radiosensitivity.
In the development of a prediction model, a critical step is how to select important markers (predictors), especially when handling high-throughput data. A popular approach is to perform a stepwise regression using forward selection and/or backward elimination, and to evaluate the performance by the Akaike information criterion (AIC). However, such methods are time-consuming and inefficient, because the number of possible combinations of predictors is large. We tried such an approach in this study, but the predictors identified from the AIC model changed substantially upon the addition of only 1 or 2 pairs of genes and lncRNAs. The GA was used instead due to its low computational complexity and high efficiency in achieving good accuracy. In this study, the number of possible combinations of the original pool in this study (34 genes and 7 lncRNAs) is C2041=269128937220, which is impossible to test in a reasonable amount of time. In addition, it is well-known that radiosensitivity differs greatly across different tissues. Therefore, we analyzed the NCI-60 cell lines to develop the prediction model because they were from distinct tissue types. One major advantage of using the GA in this situation is that it can identify the best combination of predictors (genes and lncRNAs) through random selection over many generations, taking the heterogeneity of radiosensitivity in different tissue types into consideration. The randomness in the initial step of GA is compensated for by the reproducibility of the consecutive generations. Although it is desirable to minimize the number of predictors in the GA model in order to minimize costs and simplify experimental procedures, the prediction performance of the GA was poor and unstable when the number of predictors selected was <20. Therefore, the number of predictors in the GA was set at 20 to achieve the highest possible prediction accuracy with the minimum number of predictors.
The present invention is to identify markers (predictors) based on NCI-60 cell lines and then to validate their performance in patient cohorts. Therefore, we developed the prediction model for radiosensitivity based on gene expression data. Furthermore, it is not difficult to get lymphocytes specimen from patients with glioblastoma; therefore peripheral blood samples are suggested as the source for future applications. Third, the prediction outcomes were dichotomized into RS and RR groups by SVM algorithm based on the SF2 parameter.
Accordingly, in the present invention, the expression levels of both genes and long non-coding RNAs (lncRNAs) were used to build such a prediction model. Analysis of variance and Tukey's honest significant difference tests (P<0.001) were utilized in immortalized B cells (GSE26835) to identify differentially expressed genes and lncRNAs after irradiation. A total of 41 genes and lncRNAs associated with radiation exposure were revealed by a network analysis algorithm. To develop a predictive model for radiosensitivity, the expression profiles of NCI-60 cell lines along, with their radiation parameters, were analyzed. A genetic algorithm was proposed to identify 20 markers, and the support vector machine algorithm was used to evaluate their prediction performance. The model was applied to 2 datasets of glioblastoma, The Cancer Genome Atlas and GSE16011, and significantly better survival was observed in patients with greater predicted radiosensitivity
While the invention has explained in relation to its preferred embodiments, it is well understand that various modifications thereof will become apparent to those skilled in the art upon reading the specification. Therefore, the invention disclosed herein intended to cover such modifications as fall within the scope of the appended claims.