Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Natural Cubic Spline Regression Modeling Followed by Dynamic Network Reconstruction for the Identification of Radiation-Sensitivity Gene Association Networks from Time-Course Transcriptome Data

  • Agata Michna,

    Affiliation Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, Neuherberg, Germany

  • Herbert Braselmann,

    Affiliations Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, Neuherberg, Germany, Clinical Cooperation Group "Personalized Radiotherapy in Head and Neck Cancer", Helmholtz Zentrum München, Neuherberg, Germany

  • Martin Selmansberger,

    Affiliation Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, Neuherberg, Germany

  • Anne Dietz,

    Affiliation Department of Radiation Protection and Health, Federal Office for Radiation Protection, Neuherberg, Germany

  • Julia Hess,

    Affiliations Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, Neuherberg, Germany, Clinical Cooperation Group "Personalized Radiotherapy in Head and Neck Cancer", Helmholtz Zentrum München, Neuherberg, Germany

  • Maria Gomolka,

    Affiliation Department of Radiation Protection and Health, Federal Office for Radiation Protection, Neuherberg, Germany

  • Sabine Hornhardt,

    Affiliation Department of Radiation Protection and Health, Federal Office for Radiation Protection, Neuherberg, Germany

  • Nils Blüthgen,

    Affiliation Institute of Pathology, Charité—Universitätsmedizin Berlin, Berlin, Germany

  • Horst Zitzelsberger,

    Affiliations Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, Neuherberg, Germany, Clinical Cooperation Group "Personalized Radiotherapy in Head and Neck Cancer", Helmholtz Zentrum München, Neuherberg, Germany

  • Kristian Unger

    unger@helmholtz-muenchen.de

    Affiliations Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, Neuherberg, Germany, Clinical Cooperation Group "Personalized Radiotherapy in Head and Neck Cancer", Helmholtz Zentrum München, Neuherberg, Germany

Abstract

Gene expression time-course experiments allow to study the dynamics of transcriptomic changes in cells exposed to different stimuli. However, most approaches for the reconstruction of gene association networks (GANs) do not propose prior-selection approaches tailored to time-course transcriptome data. Here, we present a workflow for the identification of GANs from time-course data using prior selection of genes differentially expressed over time identified by natural cubic spline regression modeling (NCSRM). The workflow comprises three major steps: 1) the identification of differentially expressed genes from time-course expression data by employing NCSRM, 2) the use of regularized dynamic partial correlation as implemented in GeneNet to infer GANs from differentially expressed genes and 3) the identification and functional characterization of the key nodes in the reconstructed networks. The approach was applied on a time-resolved transcriptome data set of radiation-perturbed cell culture models of non-tumor cells with normal and increased radiation sensitivity. NCSRM detected significantly more genes than another commonly used method for time-course transcriptome analysis (BETR). While most genes detected with BETR were also detected with NCSRM the false-detection rate of NCSRM was low (3%). The GANs reconstructed from genes detected with NCSRM showed a better overlap with the interactome network Reactome compared to GANs derived from BETR detected genes. After exposure to 1 Gy the normal sensitive cells showed only sparse response compared to cells with increased sensitivity, which exhibited a strong response mainly of genes related to the senescence pathway. After exposure to 10 Gy the response of the normal sensitive cells was mainly associated with senescence and that of cells with increased sensitivity with apoptosis. We discuss these results in a clinical context and underline the impact of senescence-associated pathways in acute radiation response of normal cells. The workflow of this novel approach is implemented in the open-source Bioconductor R-package splineTimeR.

Introduction

In general terms, the expression of genes can be studied from a static or temporal point of view. Static microarray experiments allow measuring gene expression responses only at one single time point. Therefore, data obtained from those experiments can be considered as more or less randomly taken snapshots of the molecular phenotype of a cell. However, biological processes are dynamic and thus, the expression of a gene is a function of time [1]. To be able to understand and model the dynamic behavior and association of genes, it is important to study gene expression patterns over time.

However, compared to static microarray data, the analysis of time-course data introduces a number of new challenges. First, the experimental costs for the generation of data as well as the computational cost increases with the increase in the number of introduced time points. Second, hidden correlation caused by co-expression of genes makes the data linearly dependent [2]. Finally, one has to be aware of additional correlations existing between neighboring time points clearly revealed in published gene expression profiles [3].

Several different algorithms have been suggested to analyze gene time-course microarray data with regard to differential expression in two or more biological groups (e.g. exposed to radiation vs. non-exposed) [47]. Nevertheless solitary identification of differentially expressed genes does not help to determine the molecular mechanisms in the investigated biological groups. Therefore, it is not only important to know differentially expressed genes per se, but also how those genes interact and regulate each other in order to determine specifically deregulated molecular networks.

Currently, many different algorithms including cluster analysis [813] and supervised classification [1416] are used to identify relationships between genes. However, both of these methods suffer from serious limitations. First, the timing information of the measurements is not incorporated and, therefore, the intrinsic temporal structure of the time-course data is neglected. Second, the available standard clustering and classification methods are not designed to measure statistical significance of the results based on a statistical hypothesis test. By nature of these methods, clusters or classes of genes with similar expression patterns will always be identified but they do not provide a measure of how reliable this information is. For this reason, we preferred usage of a dynamic network modeling approach that allows delineation of relationships between genes along with providing statistical significance for these relationships.

The aim of the present study was to identify and compare signaling pathways involved in the radiation responses of normal cells differing in their radiation sensitivity that could be used to modulate cell sensitivity to ionizing radiation. For this, we propose an approach that combines the detection of genes differentially expressed over time based on statistics determined by natural cubic spline regression (NCSRM) [17] followed by dynamic gene association network (GAN) reconstruction based on a regularized dynamic partial correlation as implemented in the GeneNet R-package [18].

Most exploratory gene expression studies focus only on the identification of differentially expressed genes by treating them as independent events and do not seek to study the interplay of identified genes. This makes it difficult to tell which genes are part of the interaction network causal of the studied phenotype and which are the most “important” with regard to the context of the investigation. The herein present approach combines the identification of differentially expressed genes and reconstruction of possible associations between them. Further analysis of identified GANs then allows hypothesizing which genes may play a crucial role in the investigated processes. This should markedly increase the likelihood to find meaningful results from an initial observation and help to understand the underlying molecular mechanisms. We applied our workflow on time-course transcriptome data of two normal and well-characterized lymphoblastoid cell lines with normal (20037–200) and increased radiation sensitivity (4060–200), in order to identify molecular mechanisms and potential key players responsible for different radiation responses [19, 20]. Our exploratory approach provides novel and informative insights in the biology of radiation sensitivity of non-tumor cells after exposure to ionizing radiation with regard to the identified signaling pathways and their key drivers. Moreover, we could demonstrate that spline regression in differential gene expression analysis for the purpose of prior selection in gene-association network reconstruction outperforms another commonly used existing approach for time-course gene expression analysis.

Results

The schematic workflow of the presented novel approach for time-course gene expression data analysis is presented in Fig 1.

thumbnail
Fig 1. Schematic workflow of the analysis of gene expression time-course data.

Samples were collected 0.25, 0.5, 1, 2, 4, 8 and 24 hours after sham or actual irradiation. Transcriptional profiling was performed using Agilent gene expression microarrays and comprises three major steps: the identification of differentially expressed genes from time-course expression data by employing a natural cubic spline regression model; the use of regularized dynamic partial correlation method to infer gene associations networks from differentially expressed genes and the topological identification and functional characterization of the key nodes in the reconstructed networks.

https://doi.org/10.1371/journal.pone.0160791.g001

Identification of ionizing radiation-responsive genes using NCSRM method

A fraction of the probes was removed due to low expression levels, with not detectable signal intensities as described in [21]. Table 1 shows the number of probes remained after quality filtering from the total number of 25220 unique probes representing HGNC annotated genes. Differential analysis was performed relative to the corresponding sham irradiated cells as a reference. In general, more genes were detected as differentially expressed in the cells with increased radiation sensitivity compared to cells with normal radiation sensitivity after each dose of gamma irradiation (Table 1). The most prominent difference was observed when comparing the responses after 1 Gy irradiation. In the cells with increased radiation sensitivity 2335 genes showed differential expression compared to only seven genes in cells with normal radiation sensitivity. We observed the same trend after irradiation with 10 Gy where the cells with increased sensitivity showed 6019 and the normal sensitive cells 3892 differentially expressed genes.

thumbnail
Table 1. Number of detected and differentially expressed genes for each dose and cell lines for NCSRM and BETR methods.

https://doi.org/10.1371/journal.pone.0160791.t001

Pathway enrichment analysis of NCSRM identified genes

Pathway enrichment analysis was performed on differentially expressed genes to identify over-represented biological pathways. The analysis on genes identified with NCSRM revealed 634 and 964 significantly enriched pathways for the cells with increased radiation sensitivity after 1 Gy and 10 Gy irradiation dose, respectively and 758 pathways for the normal sensitive cell line after 10 Gy irradiation. For the seven differentially expressed genes (i.e. FDXR, BBC3, VWCE, PHLDA3, SCARF2, HIST1H4C, PCNA) of the cell line with normal radiation sensitivity after 1 Gy dose of irradiation we did not find any significantly enriched pathways. A summary of the pathway enrichment results can be found in S2 Table.

Gene association network reconstruction

None of the edge probabilities calculated for the seven differentially expressed genes in the cell line with normal radiation sensitivity after 1 Gy irradiation exceeded the considered significance threshold and hence no network was obtained. For the remaining conditions we were able to obtain association networks as presented in Table 2. Obtained networks are provided as igraph R-objects in the supplementary data (S1 File). The graph densities for all resulting networks were in the same range as the density of the Reactome interaction network (Table 2).

thumbnail
Table 2. Number of genes subjected to GAN reconstruction and properties of resulted GANs.

https://doi.org/10.1371/journal.pone.0160791.t002

Identification and functional characterization of the most important genes in the reconstructed association networks

The combined topological centrality measure was used to characterize the biological importance of nodes (genes) in the reconstructed association networks. The 5% of the highest ranked genes listed in supplementary S3 Table were mapped to Reactome pathways in order to further evaluate their biological roles. The top 10 most relevant pathways according to the FDR values are shown in Table 3. For the cell line with increased radiation sensitivity after irradiation with 1 Gy and for the normal sensitive cell line after 10 Gy the induction of pathways associated with senescence response was detected. For the cell line with increased radiation sensitivity after 10 Gy of irradiation we mostly observed pathways associated with apoptosis. All pathways are listed in supplementary S4 Table.

thumbnail
Table 3. Comparison of NCSRM and BETR methods with respect to the top 10 pathways after mapping of 5% highest ranked genes from the reconstructed gene association networks.

https://doi.org/10.1371/journal.pone.0160791.t003

False detected differentially expressed genes between technical replicates

In order to assess the false positive rate, the spline regression based differential analyses between technical replicates of each treatment conditions and cell lines were performed. Here, we can state that the null-hypothesis of no differential expression is true for all genes. Then the q*-level of 0.05 for Benjamini-Hochberg method controls also the FWER at alpha-level equal to 0.05 (type I error) [22]. For all compared technical replicates not more than 3% rejections of null hypothesis were detected, which is in good accordance to the expected or nominal type I error.

Evaluation of spline regression model in comparison to BETR method

Table 1 compares the numbers of differentially expressed genes obtained from both methods applied on the same gene expression data set and FDR thresholds. For almost all treatment conditions the BETR method detected less differentially expressed genes in comparison to NCSRM. Only for the normal cell line after irradiation with 1 Gy BETR identified 12 genes whereas NCSRM identified only 7 genes. As a consequence of the lower numbers of detected differentially expressed genes with BETR, the obtained networks are smaller than those obtained after spline regression. The detailed comparison results including numbers of detected differentially expressed genes and the sizes of reconstructed association networks are presented in the Table 2. The lists of differentially expressed genes obtained with the two methods are shown in supplementary S1 Table. The top 10 pathways to which the 5% of the most important genes in the reconstructed association networks where mapped to are shown in Table 3. With NCSRM we were not only able to detect almost all genes that were detected also by BETR (Table 1), but also an additional set of genes resulting in almost twice the number of genes compared to BETR. Nevertheless, the top 5% hub genes of the networks derived from the differentially expressed genes defined by BETR were associated with similar biological processes as those from the spline differential expression analysis derived networks. The numbers and names of overlapping hub genes in the GANs are presented in Table 4 and in supplementary S3 Table, respectively.

thumbnail
Table 4. Comparison of hub genes in networks resulting from different methods.

https://doi.org/10.1371/journal.pone.0160791.t004

Evaluation of reconstructed networks

The evaluation of the two networks derived after 1 Gy irradiation of the cell line with increased sensitivity showed that the network reconstructed with the differentially expressed genes determined using BETR did not contain significantly more common edges than random networks (p = 0.529), whereas the network reconstructed with the differentially expressed genes determined by NCSRM did (p = 0.048). The networks derived after 10 Gy irradiation of the cell line with increased sensitivity and 10 Gy irradiation of the normal sensitive cell line contained significantly more edges that were common with the Reactome network compared to random networks for both methods.

Discussion

The success of tumor radiation therapy predominantly depends on the total applied radiation dose, but also on the tolerance of the tumor surrounding normal tissues to radiation. Toxicity towards radiation, which greatly varies on an individual level due to inherited susceptibility, is one of the most important limiting factors for dose escalation in radiooncology treatment [23, 24]. To account for radiation sensitivity of normal tissue in personalized treatment approaches the underlying molecular mechanisms need to be thoroughly understood in order to identify molecular targets for the modulation of radiation sensitivity and molecular markers for the stratification of patients with different intrinsic radiation sensitivity. In the present study we identified significantly differentially expressed genes over time between the radiation-treated group and the control group to be used as prior genes for GAN reconstruction. Two doses of gamma irradiation were used to characterize the differences in radiation response of the two lymphoblastoid cell lines with known differences in radiation sensitivity. The dose of 10 Gy was selected following the fact that the same dose has been applied in a previous research project examining the radiation sensitivity of the same lymphoblastoid cell lines analyzed in the study at hand [20]. The dose of 1 Gy reflects the dose that is delivered as part of the so called “low-dose bath” to the tumor-surrounding tissue during the radiotherapy of the tumors [25].

Here, we conducted time-resolved transcriptome analysis of radiation-perturbed cell culture models of non-tumor cells with normal and with increased radiation sensitivity in order to work out the molecular phenotype of radiation sensitivity in normal cells. Moreover, we present an innovative approach for the identification of GANs from time-course perturbation transcriptome data. The approach comprises three major steps: 1) the identification of differentially expressed genes from time-course gene expression data by employing a natural cubic spline regression model (NCSRM); 2) the use of a regularized dynamic partial correlation method to infer gene associations network from differentially expressed genes; 3) the identification and functional characterization of the key nodes (hubs) in the reconstructed gene dependency network (Fig 1).

Our proposed method for the detection of differentially expressed genes over time is based on NCSRM with a small number of basis functions. A relatively low number of basis functions generally results in a good fit of data and, at the same time, reduces the complexity of the fitted models. Treating time in the model as a continuous variable, a non-linear behavior of gene expressions was approximated by spline curves fitted to the experimental time-course data. Considering temporal changes in gene expression as continuous curves and not as single time points greatly decreases the dimensionality of the data and thereby decreases computational cost. In addition, the proposed NCRSM does not require identical sampling time points for the compared treatment conditions. Furthermore, no biological replicates are needed. Therefore, the method is applicable to data generated according to a tailored time-course differential expression study design and to data that were not specifically generated for time-course differential expression analysis, e.g. existing/previously generated data from clinical samples. Thus, the adaption of the method to differential expression analysis comprises the potential to reanalyze existing data, address new questions in silico and thereby potentially add new or additional value to existing data. Incomplete time-course data, e.g. due to the exclusion of samples for technical reasons, that often create major problems for the estimation of the model, are also suitable for fitting the spline regression model as long as enough data points remain in the data set. This is especially valuable when data on certain time points, derived from a very limited sample source, have been excluded from a time-course data set and cannot be repeatedly generated.

Since gene expression is not only dynamic in the treatment group but also in the control group, the inclusion of the time-course control data greatly improves the ability to detect truly differentially expressed genes, as the gene expression values are not referred to a single time point with static gene expression levels only. Comparing a treatment group to time point zero does not provide a proper control over the entire time-course, although it is widely practiced [2628]. The proposed workflow is implemented in an open-source R-package splineTimeR and is available through Bioconductor (https://www.bioconductor.org).

Amongst a panel, the two lymphoblastoid cell lines that were different with regard to radiation sensitivity after irradiation with 10 Gy [20], also responded differently with regard to the quantity of differentially expressed genes. Interestingly, cells with normal radiation sensitivity barely responded to 1 Gy irradiation at the transcriptome level. Only seven genes (FDXR, BBC3, VWCE, PHLDA3, SCARF2, HIST1H4C, PCNA) were identified as differentially expressed, whereas for the cell line with increased sensitivity 2335 differentially expressed genes were detected after exposure to the same dose. A similar behavior was observed for those two cell lines after irradiation with 10 Gy. We detected 6019 and 3892 genes as differentially expressed in the sensitive and normal cell lines, respectively (Table 2). Those results are in a good agreement with the previous proteomic study where more differentially expressed proteins were detected for the same sensitive cell line compare to the cell line with normal radiation sensitivity 24 hours after irradiation with 10 Gy [29]. Thus, for both applied doses, the radiation sensitive cells exhibited much more pronounced transcriptional response compared to the cells with normal radiation sensitivity and thereby underlines the expected radiation response of those two cell lines.

Concerning qualitative differences in the transcriptomic response of normal sensitive cells and cells with increased sensitivity after treatment with 1 Gy and 10 Gy pathway enrichment analysis was performed. Differentially expressed genes identified for all considered treatment conditions except for the normal sensitive cells after exposure to 1 Gy radiation showed statistically significant enrichment of pathways. Most of which were in agreement with known radiation responses such as DNA repair, cell cycle regulation, oxidative stress response or pathways related to apoptosis (S2 Table) [3032]. Therefore, the pathway enrichment analysis results suggest plausibility of generated data and, more importantly, underline the meaningfulness of our suggested approach based on cubic spline regression for differential gene expression analysis of time-course data. However, differential expression analysis alone followed by pathway enrichment analysis does not provide any mechanistic insights. For this reason we performed GAN reconstruction using identified differentially expressed genes. Based on the assumption that the expression levels of functionally related genes are highly correlated, partial correlation was used for GAN reconstruction. In simple correlation, the strength of the linear relationship between two genes is measured, without taking into account that those genes may be actually influenced by other genes. Partial correlation eliminates the influence of other genes when one particular relationship between pair of genes is considered. Network reconstruction was performed separately for the cell line with increased radiation sensitivity after 1 Gy and 10 Gy and for the cell line with normal radiation sensitivity after 10 Gy of radiation dose. Due to the sparseness of the set of genes differentially expressed after irradiation of the normal-sensitive cell line with 1 Gy, no GAN was obtained.

Subsequently, we identified the network hubs (i.e. most important genes) of the GANs by combining three network centrality measures: degree, closeness and shortest path betweenness [33]. Combining different centrality measures is a widely used approach to identify nodes that are likely to control the network [34]. Also, this approach allows identification of nodes that are connected to the central nodes at the same time which can be informative for the interpretation of the whole GAN or single modules making up the network [33, 34].

Identification of key pathways associated with radiation sensitivity

In order to get functional insights into the reconstructed GANs the 5% top important nodes were identified after a ranking with the combined centrality measure and mapped to the pathways from the interactome database Reactome [35]. The obtained results revealed different pathways considered as the most important in cells with different radiation sensitivity after different doses of ionizing radiation. For the radiation sensitive cell line 4060–200 and 1 Gy irradiation, we mainly detected pathways associated with senescence (Table 3).

A different outcome was observed after irradiation with 10 Gy. For the radiation sensitive cells three out of the ten top pathways were linked to apoptotic processes with the genes BBC3, BCL2, TP53 as key players, whereas for the normal sensitive cell line we mainly observed the induction of senescence related pathways. This indicates that different doses are necessary to induce a similar response in the two cell lines. The activation of senescence genes is a damage response mechanism, which stably arrests proliferating cells and protects them from apoptotic cell death [36]. Together with the senescence pathway we observed increased levels of chemokine, cytokine and interleukin genes that are known to activate an immune response and signal transduction pathways in response to irradiation.

Although the senescence-associated pathways were not seen as the most important ones for the treatment condition 10 Gy/increased sensitivity, they were significantly enriched in the GANs of the three conditions 1 Gy/increased sensitivity, 10 Gy/ increased sensitivity and 10 Gy/normal sensitivity. All differentially expressed genes that related to senescence-associated pathways are shown in supplementary S5 Table. The observation that cells with increased radiation sensitivity compared to cells with normal sensitivity, become senescent after exposure to doses in the range of 1 Gy, rises the question whether this has a positive or negative influence on the tumor therapy. On the one hand side, senescent cell may secret the so-called SASP (“senescence-associated secretory phenotype“) factors, including growth factors, chemokines and cytokines, which participate in intercellular signaling leading to the attraction of immune cells to the tumor location that, in turn, eliminate the tumor cells and, thereby, positively contribute to the tumor therapy [37, 38]. On the other hand side, senescent cells and the SASP are reported to promote proliferation, survival, invasion and migration of neighboring cells by the release of pro-inflammatory cytokines leading to sustained inflammation [36]. In this way senescence cells can damage their local environment and stimulate angiogenesis and tumor progression [39, 40]. Besides, there are some evidences that the induction of senescence in surrounding normal tissue may lead to an increased radio-tolerance or even radioresistance of the tumor and is, therefore, not desirable and negatively influences the tumor radiotherapy [41]. Thus, it might be beneficial to block senescence in order to prevent the radio-hyposensibilization of tumor cells. Therefore, we suggest a detailed investigation of the consequences of senescent non-tumor cells with the aim to improve the radiotherapy of tumors in radiosensitive patients.

Identification of senescence associated genes involved in cell radiation responses

CDKN1A gene was identified as one of the most important key players linked to the identified senescence associated pathways for both 1 Gy/sensitive and 10 Gy/normal treatment conditions. For both conditions the expression of the CDKN1A was up-regulated for all considered time points. CDKN1A is a well-known damage response gene for which aberrant transcriptional response has been associated with abnormal sensitivity to ionizing radiation [42, 43]. The study by Badie et al. (2008) has shown that a subgroup of breast cancer patients, who developed severe reactions to radiation therapy, could be identified by aberrant overexpression of CDKN1 in peripheral blood lymphocytes [43].

LMNB1 is another genes we identified as a response hub gene after irradiation of sensitive cell line with 1 Gy radiation dose that is associated with senescence. Although the LMNB1 gene was not identified as hub gene in the GAN of the 10 Gy/normal treatment condition, it was still differentially expressed. For both treatment conditions we observed significant downregulation of this gene 24 hours after irradiation. Shah et al (2013) has suggested that downregulation of LMNB1 in senescence is a key trigger of chromatin changes affecting gene expression [44]. In fact also in our data we observed strong downregulation of a group of histone genes associated with senescence (S5 Table) for the treatment conditions 1 Gy/increased sensitivity and 10 Gy/normal sensitivity. Furthermore, Lee et al. (2012) has shown that histone protein modification may have an impact on the radiation sensitivity of a tissue [45]. Moreover, evidence has been provided that mutation of LMNA can cause increased sensitivity to ionizing radiation [46], however, to our knowledge there are no data showing the role of LMNB gene in the context of radiation sensitivity.

Another potential therapeutic candidate associated with senescence that was identified for the 10 Gy/normal sensitivity treatment condition was MRE11A for which cell culture data suggest that treatment of cells with Mre11 siRNA increases radiation sensitivity and reduces heat-induced radiosensitization [47, 48]. However, the clinical applicability of MRE11, has not been confirmed [49].

Assessment of the false positive rate and validation of the NCSRM method

The spline regression based differential analyses between technical replicates were performed in order to estimate the extent of random fluctuations of gene expression values. The detected 3% rejections of the overall null hypothesis of no differential gene expression are in accordance with the alpha-level of 5% of the familywise error rate (FWER) and can be considered as false positives. On the other hand, it shows that type I error, due to technical variation, is covered by the model and test assumptions (moderated F-test, [50]) so that it was not necessary to include an extra parameter for technical replicates into the model.

In order to validate the previously mentioned biological results using NCSRM, we performed the differential expression analysis with another established method for time-course data analysis called BETR (Bayesian Estimation of Temporal Regulation) [6]. The number of genes detected by BETR was considerably lower compared to NCRSM (Table 1), however the majority of which were also detected with NCSRM (S1 Table). This is in line with the calculations on the false positive rates that have been conducted on the simulated data presented in the BETR study. In an analysis of the simulated data set, 65% of truly differentially expressed genes have been identified after accepting a false positive rate of 5% [6]. This means that a substantial proportion of differentially expressed genes remained undetected, which is likely to be also the case for the herein analyzed data with BETR. Although the numbers of differentially expressed genes and genes remained in the reconstructed networks greatly differ (Table 1), the qualitative results are well comparable (Table 3). For all treatment conditions where for which we were able to reconstruct GANs, we observed a great overlap of pathways where the 5% of hub genes were mapped to (Table 3). The detection of a higher number of differentially expressed genes with NCSRM resulted in larger GANs with additional information compared to the smaller GANs that were reconstructed on the basis of genes detected with BETR. This is underlined by the results of the conducted evaluation of GANs. Except one network based on the differentially expressed genes using BETR, all investigated networks consist significantly more common edges with the Reactome reference network compared to random networks with identical network topology and genes. This shows that the additionally detected genes with NCSRM add additional information rather than adding false positives or noise to the set of differentially expressed genes. Moreover the spline regression method is much more flexible and allows for more freedom during the data collection process. As already mentioned, NCSRM does not require the same sampling time for treated and control groups and can easily deal with incomplete data, whereas BETR method is not able to overcome or bypass those limitations. Thus, NCSRM is very robust against the frequently occurring shortcomings in study design and subsequent data generation occurring in life sciences.

Conclusion

Prospectively, we suggest and plan a detailed in silico and in vitro analysis of the interactions in the proposed gene association networks in order to add meaningful knowledge to the mechanism of radiosensitivity at the experimental level. This novel knowledge has the potential to improve cancer radiation therapy by preventing or lowering the acute responses of normal cells resulting from radiation therapy. The results add novel information to the understanding of mechanisms that are involved in the radiation response of human cells, with the potential to improve tumor radiotherapy. Besides, the presented workflow is not limited to presented study only, but may be applied in other special fields with different biological questions to be addressed.

The software is provided as R-package “splineTimeR” and freely available via the Bioconductor project at http://www.bioconductor.org.

Material and Methods

Cell culture

Experiments were conducted with two monoclonal lymphoblastoid Epstein-Barr virus-immortalized cell lines (LCL) obtained from young lung cancer patients of the LUCY study (LUng Cancer in Young) that differ in radiosensitivity, as tested with Trypan Blue and WST-1 assays [19, 20]. The non-cancer cell lines LCL 4060–200 with increased radiation sensitivity and LCL 20037–200 with normal radiation sensitivity were cultured at 37°C/5% CO2 in RPMI 1640 medium (Biochrom) supplemented with 10% fetal calf serum (FCS; PAA). Mycoplasma contamination was routinely tested using luminescence-based assays (MycoAlert, Lonza).

Irradiation and sample preparation

The cells were seeded in 75 cm2 flasks at a concentration of 0.5 x 106 cells/ml in a total volume of 60 ml. Exponentially growing cells were irradiated with sham, 1 Gy and 10 Gy of gamma-irradiation (137Cs-source HWM-D 2000, Markdorf, Germany) at a dose rate of 0.49 Gy/min. Samples were collected 0.25, 0.5, 1, 2, 4, 8 and 24 hours after sham or actual irradiation. Between the time of collection cells were kept in the incubator. Collected cells were washed with PBS and frozen at -80°C. Total RNA was isolated from frozen cell pellets obtained from two independent experiments using the AllPrep DNA/RNA/miRNA Universal Kit (Qiagen) including an DNase digestion step, according to the manufacturer's protocol. The concentration of RNA was quantified with a Qubit 2.0 Fluorometer (Life Technologies), and integrity was determined using a Bioanalyzer 2100 (Agilent Technologies). RNA samples with a RNA integrity number (RIN) greater than 7 indicated sufficient quality to be used in subsequent RNA microarray analysis.

Gene expression profiling

Transcriptional profiling was performed using SurePrint G3 Human Gene Expression 8x60k V2 microarrays (Agilent Technologies, AMADID 39494) according to the manufacturer’s protocol. 75 ng of total RNA was used in labeling using the Low Input Quick Amp Labeling Kit (one-color, Agilent Technologies). Raw gene expression data were extracted as text files with the Feature Extraction software 11.0.1.1 (Agilent Technologies). The expression microarray data were uploaded to ArrayExpress (www.ebi.ac.uk/arrayexpress/) and the data set is available under the accession number E-MTAB-4829. All data analysis was conducted using the R statistical platform (version 3.2.2, www.r-project.org) [51]. Data quality assessment, filtering, preprocessing, normalization, batch correction based on nucleic acid labeling batches and data analyses were carried out with the Bioconductor R-packages limma, Agi4x44PreProcess and the ComBat function of the sva R-package [4, 21, 52]. All quality control, filtering, preprocessing and normalization thresholds were set to the same values as suggested in Agi4x44PreProcess R-package user guide [21]. Only HGNC annotated genes were used in the analysis. For multiple microarray probes representing the same gene the optimal probe was selected according to the Megablast score of probe sequences against the human reference sequence (http://www.ncbi.nlm.nih.gov/refseq/) [53]. If the resulted score was equal for two or more probes, the probe with the lowest differential gene expression FDR value was kept for further analyses since only one expression value per gene was allowed in subsequent GAN reconstruction analysis.

Spline regression model for two-way experimental design

A natural cubic spline regression model (NCSRM) with three degrees of freedom for an experimental two-way design with one treatment factor and time as a continuous variable was fitted to the experimental time-course data. The mathematical model is defined by the following eq (1): where b0, b1, …, bm are the spline coefficients in the control group and d0, d1, …, dm are differential spline coefficients between the control and the irradiated group. B1(t-t0), B2(t-t0), …, Bm(t-t0) are the spline base functions and t0 is the time of the first measurement. For x = 0, y = ycontrol and for x = 1, y = yirradiated. For three degrees of freedom (df = 3), m = 3.

Depending on the number of degrees of freedom, two boundary knots and df-1 interior knots are specified. The interior knots were chosen at values corresponding to equally sized quantiles of the sampling time from both compared groups. For example, for df = 3 interior knots correspond to the 0.33- and 0.66-quantiles. The spline function is cubic on each defined by knots intervals, continuous at each knot and has continuous derivatives of first and second orders.

Time-course differential gene expression analysis

The time-course differential gene expression analyses were conducted between irradiated and control cells (sham-irradiated). Analyses were performed on the normalized gene expression data using NCSRM with three degrees of freedom. The splines were fitted to the real time-course expression data for each gene separately according to eq (1). The example of spline regression model fitted to the measured time-course data for one selected gene is shown on the Fig 2.

thumbnail
Fig 2. Example of fitted spline regression models.

The plot shows spline regression models fitted to the measured time-course expression data of an arbitrary chosen gene (BBC3). The blue line represents the fitted model for the control (0 Gy) and read line that for the irradiated group (1 Gy). Blue and red dots represent the measured expression levels of the biological replicates. Vertical lines represent the endpoints and interior knots correspond to the 0.33- and 0.66-quantiles.

https://doi.org/10.1371/journal.pone.0160791.g002

Time dependent differential expression of a gene between the irradiated and corresponding control cells was determined by the application of empirical Bayes moderated F-statistics [50] on the differential coefficients values in eq (1). In order to account for the multiple-testing error, corresponding p-values were adjusted by the Benjamini-Hochberg method for false discovery [22]. Genes with an adjusted p-value (FDR, false discovery rate) lower than 0.05 were considered as differentially expressed and associated with radiation response.

Assessment of the false positive rate of the NCSRM

Additionally, in order to assess the false positive rate (statistical type I error, also called familywise error rate or FWER) we applied differential gene expression analysis using NCSRM between two technical replicates for all treatment groups. Because only two technical replicates were generated for each time point and treatment, we could not use the same approach to assess the technical variability for the BETR method, as it requires at least two replicates in each compared groups.

Gene association network reconstruction from prior selected differentially expressed genes

Differentially expressed genes were subjected to gene association network reconstruction from time-course data using a regularized dynamic partial correlation method [54]. Pairwise relationships between genes over time were inferred based on a dynamic Bayesian network model with shrinkage estimation of covariance matrices as implemented in the GeneNet R-package available from CRAN [18]. Analyses were conducted with a posterior probability of 0.95 for each potential edge. Edge directions were not considered. In order to assess the complexity of the resulting networks, the density of each network was compared to the density of the Reactome functional interaction network [35, 55].

Identification of important nodes in the network

Graph topological analyses based on centrality measures were applied in order to determine the importance of each node in the reconstructed association networks [56]. Three most commonly used centrality measures: degree, shortest path betweenness and closeness were combined into one cumulative centrality measure [34]. For each gene the three centrality values where ranked. The consensus centrality measure for each node was defined as the mean of the three independent centrality ranks. Combining centrality measures supports the identification of the nodes that are central themselves and also connected to direct central nodes, which demonstrates strategic positions for controlling the network.

Pathway enrichment analysis

The Reactome pathway database was used to conduct the pathway enrichment analysis in order to further investigate the functions of the selected sets of differentially expressed genes [35]. Statistical significance of enriched pathways was determined by one-sided Fisher's exact test. The resulting p-values were adjusted for FDR using the Benjamini-Hochberg method. Pathways with FDR<0.05 were considered statistically significant and pathways were ranked according to ascending FDRs.

Evaluation of NCSRM approach

Since we decided to use the set of genes that appeared to be differentially expressed we assessed the performance of the herein used NCSRM approach in comparison to the BETR approach implemented in the R/Bioconductor package betr [6]. BETR is a well-established algorithm that has been previously compared to limma, MB-statistic and EDGE methods and showed the best performance [6]. The results of spline and BETR methods were compared using the same initial microarray gene expression data set. The probabilities of each gene to be differentially expressed obtained with BETR method, were transformed to p-values as described in the original paper. Genes were considered significantly differentially expressed if the Benjamini-Hochberg adjusted p-value was lower than 0.05 (FDR<0.05). This transformation allowed us to compare the outcomes of both methods based on the FDR values for differential expression. The resulting differentially expressed genes using BETR were analyzed and subjected to network reconstruction as described above for the differentially expressed genes obtained using NCSRM. Outcomes of both obtained association networks were compared to each other and to the a priori known biological network provided by the Reactome database [35].

Evaluation of reconstructed gene association networks

In order to assess the quality of the de novo reconstructed gene association networks (GANs), we developed a novel method that compares the interactions in the reconstructed network to the experimentally validated interactions present in the Reactome interaction network. For this purpose we used the Reactome reference network, consisting of protein-protein interaction pairs stored in the Reactome database (http://www.reactome.org/pages/download-data/). For the comparison, sub-networks of reconstructed networks consisting only of genes overlapping with the Reactome network were built. The number of common edges between these two sub-networks was determined and referred to the total number of edges in the reconstructed network (percentage of common edges in the reconstructed network). Further, a permutation test was performed to assess whether the number of common edges in the reconstructed network was significantly higher than in randomized networks with the same genes. Random networks were generated by permutation of the node names in the network, while preserving the reconstructed sub-network topology. After each permutation (n = 1000) the number of common edges with the reference Reactome sub-network was determined. The reconstructed network was considered significantly better than random, if more than 90% of the random sub-networks contained lower numbers of edges common with the Reactome network than the reconstructed sub-network (p-value < 0.1). All networks reconstructed with the genes determined as differentially expressed from the herein presented spline regression method and the BETR method were evaluated.

Supporting Information

S1 File. Reconstructed gene association networks.

All obtained gene association networks are provided as R-objects of type igraph.

https://doi.org/10.1371/journal.pone.0160791.s001

(RDATA)

S1 Table. Lists of differentially expressed genes.

Table includes differentially expressed genes identified by spline regression and BETR methods. Additionally, a list of overlapping differentially expressed genes between both methods is included.

https://doi.org/10.1371/journal.pone.0160791.s002

(XLSX)

S2 Table. Lists of significantly enriched pathways using differentially expressed genes identified by spline regression method.

Four lists of significantly enriched pathways correspond to each used treatment condition. Lists include total numbers of known genes in the pathways, numbers of differentially expressed genes that belong to a single pathway (matches), percentages of differentially expressed genes in comparison to the total number of know genes in the pathway (% match), p-values, FDRs and names of pathways related differentially expressed genes.

https://doi.org/10.1371/journal.pone.0160791.s003

(XLSX)

S3 Table. Lists of 5% of most important genes identified by centrality measures.

Lists of 5% highest ranked genes from the reconstructed gene association networks using spline regression and BETR methods. Overlap represents common most important genes identified in networks from compared methods.

https://doi.org/10.1371/journal.pone.0160791.s004

(XLSX)

S4 Table. Lists of pathways after mapping of 5% highest ranked genes from the reconstructed gene association networks.

Lists include names of pathways together with names of mapped most important genes.

https://doi.org/10.1371/journal.pone.0160791.s005

(XLSX)

S5 Table. Significantly enriched senescence associated pathways with corresponding differentially expressed genes.

Table presents the names of significantly enriched (FDR<0.05) senescence associated pathways with corresponding differentially expressed genes for all treatment conditions.

https://doi.org/10.1371/journal.pone.0160791.s006

(XLSX)

Acknowledgments

We thank Aaron Selmaier from the Research Unit Radiation Cytogenetics for technical support. This study was supported by the Federal Ministry of Education and Research, ZiSS project, 02NUK024B.

Author Contributions

  1. Conceived and designed the experiments: AD AM KU HZ NB JH MG SH.
  2. Performed the experiments: AM AD.
  3. Analyzed the data: AM HB MS.
  4. Contributed reagents/materials/analysis tools: AM AD MS KU MG SH HZ.
  5. Wrote the paper: AM KU MS.
  6. Biological data interpretation: AM KU HZ JH MS MG SH NB. Designed the software used in analysis: AM HB. Revision of the manuscript: KU HZ JH MG SH NB AD HB. Final approval of the version to be published: AM HB MS AD JH MG SH NB HZ KU.

References

  1. 1. Bar-Joseph Z, Gitter A, Simon I. Studying and modelling dynamic biological processes using time-series gene expression data. Nature Reviews Genetics. 2012;13:552–64. pmid:22805708
  2. 2. Bandyopadhyay S, Bhattacharyya M. A biologically inspired measure for coexpression analysis. IEEE/ACM Transactions on Computational Biology and Bioinformatics. 2011;8(4):929–42. pmid:21566252
  3. 3. Yuan M, Kendziorski C. Hidden Markov models for microarray time course data in multiple biological conditions. Journal of the American Statistical Association. 2006;101(476):1323–32.
  4. 4. Smyth GK. Limma: linear models for microarray data. In: 'Bioinformatics and Computational Biology Solutions Using {R} and Bioconductor'. New York: Springer, p. 397–420; 2005.
  5. 5. Leek JT, Monsen E, Dabney AR, Storey JD. EDGE: extraction and analysis of differential gene expression. Bioinformatics. 2006;22(4):507–8. pmid:16357033
  6. 6. Aryee MJ, Gutiérrez-Pabello JA, Kramnik I, Maiti T, Quackenbush J. An improved empirical bayes approach to estimating differential gene expression in microarray time-course data: BETR (Bayesian Estimation of Temporal Regulation). BMC Bioinformatics. 2009;10(409).
  7. 7. Conesa A, Nueda MJ, Ferrer A, Talón M. maSigPro: a method to identify significantly differential expression profiles in time-course microarray experiments. Bioinformatics. 2006;22(9):1096–102. pmid:16481333
  8. 8. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006;7(191).
  9. 9. Schliep A, Steinhoff C, Schönhuth A. Robust inference of groups in gene expression time-courses using mixtures of HMMs. Bioinformatics. 2004;20:i283–i9. pmid:15262810
  10. 10. Magni P, Ferrazzi F, Sacchi L, Bellazzi R. TimeClust: a clustering tool for gene expression time series. Bioinformatics. 2008;24(3):430–2. pmid:18065427
  11. 11. Sivriver J, Habib N, Friedman N. An integrative clustering and modeling algorithm for dynamical gene expression data. Bioinformatics. 2011;27(13): i392–i400. pmid:21685097
  12. 12. Sinha A, Markatou M. A Platform for Processing Expression of Short Time Series (PESTS). BMC Bioinformatics. 2011;12(13).
  13. 13. Ramoni MF, Sebastiani P, Kohane IS. Cluster analysis of gene expression dynamics. Proc Natl Acad Sci USA. 2002;99:9121–6. pmid:12082179
  14. 14. Lin T, Kaminski N, Bar-Joseph Z. Alignment and classification of time series gene expression in clinical studies. Bioinformatics. 2008;24:i147–i55. pmid:18586707
  15. 15. Costa IG, Schönhuth A, Hafemeister C, Schliep A. Constrained mixture estimation for analysis and robust classification of clinical time series. Bioinformatics. 2009;25:i6–i14. pmid:19478017
  16. 16. Hafemeister C, Costa IG, Schönhuth A, Schliep A. Classifying short gene expression time-courses with Bayesian estimation of piecewise constant functions. Bioinformatics. 2011;27(7):946–52. pmid:21266444
  17. 17. Storey JD, Xiao W, Leek JT, Tompkins RG, Davis RW. Significance analysis of time course microarray experiments. PNAS. 2005;102(36):12837–42. pmid:16141318
  18. 18. Opgen-Rhein R, Strimmer K. From correlation to causation networks: a simple approximate learning algorithm and its application to high-dimensional plant gene expression data. BMC Syst Biol. 2007;1:37. pmid:17683609
  19. 19. Rosenberger A, Rossler U, Hornhardt S, Sauter W, Bickeboller H, Wichmann HE, et al. Validation of a fully automated COMET assay: 1.75 million single cells measured over a 5 year period. DNA Repair (Amst). 2011;10(3):322–37.
  20. 20. Guertler A, Kraemer A, Roessler U, Hornhardt S, Kulka U, Moertl S, et al. The WST survival assay: an easy and reliable method to screen radiation-sensitive individuals. Radiat Prot Dosimetry. 2011;143(2–4):487–90. pmid:21183542
  21. 21. Lopez-Romero P. Agi4x44PreProcess: PreProcessing of Agilent 4x44 array data. R package version 1.16.0.
  22. 22. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society. 1995;57:289–300.
  23. 23. Alsner J, Andreassen CN, Overgaard J. Genetic markers for prediction of normal tissue toxicity after radiotherapy. Seminars in Radiation Oncology. 2008;18:126–35. pmid:18314067
  24. 24. Andreassen CN. Can risk of radiotherapy-induced normal tissue complications be predicted from genetic profiles? Acta Oncologica. 2005;44:801–15. pmid:16332587
  25. 25. Liu H, Andrews DW, Evans JJ, Werner-Wasik M, Yu Y, Dicker AP, et al. Plan quality and treatment efficiency for radiosurgery to multiple brain metastases: non-coplanar RapidArc vs. Gamma Knife. Frontiers in Oncology. 2016;6(26).
  26. 26. Pauli A, Valen E, Lin MF, Garber M, Vastenhouw NL, Levin JZ, et al. Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Research. 2012;22:577–91. pmid:22110045
  27. 27. Lienau J, Schmidt-Bleek K, Peters A, Weber H, Bail HJ, Duda HN, et al. Insight into the molecular pathophysiology of delayed bone healing in a sheep model. Tissue Engineering Part A. 2010;16(1):191–9. pmid:19678759
  28. 28. Schell H, Thompson MS, Bail HJ, Hoffmanna J-E, Schilla A, Dudaa GN, et al. Mechanical induction of critically delayed bone healing in sheep: radiological and biomechanical results. Journal of Biomechanics. 2008;41(14):3066–72. pmid:18778822
  29. 29. Gürtler A, Hauptmann M, Pautz S, Kulka U, Friedl AA, Lehr S, et al. The inter-individual variability outperforms the intra-individual variability of differentially expressed proteins prior and post irradiation in lymphoblastoid cell lines. Arch Physiol Biochem. 2014;120(5):198–207. pmid:25174346
  30. 30. Azzama EI, Jay-Gerinb J-P, Pain D. Ionizing radiation-induced metabolic oxidative stress and prolonged cell injury. Cancer Letters. 2012;327(1–2):48–60. pmid:22182453
  31. 31. Li L, Story K, Legerski RJ. Cellular responses to ionizing radiation damage. International Journal of Radiation Oncology, Biology, Physics. 2001;49(4):1157–62. pmid:11240259
  32. 32. Jung M, Dritschilo A. Signal transduction and cellular responses to ionizing radiation. Seminars in Radiation Oncology. 1996;6(4):268–72. pmid:10717184
  33. 33. Koschützki D. Network Centralities. In: Junker BH, Schreiber F, editors. Analysis of Biological Networks: Wiley; 2007. p. 65–84.
  34. 34. Abbasi A, Hossain L. Hybrid Centrality Measures for Binary and Weighted Networks. In: Menezes R, Evsukoff A, González MC, editors. Complex Networks. 424. Berlin: Springer; 2013. p. 1–7.
  35. 35. Matthews L, Gopinath G, Gillespie M, Caudy M, Croft D, de Bono B, et al. Reactome knowledgebase of human biological pathways and processes. Nucleic Acids Res. 2009;37(Database issue):D619–D22. pmid:18981052
  36. 36. Sabin RJ, Anderson RM. Cellular Senescence—its role in cancer and the response to ionizing radiation. Genome Integr 2011; 2: 7 2011;2(1):7. pmid:21834983
  37. 37. Meng Y, Efimova EV, Hamzeh KW, Darga TE, Mauceri HJ, Fu YX, et al. Radiation-inducible immunotherapy for cancer: senescent tumor cells as a cancer vaccine. Molecular Therapy. 2012;20(5):1046–55. pmid:22334019
  38. 38. Freund A, Orjalo AV, Desprez PY, Campisi J. Inflammatory networks during cellular senescence: causes and consequences. Trends Mol Med. 2010;16(5):238–46. pmid:20444648
  39. 39. Nelson G, Wordsworth J, Wang C, Jurk D, Lawless C, Martin-Ruiz C, et al. A senescent cell bystander effect: senescence-induced senescence. Aging Cell. 2012;11(2):345–9. pmid:22321662
  40. 40. Wu PC, Wang Q, Grobman L, Chu E, Wu DY. Accelerated cellular senescence in solid tumor therapy. Experimental Oncology. 2012;34(3):298–305. pmid:23070015
  41. 41. Tsai KK, Stuart J, Chuang YY, Little JB, Yuan ZM. Low-dose radiation-induced senescent stromal fibroblasts render nearby breast cancer cells radioresistant. Radiation Research. 2009;172(3):306–13. pmid:19708779
  42. 42. Amundson SA, Grace MB, McLeland CB, Epperly MW, Yeager A, Zhan Q, et al. Human in vivo radiation-induced biomarkers: gene expression changes in radiotherapy patients. Cancer Research. 2004;64(18):6368–71. pmid:15374940
  43. 43. Badie C, Dziwura S, Raffy C, Tsigani T, Alsbeih G, Moody J, et al. Aberrant CDKN1A transcriptional response associates with abnormal sensitivity to radiation treatment. British Journal of Cancer. 2008;98(11):1845–51. pmid:18493234
  44. 44. Shah PP, Donahue G, Otte GL, Capell BC, Nelson DM, Cao K, et al. Lamin B1 depletion in senescent cells triggers large-scale changes in gene expression and the chromatin landscape. Genes & Development. 2013;27:1787–99.
  45. 45. Lee M, Urata SM, Aguilera JA, Perry CC, Milligan JR. Modeling the influence of histone proteins on the sensitivity of DNA to ionizing radiation. Radiation Research. 2012;177(2):152–63. pmid:22103271
  46. 46. di Masi A, D'Apice MR, Ricordy R, Tanzarella C, Novelli G. The R527H mutation in LMNA gene causes an increased sensitivity to ionizing radiation. Cell Cycle. 2008;7(13):2030–7. pmid:18604166
  47. 47. Xu M, Myerson R, Hunt C, Kumar S, Moros E, Straube B, et al. Treatment of cells with Mre11 siRNA increases radiation sensitivity and reduces heat induced radiosensitization. International Journal of Radiation Oncology Biology Physics. 2003;57(2):144–5.
  48. 48. Xu M, Myerson RJ, Hunt C, Kumar S, Moros EG, Straube WL, et al. Transfection of human tumour cells with Mre11 siRNA and the increase in radiation sensitivity and the reduction in heat-induced radiosensitization. International Journal of Hyperthermia. 2004;20(2).
  49. 49. Söderlund K, Stål O, Skoog L, Rutqvist LE, Nordenskjöld B, Askmalm MS. Intact Mre11/Rad50/Nbs1 complex predicts good response to radiotherapy in early breast cancer. International Journal of Radiation Oncology Biology Physics. 2007;68(1):50–8.
  50. 50. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology. 2004;3(1):Article 3.
  51. 51. R Core Team. R: A language and environment for statistical computing 2013.
  52. 52. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5(10):R80. pmid:15461798
  53. 53. Li Q, Birkbak NJ, Gyorffy B, Szallasi Z, Eklund AC. Jetset: selecting the optimal microarray probe set to represent a gene. BMC Bioinformatics. 2011;12(474).
  54. 54. Opgen-Rhein R, Strimmer K. Using regularized dynamic correlation to infer gene dependency networks from time-series microarray data. The 4th International Workshop on Computational Systems Biology, WCSB. 2006.
  55. 55. Wasserman S, Faust K. Social Network Analysis: Methods and Applications: Cambridge: Cambridge University Press; 1994.
  56. 56. Koschützki D, Schreiber F. Centrality analysis methods for biological networks and their application to gene regulatory networks. Gene Regulation and Systems Biology. 2008;2:193–201. pmid:19787083