BACKGROUND
the complex nature of human disease has long been recognized and, with the exception of a limited number of examples which follow the rules of mendelian inheritance patterns, common disease results from the poorly understood interaction of genetic and environmental factors  <cit> .
at the same time, gene-gene interactions that do not result in linearity between genotype and phenotype , may involve several genes at time, dramatically increasing the complexity of the phenomenon.
epistasis can either be defined from a biological point of view as deviations from the simple inheritance patterns observed by mendel  <cit>  or, from a mathematical point of view, as deviations from additivity in a linear statistical model  <cit> .
the study of statistical epistasis by traditional parametric models is challenging and hindered by several limitations.
these include, the problem of the sparseness of data into the multidimensional space  <cit> , the loss of power when adjusting for multiple testing to decrease type i error  <cit> , the loss of power in presence of multicollinearity  <cit>  or genetic heterogeneity  <cit> .
to address these issues, several non-parameteric multi-locus methods, essentially based on machine-learning techniques, have been developed and/or applied to genetic association studies with positive results  <cit> .
the application of data mining algorithms to detect non-linear high-order interactions in the context of survival analysis is more complex and thus far limited to a few examples  <cit> .
however, the effective ability of these algorithms to model gene-gene interactions and their power to detect epistasis in survival analysis has yet to be determined.
at least two points in modelling non-linear interactions in survival analysis should be taken into account.
the first, is the proper way to handle censored data, that is those cases for whom the outcome has not yet happened at the end of the observation time  or who did not have the event until the end of study , which are commonly referred to as right-censored cases  <cit> .
the second, is the optimal performance measure to be used in assessing a learned model in survival analysis.
in this paper we present an extension of the multifactor dimensionality reduction  algorithm  <cit> , to detect and characterize epistatic interactions in the context of survival analysis which was specifically designed to address the abovementioned issues.
censored data were directly handled by estimating individual multilocus cells survival functions by the kaplan-meier method  <cit> .
multilocus genotypes were then pooled into high-risk and low-risk groups whose predictive accuracy was evaluated by the brier score for censored samples proposed by graf et al  <cit> .
the power of the method we propose was at first evaluated in lifetime simulated datasets with epistatic effects which belonged to the most common survival distributions and with different degrees of right-censorship.
the method was then applied to identification of single-nucleotide polymorphisms  associated with responses to anti-tumor necrosis factor  agents in patients with rheumatoid arthritis  and active disease.
the notion of pharmacogenetics is not anew in ra and several candidate-gene studies have demonstrated a genetically-based individual variability to treatment with methotrexate or anti-tnf therapy  <cit> .
however, there is no consensus at present as to whether pharmacogenomics will allow prediction of anti-tnf therapy efficacy in ra.
so far pharmacogenomics studies in ra have produced conflicting results and population stratification and linkage disequilibrium have been cited as potential causes for the inability to replicate results of genetic association studies  <cit> .
yet, as demonstrated by greene et al  <cit>  when main effects fail to replicate, gene-gene interaction analysis should also be considered as a potential source of variance.
methods
description of the survival dimensionality reduction  algorithm
the core of the sdr algorithm is the classification procedure used to label as "high-risk" or "low-risk" the multilocus cells that result from gene-gene interaction.
this procedure will be used both for feature selection and for model validation as described in the forthcoming sections.
sdr assignments and evaluation
the sdr procedure for classification is illustrated in figure  <dig> and it involves  <dig> steps.
step 1: we firstly calculate by the kaplan-meier method  <cit>  the survival estimates Ŝ for the whole population in a dataset or dataset partition:
where, ni is the number of cases "at risk" of an event prior to ti, and δi, is the number of events at time ti.
step 2: we then select n discrete variables from the dataset and represent all the possible multidimensional cells resulting from their interaction.
for each multidimensional cell, survival estimates at each time interval Ŝc are calculated as described above.
step 3: the difference dc between the multilocus cell survival estimates and the whole population survival estimates, is calculated for each time interval ti:
step 4: all the dc for each multilocus cells are then averaged.
as the product-limit estimator is de facto a geometric progression, dc are averaged via the geometric mean  rather than via the algebraic mean.
as it is impossible to calculate the geometric with zero or negative data point values, these should be transformed to a meaningful equivalent positive number; being - <dig> <dc <  <dig>  transformation is made adding  <dig> to any dc value.
considering a finite number n of time intervals and denoting tn as the survival time at the nth time-interval we thus have:
step 5: cells with gmc ≤  <dig> are classified as "high-risk" and cells with gmc >  <dig> are classified as "low-risk".
examples from high-risk cells are pooled into one group and those from low-risk cells into another.
once dimensionality has been reduced to one dimension, sdr predictions may be evaluated via different fitness measure.
herein, we employed the brier score for censored data.
the brier score is a metric widely used for predicting the inaccuracy of a model and in its modified version proposed by graf et al  <cit> , it can incorporate censored samples.
the brier score bs for censored samples for a given t >  <dig>  is defined as:
where Ĝ denotes the kaplan-meier estimate of the censoring distribution g which is based on the observations  and i stands for the indicator function.
bs depends on time t, hence it makes sense to use the integrated brier score  as an overall measure for the prediction of the model at all times:
the lower the ibs the less inaccurate or, conversely, the more precise the prediction is.
feature selection and model validation
from a lifetime dataset, relevant features are extracted and validated via the k-fold cross-validation method, as described by ritchie and co-workers  <cit> .
for feature selection, the dataset is equally partitioned in k mutually-exclusive testing sets; one k set is retained for model validation whilst the remaining k- <dig> parts of the dataset are used as a training set.
this process is then repeated k-times with each of the k testing samples never included in the feature selection process.
in every training set, all the possible combinations of n-variables and the multilocus cells that result from their interactions are represented into the multidimensional space.
sdr models are then iteratively built for each combination of n-variables and training ibs scores are calculated.
for each n-combination of variables, the k training ibs are then averaged and the n-combination yielding the lowest mean ibs is selected and considered for model validation.
in the model validation phase, sdr assignments for the best n-combination of variables are determined in every training set.
on the basis of training assignments, the instances from the corresponding k testing sets are labelled as "high-risk" or "low-risk".
from these labels, a cross-validated ibs is then calculated; for this purpose, instead of mathematically averaging the k testing ibs, we merged the testing instances in a meta-analysis-based on individual patient data-fashion  <cit> .
let t <dig>  t <dig>  --- tk be the k testing sets, the individual patients data together by their assigned labels are merged sort to produce a larger tm testing set.
the ibs  is computed in the tm set and the n-combination yielding the lowest meta-ibs is then chosen as the final model.
this way both feature selection and model validation are used to determine the best epistatic model in lifetime datasets.
working on the merged dataset , we: 1) still ensure the independence of the testing sets, as testing assignments are calculated during the feature selection phase; 2) reduce the bias in the calculation of the bs, as the Ĝ and g weights used in the bs formula would otherwise be unreliably estimated in testing sets of limited size; 3) avoid to solely rely on a measure of central tendency  to estimate the predictive accuracy of our model, utterly ignoring any measure of variance.
data simulation and power calculation
simulated epistatic datasets were modelled upon five different survival distributions, described by the logistic-exponential equation  <cit> ; exponential , bathtub-shaped failure rate , upside-down bathtub-shaped failure rate , decreasing failure rate  and increasing failure rate :
where, s is the logistic-exponential survival distribution, t is the survival time, λ is a positive scale parameter and κ is a positive shape parameter and θ is a ≥  <dig> parameter that shifts the distribution to the left.
λ, κ and θ were adjusted so that the cumulative prevalence of the event at the end of the observation time tn or k, was equal to an arbitrary value of  <dig>  .
for data simulation tn was set to  <dig> time units.
according to culverhouse et al  <cit> , we then generated different epistatic models for two biallelic snps a and b, both in hardy-weinberg equilibrium  and with minor allele frequency  =  <dig> , so that k = ka = kb, where ka and kb are the marginal penetrances for snp a and snp b. as these models were adjusted to fit the cumulative prevalence at tn, their broad-sense heritability  was considered to be a cumulative estimate of h <dig> at tn or h <dig>  for data simulation h <dig> was set to  <dig> ,  <dig> ,  <dig>  or  <dig> .
time-point multilocus genotype penetrances- henceforth labelled as ti -, were adjusted so that the cumulative prevalence at each time-point k would fit the survival distribution.
we assumed that the relation between snp a and snp b is purely epistatic at each time-interval ti, that is k = ka = kb with  <dig> <ti ≤ tn, and that the two-locus model is always proportional at the different ti.
hence, applying the product-limit estimation, we obtain:
where kti is the time-point prevalence at ti and; ti ≤ tn.
from kti we can calculate ti from the cumulative multilocus genotype penetrances  previously used to compute k and h2:
similarly, time-point cumulative estimated multilocus genotype penetrances  are proportional to k. these penetrances can be used to derive the time-point cumulative estimated h <dig> or h <dig>
once ti values had been calculated, a population of  <dig> individuals was built considering all the time-intervals  <dig> < ti ≤ tn.
herein, from  <dig> survival distributions with  <dig> different h <dig> and two tn we obtained  <dig> populations where the outcome was related to the epistatic interaction between snp a and snp b .
to each of these populations  <dig> unrelated snps in hwe, with maf ranging from  <dig>  to  <dig>  were added.
an additional 5% censoring/year was also added to account for hypothetical non-event related causes of withdrawal from observation.
from the simulated populations we finally randomly draw  <dig> samples of  <dig> cases and  <dig> controls  or  <dig> cases and  <dig> controls .
a total of  <dig>  datasets were then generated for simulation.
power was estimated as the number of times sdr correctly identified the two functional snps out of  <dig> datasets/model/degree-of-censorship.
datasets can be obtained upon request from the authors.
application of the sdr algorithm to the ra dataset
the sdr algorithm was then tested in a real-world dataset which consists of previously unpublished data about  <dig> dutch patients with  a diagnosis of ra according to acr criteria  <cit> ,  a disease activity score  > <dig>   <cit>  and  previous treatment with at least two other anti-rheumatics including methotrexate  at an optimal dose  or intolerance for mtx, that underwent treatment with anti-tnfα agents.
these patients were extracted from the dutch rheumatoid arthritis monitoring  registry  <cit>  and genotyped for  <dig> candidate snps and evaluated every  <dig> months to ascertain whether they had reached a clinical remission, defined as a das <dig> ≤  <dig>   <cit> .
genotyping details can be found in: pavy et al  <cit> , coenen et al  <cit> , toonen et al  <cit>  and alizadeh et al  <cit> .
the choice of the studied snps was motivated by results from previous association and/or pharmacogenomics studies in ra  <cit> .
a detailed list of the analysed genetic variants is provided in additional file  <dig>  table s <dig>  the k-nearest-neighbour method was used to impute genotypes with missing data <10%  <cit> .
the open-source orange data mining software  was used for imputation.
overall the right-censorship of this dataset was 68%; 5-fold cross-validation was used for the sdr analysis.
an empirical p-value for the sdr results was calculated by performing 100-fold permutation testing  <cit> .
a whole sdr analysis, up to the 3rd dimension was conducted in the permutated datasets during the permutation procedure.
for all the analyses a modified version of the freely available sdr algorithm written in python http://sourceforge.net/projects/sdrproject/ was used.
RESULTS
simulated datasets
the penetrance functions for the simulated datasets are reported in additional file 1; as it can be observed, time-point h <dig> across the different models were consistently low, with a median of  <dig>  .
power for the sdr algorithm to correctly identify the causative pair of snps in the  <dig> simulated survival models with  <dig> different degrees of censorship is reported in table  <dig>  overall, the median power across all the datasets was  <dig> % .
the relationship between h <dig> and power followed a direct logistic distribution as shown in figure  <dig>  .
moreover, the power resulted to be independent from the survival distribution the sample was withdrawn from whilst it was inversely related to the degree of censorship.
the h <dig> values we employed correspond to mild-to-moderate hazard ratios  for high-risk vs low-risk combinations in the different epistatic models: hr =  <dig>  for h <dig> =  <dig> , hr =  <dig>  for h <dig> =  <dig> , hr =  <dig>  for h <dig> =  <dig>  and hr =  <dig>  for h <dig> =  <dig> .
altogether these results suggest that sdr has a satisfactory power to identify a pair of causative genes in purely epistatic lifetime models with mild-to-moderate effect size and for which the proportionality of hazards holds.
power of the sdr algorithm in simulated datasets modelled using the  <dig> classes of survival distribution described by the logistic-exponential equation.
different degrees of right-censoring rate , and cumulative heritability at the survival time  are considered.
ubt, upside-down bathtub-shaped failure rate; dfr, decreasing failure rate; ifr, increasing failure rate; bt bathtub-shaped failure rate; exp, exponential.
ra dataset
the ra datasets comprises  <dig> anti-tnfα-treated patients, aging  <dig> ±  <dig>   at the onset of disease;  <dig> % of patients tested positive for the rheumatoid factor and  <dig>  were females.
one-hundred-forty-five patients  were treated with adalimumab,  <dig>  with infliximab and  <dig>  with etanercept; overall  <dig> patients  were treated with anti-tnf antibodies .
das <dig> at the beginning of therapy was  <dig>  ±  <dig> .
clinical remission, based on das <dig>  was observed in  <dig> cases .
details about the genotyped snps along with their frequencies in the studied population are reported in additional file  <dig>  table s1; all the snps were in hwe.
none of the single snps showed a statistically significant association with response to anti-tnf agents, either under a dominant or recessive model, as illustrated in additional file  <dig>  table s <dig> .
the sdr algorithm sorted out two-way interaction model, involving the rs <dig>  and the rs <dig>  snps, as the most predictive for responses to anti-tnf therapy in patients with active ra.
table  <dig> shows the full analysis conducted by the sdr algorithm in the ra datasets.
as expected, the model overfits in the training population as the number of snps included in the model increases.
yet, cross-validation prevented over-fitting as the minimum meta-ibs was observed for the 2-way interaction, that was thus chosen as the best epistatic model.
this model was significant at the  <dig>  threshold after 100-fold permutation testing.
figure 3a summarizes the multilocus cells for the rs <dig> × rs <dig> interaction along with sdr "high-" and "low-risk" assignments .
this interaction had the typical non-linear behaviour of epistatic model.
plotting this sdr assignments we can observe that patients labelled as "responders" achieved earlier and higher rates of clinical remission after anti-tnf therapy compared to patients labelled as "non-responders"   <cit> .
selection of the best combination of attributes by the sdr method.
the model with the minimum testing ibs value in the cross-validated testing sets is indicated in boldface type.
p values associated with the log-rank test statistic calculated by the 100-fold permutation test.
snp, single nucleotide polymorphism.
discussion
in the present paper we introduce sdr, an algorithm specifically conceived to detect non-linear gene-gene interactions in presence of right-censored data.
the need for such a bioinformatics tool comes from the observation that several studies in the medical field deal with the loss of data during the period of study, and that methods that do not take into account censored data give upwardly biased estimates of failure and/or suffer from information loss due to reduced sample size.
cox regression  <cit> , the most popular statistical technique used to analyse time-to-event multivariate model, is not adequate to detect non-linearities either.
indeed, to properly model interaction in cox regression, the user should have a priori knowledge of the variable relationships and may need to enter nonlinear transforms of the predictors, but this is often a trial and error approach.
also, the number of polynomial terms needed to model complex interactions is inflated as the number of predictors grows, increasing standard errors and, thus, type i error  <cit> .
using simulated epistatic lifetime datasets, we demonstrated that the sdr algorithm retains a fully satisfactory power to sort out a set of causative genes with mild-to-moderate epistatic effect size from a pool of candidate genes.
these results have been accomplished by the intrinsic properties of the sdr methodology.
firstly, sdr is non-parametric, in the sense that it is not necessary to make a priori assumptions about the underlying interaction model.
also, sdr requires no assumptions concerning the nature or shape of the underlying survival distribution.
secondly, sdr performs well even in datasets where the right-censorship rate is high .
notably, when we run mdr ignoring censorship on the same simulated datasets we observed a dramatic reduction in the power to detect the causative pair of genes, that ranged from 5% to more than 90% .
beside power, an additional advantage of sdr is that, combining cross-validation with permutation testing, the chance of false-positive findings is minimized  <cit> .
yet, to generate the simulated datasets we required the hazards among multilocus cells to be proportional along the lifetime distribution.
hence, further simulations are needed to ascertain whether sdr is suitable for situations where this assumption does not hold, as for instance in case of additive hazard models.
similarly, further simulations are needed to establish sdr performance in larger datasets, in presence of genetic heterogeneity, linkage disequilibrium, or different ranges of maf.
sdr is, de facto, an extension of the mdr algorithm optimized to analyse lifetime distributions, hence it suffers from similar limitations  <cit>  and it shares some peculiarities with the latter.
namely, the power of sdr is influenced by the epistatic effect size, which is strictly related to the  heritability of the model  <cit> .
moreover, as with mdr, the biological significance of sdr models may be difficult to interpret due to the non-linear distribution of high-risk and low-risk cells across the multidimensional space  <cit> .
finally, sdr cannot make predictions when multilocus cells contain no data and gm estimates may be upwardly or downwardly inflated when multilocus cells contain few data.
having demonstrated the capability of sdr to detect gene-gene interactions in lifetime datasets, we applied the algorithm to a population of patients with active ra to identify epistatic interactions that may affect time-related responses to anti-tnf biological agents.
we did show that among a set of  <dig> candidate-gene loci, none of which had a detectable marginal effect on the outcome variable, the non-linear interaction between the rs <dig>  and the rs <dig>  snps significantly predicted the responses to anti-tnf therapy.
whilst it is difficult to dissect the biological meaning of statistical epistatic models  <cit> , it should be noted that several lines of evidence support a role for fcγ receptors and irf <dig> in rheumatoid arthritis and tnf driven processes.
associations between irf <dig> polymorphisms and ra have been described in different populations  <cit>  and a recent genome-wide association study  showed that the rs <dig> variant in the irf <dig> gene, in tight linkage disequilibrium with rs <dig>  is strongly associated with ra susceptibility  <cit> .
of interest, functional experiments demonstrated that the rs <dig> snp significantly alters irf <dig> mrna expression  <cit> .
the association between irf <dig> gene and ra may be linked, at least in part, to the ability of irf <dig> to regulate the secretion of pro-inflammatory cytokines.
indeed, takaoka et al  <cit>  using mouse models deficient in the irf <dig> gene, showed that irf <dig> is generally involved downstream of the toll-like receptor -myd <dig> signalling pathway for gene induction of tnfα and other cytokines relevant to the pathogenesis of ra.
of interest, the use of anti-tnf agents was shown to decrease tlrs expression on different cellular types  <cit> .
similarly to irf <dig>  the fcγriia is involved in tnfα production in the rheumatoid synovia, as observed by clavel and co-workers  <cit> .
the interaction between the genetic variants of irf <dig> and fcγ receptors could thus influence tnfα production and/or availability, affecting the clinical response to anti-tnf agents.
additionally, as postulated by cañete at al  <cit> , polymorphisms of the fcγriia may alter the clearance rate of anti-tnf antibodies modulating plasma concentrations and consequently their biological effect in subjects with active ra.
theoretically, this effect should not restricted only to anti-tnf antibodies, as also anti-tnfα receptors  contain a fc portion of igg <dig> capable of binding to fcγ receptors to produce biological effects, such as antibody-dependent cell-mediated cytotoxicity  <cit> .
CONCLUSIONS
herein we introduced sdr, an innovative algorithm to detect epistasis in lifetime datasets.
simulation studies and application in a real-world setting, demonstrate the capability of sdr to detect non linear gene-gene interactions in studies aimed at evaluating the effect of candidate genes on time-dependent outcomes.
further studies are necessary to evaluate its applicability in large-scale datasets as well.
supplementary material
additional file 1
epistatic models and simulation specifics.
the file contains the settings used to generate the five survival distributions upon which epistatic models were modelled.
for each epistatic model, ttime-point and cumulative multilocus genotype penetrances are reported, along with the time-point and the cumulative broad-sense heritability and prevalence of the event.
click here for file
additional file 2
genotype frequencies and univariate analysis in the rheumatoid arthritis  dataset.
the file lists in tabular form the single nucleotide polymorphisms  included in the ra case-control study.
it also reports the associations by the log-rank test statistics under the dominant and recessive model between the studied snps and the occurrence of clinical remission after therapy with anti-tnf agents.
click here for file

 the complex nature of human disease has long been recognized and, with the exception of a limited number of examples which follow the rules of mendelian inheritance patterns, common disease results from the poorly understood interaction of genetic and environmental factors. at the same time, gene-gene interactions that do not result in linearity between genotype and phenotype, may involve several genes at time, dramatically increasing the complexity of the phenomenon. epistasis can either be defined from a biological point of view as deviations from the simple inheritance patterns observed by mendel or, from a mathematical point of view, as deviations from additivity in a linear statistical model. the study of statistical epistasis by traditional parametric models is challenging and hindered by several limitations. these include, the problem of the sparseness of data into the multidimensional space, the loss of power when adjusting for multiple testing to decrease type i error, the loss of power in presence of multicollinearity or genetic heterogeneity. to address these issues, several non-parameteric multi-locus methods, essentially based on machine-learning techniques, have been developed and/or applied to genetic association studies with positive results. the application of data mining algorithms to detect non-linear high-order interactions in the context of survival analysis is more complex and thus far limited to a few examples. however, the effective ability of these algorithms to model gene-gene interactions and their power to detect epistasis in survival analysis has yet to be determined. at least two points in modelling non-linear interactions in survival analysis should be taken into account. the first, is the proper way to handle censored data, that is those cases for whom the outcome has not yet happened at the end of the observation time or who did not have the event until the end of study, which are commonly referred to as right-censored cases. the second, is the optimal performance measure to be used in assessing a learned model in survival analysis. in this paper we present an extension of the multifactor dimensionality reduction algorithm, to detect and characterize epistatic interactions in the context of survival analysis which was specifically designed to address the abovementioned issues. censored data were directly handled by estimating individual multilocus cells survival functions by the kaplan-meier method. multilocus genotypes were then pooled into high-risk and low-risk groups whose predictive accuracy was evaluated by the brier score for censored samples proposed by graf et al. the power of the method we propose was at first evaluated in lifetime simulated datasets with epistatic effects which belonged to the most common survival distributions and with different degrees of right-censorship. the method was then applied to identification of single-nucleotide polymorphisms associated with responses to anti-tumor necrosis factor agents in patients with rheumatoid arthritis and active disease. the notion of pharmacogenetics is not anew in ra and several candidate-gene studies have demonstrated a genetically-based individual variability to treatment with methotrexate or anti-tnf therapy. however, there is no consensus at present as to whether pharmacogenomics will allow prediction of anti-tnf therapy efficacy in ra. so far pharmacogenomics studies in ra have produced conflicting results and population stratification and linkage disequilibrium have been cited as potential causes for the inability to replicate results of genetic association studies. yet, as demonstrated by greene et al when main effects fail to replicate, gene-gene interaction analysis should also be considered as a potential source of variance. methods description of the survival dimensionality reduction algorithm the core of the sdr algorithm is the classification procedure used to label as "high-risk" or "low-risk" the multilocus cells that result from gene-gene interaction. this procedure will be used both for feature selection and for model validation as described in the forthcoming sections. sdr assignments and evaluation the sdr procedure for classification is illustrated in figure and it involves steps. step 1: we firstly calculate by the kaplan-meier method the survival estimates Ŝ for the whole population in a dataset or dataset partition: where, ni is the number of cases "at risk" of an event prior to ti, and δi, is the number of events at time ti. step 2: we then select n discrete variables from the dataset and represent all the possible multidimensional cells resulting from their interaction. for each multidimensional cell, survival estimates at each time interval Ŝc are calculated as described above. step 3: the difference dc between the multilocus cell survival estimates and the whole population survival estimates, is calculated for each time interval ti: step 4: all the dc for each multilocus cells are then averaged. as the product-limit estimator is de facto a geometric progression, dc are averaged via the geometric mean rather than via the algebraic mean. as it is impossible to calculate the geometric with zero or negative data point values, these should be transformed to a meaningful equivalent positive number; being - <dc < transformation is made adding to any dc value. considering a finite number n of time intervals and denoting tn as the survival time at the nth time-interval we thus have: step 5: cells with gmc ≤ are classified as "high-risk" and cells with gmc > are classified as "low-risk". examples from high-risk cells are pooled into one group and those from low-risk cells into another. once dimensionality has been reduced to one dimension, sdr predictions may be evaluated via different fitness measure. herein, we employed the brier score for censored data. the brier score is a metric widely used for predicting the inaccuracy of a model and in its modified version proposed by graf et al, it can incorporate censored samples. the brier score bs for censored samples for a given t > is defined as: where Ĝ denotes the kaplan-meier estimate of the censoring distribution g which is based on the observations and i stands for the indicator function. bs depends on time t, hence it makes sense to use the integrated brier score as an overall measure for the prediction of the model at all times: the lower the ibs the less inaccurate or, conversely, the more precise the prediction is. feature selection and model validation from a lifetime dataset, relevant features are extracted and validated via the k-fold cross-validation method, as described by ritchie and co-workers. for feature selection, the dataset is equally partitioned in k mutually-exclusive testing sets; one k set is retained for model validation whilst the remaining k- parts of the dataset are used as a training set. this process is then repeated k-times with each of the k testing samples never included in the feature selection process. in every training set, all the possible combinations of n-variables and the multilocus cells that result from their interactions are represented into the multidimensional space. sdr models are then iteratively built for each combination of n-variables and training ibs scores are calculated. for each n-combination of variables, the k training ibs are then averaged and the n-combination yielding the lowest mean ibs is selected and considered for model validation. in the model validation phase, sdr assignments for the best n-combination of variables are determined in every training set. on the basis of training assignments, the instances from the corresponding k testing sets are labelled as "high-risk" or "low-risk". from these labels, a cross-validated ibs is then calculated; for this purpose, instead of mathematically averaging the k testing ibs, we merged the testing instances in a meta-analysis-based on individual patient data-fashion. let t t --- tk be the k testing sets, the individual patients data together by their assigned labels are merged sort to produce a larger tm testing set. the ibs is computed in the tm set and the n-combination yielding the lowest meta-ibs is then chosen as the final model. this way both feature selection and model validation are used to determine the best epistatic model in lifetime datasets. working on the merged dataset, we: 1) still ensure the independence of the testing sets, as testing assignments are calculated during the feature selection phase; 2) reduce the bias in the calculation of the bs, as the Ĝ and g weights used in the bs formula would otherwise be unreliably estimated in testing sets of limited size; 3) avoid to solely rely on a measure of central tendency to estimate the predictive accuracy of our model, utterly ignoring any measure of variance. data simulation and power calculation simulated epistatic datasets were modelled upon five different survival distributions, described by the logistic-exponential equation; exponential, bathtub-shaped failure rate, upside-down bathtub-shaped failure rate, decreasing failure rate and increasing failure rate: where, s is the logistic-exponential survival distribution, t is the survival time, λ is a positive scale parameter and κ is a positive shape parameter and θ is a ≥ parameter that shifts the distribution to the left. λ, κ and θ were adjusted so that the cumulative prevalence of the event at the end of the observation time tn or k, was equal to an arbitrary value of. for data simulation tn was set to time units. according to culverhouse et al, we then generated different epistatic models for two biallelic snps a and b, both in hardy-weinberg equilibrium and with minor allele frequency =, so that k = ka = kb, where ka and kb are the marginal penetrances for snp a and snp b. as these models were adjusted to fit the cumulative prevalence at tn, their broad-sense heritability was considered to be a cumulative estimate of h at tn or h for data simulation h was set to,, or. time-point multilocus genotype penetrances- henceforth labelled as ti -, were adjusted so that the cumulative prevalence at each time-point k would fit the survival distribution. we assumed that the relation between snp a and snp b is purely epistatic at each time-interval ti, that is k = ka = kb with <ti ≤ tn, and that the two-locus model is always proportional at the different ti. hence, applying the product-limit estimation, we obtain: where kti is the time-point prevalence at ti and; ti ≤ tn. from kti we can calculate ti from the cumulative multilocus genotype penetrances previously used to compute k and h2: similarly, time-point cumulative estimated multilocus genotype penetrances are proportional to k. these penetrances can be used to derive the time-point cumulative estimated h or h once ti values had been calculated, a population of individuals was built considering all the time-intervals < ti ≤ tn. herein, from survival distributions with different h and two tn we obtained populations where the outcome was related to the epistatic interaction between snp a and snp b. to each of these populations unrelated snps in hwe, with maf ranging from to were added. an additional 5% censoring/year was also added to account for hypothetical non-event related causes of withdrawal from observation. from the simulated populations we finally randomly draw samples of cases and controls or cases and controls. a total of datasets were then generated for simulation. power was estimated as the number of times sdr correctly identified the two functional snps out of datasets/model/degree-of-censorship. datasets can be obtained upon request from the authors. application of the sdr algorithm to the ra dataset the sdr algorithm was then tested in a real-world dataset which consists of previously unpublished data about dutch patients with a diagnosis of ra according to acr criteria, a disease activity score > and previous treatment with at least two other anti-rheumatics including methotrexate at an optimal dose or intolerance for mtx, that underwent treatment with anti-tnfα agents. these patients were extracted from the dutch rheumatoid arthritis monitoring registry and genotyped for candidate snps and evaluated every months to ascertain whether they had reached a clinical remission, defined as a das ≤. genotyping details can be found in: pavy et al, coenen et al, toonen et al and alizadeh et al. the choice of the studied snps was motivated by results from previous association and/or pharmacogenomics studies in ra. a detailed list of the analysed genetic variants is provided in additional file table s the k-nearest-neighbour method was used to impute genotypes with missing data <10%. the open-source orange data mining software was used for imputation. overall the right-censorship of this dataset was 68%; 5-fold cross-validation was used for the sdr analysis. an empirical p-value for the sdr results was calculated by performing 100-fold permutation testing. a whole sdr analysis, up to the 3rd dimension was conducted in the permutated datasets during the permutation procedure. for all the analyses a modified version of the freely available sdr algorithm written in python http://sourceforge.net/projects/sdrproject/ was used.  simulated datasets the penetrance functions for the simulated datasets are reported in additional file 1; as it can be observed, time-point h across the different models were consistently low, with a median of. power for the sdr algorithm to correctly identify the causative pair of snps in the simulated survival models with different degrees of censorship is reported in table overall, the median power across all the datasets was %. the relationship between h and power followed a direct logistic distribution as shown in figure. moreover, the power resulted to be independent from the survival distribution the sample was withdrawn from whilst it was inversely related to the degree of censorship. the h values we employed correspond to mild-to-moderate hazard ratios for high-risk vs low-risk combinations in the different epistatic models: hr = for h =, hr = for h =, hr = for h = and hr = for h =. altogether these results suggest that sdr has a satisfactory power to identify a pair of causative genes in purely epistatic lifetime models with mild-to-moderate effect size and for which the proportionality of hazards holds. power of the sdr algorithm in simulated datasets modelled using the classes of survival distribution described by the logistic-exponential equation. different degrees of right-censoring rate, and cumulative heritability at the survival time are considered. ubt, upside-down bathtub-shaped failure rate; dfr, decreasing failure rate; ifr, increasing failure rate; bt bathtub-shaped failure rate; exp, exponential. ra dataset the ra datasets comprises anti-tnfα-treated patients, aging ± at the onset of disease; % of patients tested positive for the rheumatoid factor and were females. one-hundred-forty-five patients were treated with adalimumab, with infliximab and with etanercept; overall patients were treated with anti-tnf antibodies. das at the beginning of therapy was ±. clinical remission, based on das was observed in cases. details about the genotyped snps along with their frequencies in the studied population are reported in additional file table s1; all the snps were in hwe. none of the single snps showed a statistically significant association with response to anti-tnf agents, either under a dominant or recessive model, as illustrated in additional file table s. the sdr algorithm sorted out two-way interaction model, involving the rs and the rs snps, as the most predictive for responses to anti-tnf therapy in patients with active ra. table shows the full analysis conducted by the sdr algorithm in the ra datasets. as expected, the model overfits in the training population as the number of snps included in the model increases. yet, cross-validation prevented over-fitting as the minimum meta-ibs was observed for the 2-way interaction, that was thus chosen as the best epistatic model. this model was significant at the threshold after 100-fold permutation testing. figure 3a summarizes the multilocus cells for the rs × rs interaction along with sdr "high-" and "low-risk" assignments. this interaction had the typical non-linear behaviour of epistatic model. plotting this sdr assignments we can observe that patients labelled as "responders" achieved earlier and higher rates of clinical remission after anti-tnf therapy compared to patients labelled as "non-responders". selection of the best combination of attributes by the sdr method. the model with the minimum testing ibs value in the cross-validated testing sets is indicated in boldface type. p values associated with the log-rank test statistic calculated by the 100-fold permutation test. snp, single nucleotide polymorphism. discussion in the present paper we introduce sdr, an algorithm specifically conceived to detect non-linear gene-gene interactions in presence of right-censored data. the need for such a bioinformatics tool comes from the observation that several studies in the medical field deal with the loss of data during the period of study, and that methods that do not take into account censored data give upwardly biased estimates of failure and/or suffer from information loss due to reduced sample size. cox regression, the most popular statistical technique used to analyse time-to-event multivariate model, is not adequate to detect non-linearities either. indeed, to properly model interaction in cox regression, the user should have a priori knowledge of the variable relationships and may need to enter nonlinear transforms of the predictors, but this is often a trial and error approach. also, the number of polynomial terms needed to model complex interactions is inflated as the number of predictors grows, increasing standard errors and, thus, type i error. using simulated epistatic lifetime datasets, we demonstrated that the sdr algorithm retains a fully satisfactory power to sort out a set of causative genes with mild-to-moderate epistatic effect size from a pool of candidate genes. these results have been accomplished by the intrinsic properties of the sdr methodology. firstly, sdr is non-parametric, in the sense that it is not necessary to make a priori assumptions about the underlying interaction model. also, sdr requires no assumptions concerning the nature or shape of the underlying survival distribution. secondly, sdr performs well even in datasets where the right-censorship rate is high. notably, when we run mdr ignoring censorship on the same simulated datasets we observed a dramatic reduction in the power to detect the causative pair of genes, that ranged from 5% to more than 90%. beside power, an additional advantage of sdr is that, combining cross-validation with permutation testing, the chance of false-positive findings is minimized. yet, to generate the simulated datasets we required the hazards among multilocus cells to be proportional along the lifetime distribution. hence, further simulations are needed to ascertain whether sdr is suitable for situations where this assumption does not hold, as for instance in case of additive hazard models. similarly, further simulations are needed to establish sdr performance in larger datasets, in presence of genetic heterogeneity, linkage disequilibrium, or different ranges of maf. sdr is, de facto, an extension of the mdr algorithm optimized to analyse lifetime distributions, hence it suffers from similar limitations and it shares some peculiarities with the latter. namely, the power of sdr is influenced by the epistatic effect size, which is strictly related to the heritability of the model. moreover, as with mdr, the biological significance of sdr models may be difficult to interpret due to the non-linear distribution of high-risk and low-risk cells across the multidimensional space. finally, sdr cannot make predictions when multilocus cells contain no data and gm estimates may be upwardly or downwardly inflated when multilocus cells contain few data. having demonstrated the capability of sdr to detect gene-gene interactions in lifetime datasets, we applied the algorithm to a population of patients with active ra to identify epistatic interactions that may affect time-related responses to anti-tnf biological agents. we did show that among a set of candidate-gene loci, none of which had a detectable marginal effect on the outcome variable, the non-linear interaction between the rs and the rs snps significantly predicted the responses to anti-tnf therapy. whilst it is difficult to dissect the biological meaning of statistical epistatic models, it should be noted that several lines of evidence support a role for fcγ receptors and irf in rheumatoid arthritis and tnf driven processes. associations between irf polymorphisms and ra have been described in different populations and a recent genome-wide association study showed that the rs variant in the irf gene, in tight linkage disequilibrium with rs is strongly associated with ra susceptibility. of interest, functional experiments demonstrated that the rs snp significantly alters irf mrna expression. the association between irf gene and ra may be linked, at least in part, to the ability of irf to regulate the secretion of pro-inflammatory cytokines. indeed, takaoka et al using mouse models deficient in the irf gene, showed that irf is generally involved downstream of the toll-like receptor -myd signalling pathway for gene induction of tnfα and other cytokines relevant to the pathogenesis of ra. of interest, the use of anti-tnf agents was shown to decrease tlrs expression on different cellular types. similarly to irf the fcγriia is involved in tnfα production in the rheumatoid synovia, as observed by clavel and co-workers. the interaction between the genetic variants of irf and fcγ receptors could thus influence tnfα production and/or availability, affecting the clinical response to anti-tnf agents. additionally, as postulated by cañete at al, polymorphisms of the fcγriia may alter the clearance rate of anti-tnf antibodies modulating plasma concentrations and consequently their biological effect in subjects with active ra. theoretically, this effect should not restricted only to anti-tnf antibodies, as also anti-tnfα receptors contain a fc portion of igg capable of binding to fcγ receptors to produce biological effects, such as antibody-dependent cell-mediated cytotoxicity.  herein we introduced sdr, an innovative algorithm to detect epistasis in lifetime datasets. simulation studies and application in a real-world setting, demonstrate the capability of sdr to detect non linear gene-gene interactions in studies aimed at evaluating the effect of candidate genes on time-dependent outcomes. further studies are necessary to evaluate its applicability in large-scale datasets as well. supplementary material additional file 1 epistatic models and simulation specifics. the file contains the settings used to generate the five survival distributions upon which epistatic models were modelled. for each epistatic model, ttime-point and cumulative multilocus genotype penetrances are reported, along with the time-point and the cumulative broad-sense heritability and prevalence of the event. click here for file additional file 2 genotype frequencies and univariate analysis in the rheumatoid arthritis dataset. the file lists in tabular form the single nucleotide polymorphisms included in the ra case-control study. it also reports the associations by the log-rank test statistics under the dominant and recessive model between the studied snps and the occurrence of clinical remission after therapy with anti-tnf agents. click here for file
the sdr method was then applied to a series of dutch patients with active rheumatoid arthritis that were treated with anti-tnf biological agents. the algorithm requires neither specification about the underlying survival distribution nor about the underlying interaction model and proved satisfactorily powerful to detect a set of causative genes in synthetic epistatic lifetime datasets with a limited number of samples and high degree of right-censorship. simulation studies and application in a real-world setting support the capability of the sdr algorithm to model epistatic interactions in candidate-genes studies in presence of right-censored data. herein we introduce the survival dimensionality reduction algorithm, a non-parametric method specifically designed to detect epistasis in lifetime datasets. among a set of candidate genes, none of which showed a detectable marginal effect on anti-tnf responses, the sdr algorithm did find that the rs snp in the fcγriia gene and the rs snp in the irf gene non-linearly interact to predict clinical remission after anti-tnf biologicals. several computational approaches have been developed to model gene-gene interactions in case-control studies, however, none of them is suitable for time-dependent analysis.
 epistasis is recognized as a fundamental part of the genetic architecture of individuals. several computational approaches have been developed to model gene-gene interactions in case-control studies, however, none of them is suitable for time-dependent analysis. herein we introduce the survival dimensionality reduction algorithm, a non-parametric method specifically designed to detect epistasis in lifetime datasets.  the algorithm requires neither specification about the underlying survival distribution nor about the underlying interaction model and proved satisfactorily powerful to detect a set of causative genes in synthetic epistatic lifetime datasets with a limited number of samples and high degree of right-censorship. the sdr method was then applied to a series of dutch patients with active rheumatoid arthritis that were treated with anti-tnf biological agents. among a set of candidate genes, none of which showed a detectable marginal effect on anti-tnf responses, the sdr algorithm did find that the rs snp in the fcγriia gene and the rs snp in the irf gene non-linearly interact to predict clinical remission after anti-tnf biologicals.  simulation studies and application in a real-world setting support the capability of the sdr algorithm to model epistatic interactions in candidate-genes studies in presence of right-censored data. availability: http://sourceforge.net/projects/sdrproject/ 