The package implements a two stage method for finding significant interaction effects within a high-dimensional setting where the covariates outnumber the subjects assuming a Cox proportional hazards model.
Since the package is not yet available on CRAN, installation goes as follows:
invisible(devtools::install_github("luc-van-schijndel/twostageGWASsurvival"))
#> Using GitHub PAT from the git credential store.
#> Skipping install of 'twostageGWASsurvival' from a github remote, the SHA1 (c35d4e84) has not changed since last install.
#> Use `force = TRUE` to force installationThe two stage method is, as the name implies, based on two stages. The first stage tests all covariates individually for a significant marginal effect. The covariates that are found to be marginally significant are passed on to the next stage. The second stage then tests all pairs of marginal significant for a significant interaction effect. Since we are testing multiple hypotheses, we have to correct for this fact, which is of course done by applying a multiple hypotheses correction such as the Bonferroni correction for example. The crucial part is that the two stages, more specifically the p-values found in the two stages, are independent of each other. The proof of this is quite long, but it follows from a thorough analysis of the Fisher information matrix. By this independence, we only need to apply the multiple hypotheses correction to the hypotheses tested in the second stage.
By viewing the first stage as a filter of some sort, you get to see the power of this two stage method. Due to the filter, the number of hypotheses tested in the second stage decreases quickly as the filter becomes more strict. The number of pairs tested in the second stage grows quadratically, so if only one in one hundred covariates are found to be marginally significant, the second stage tests only one in every ten thousand pairs.
This allows for a fast computation to a relatively complex problem. Especially in genetic research, this is a big boon, since the number of genetic covariates can be as large as several million and consequently, the number of pairs of covariates can be as large as several trillion. By applying a strict filter, this number of pairs can be reduced by a large enough factor that the analysis of these interactions becomes possible.
This filter is specified by the value of the first stage threshold. This is the threshold the p-values of the first stage need to beat in order to be passed to the second stage. Under the null hypotheses, these p-values are uniformly distributed, so we expect that the fraction specified by the first stage threshold would be marginally significant and therefore passed to the second stage. This value can be tweaked by the user to achieve a suitable analysis.
By lowering the first stage threshold, less pairs of covariates will be tested in the second stage, resulting in a less strict multiple hypotheses correction and therefore increasing power. This is another benefit, aside from decreasing computation times. The filter must not be too strict however, since we do want any actual effects to be detected. Consider this value therefore properly.
A small example is given by the following datasets:
survival.dataset <- survival::Surv(c(5,5,3,3,2,2,2,1,1,1),
c(0,0,1,1,1,1,1,1,1,1))
covariate.matrix <- matrix(c(2,2,1,
2,2,1,
1,2,1,
2,1,1,
1,1,1,
1,1,1,
1,0,0,
0,1,0,
1,0,0,
0,0,0),
nrow = 10, ncol = 3, byrow = TRUE)We can see that all covariates have a positive influence on the time to event, since we find the highest values of the covariates for the highest survival times. The main function of this package, can be used to analyse this tiny dataset assuming a Cox proportional hazards model:
result <- twostagecoxph(survival.dataset, covariate.matrix,
control = twostagecoxph.control(progress = 0))print(result)
#>
#> Call:
#> twostagecoxph(survival.dataset = survival.dataset, covariate.matrix = covariate.matrix,
#> control = twostagecoxph.control(progress = 0))
#>
#> 3 covariates were marginally significant at level 0.05, resulting in 3 second stage tests.
#>
#> These are the most significant interactions found with their respective p-values:
#> 1 x 2
#> p-value epistasis 0.585673
#> amount similar results 0As we can see, the object obtained from the function has a print method associated with it, and this method provides a quick summary of the results. The object itself contains all aspects of the result, such as the vector of p-values found in the first stage, and the resulting p-values from the second stage arranged in a sparse matrix.
The function also has parallel processing capabilities, but these require a parallel back-end to be registered beforehand. Since the previous example is too small to profit from the multicore capabilities, let us introduce a toy GWAS:
str(example_survival_data)
#> 'Surv' num [1:200, 1:2] 1.42723+ 0.58487+ 0.88960 0.20406+ 3.85279+ 4.05749+ 1.01883+ 2.83626 1.46397 1.54486+ ...
#> - attr(*, "dimnames")=List of 2
#> ..$ : NULL
#> ..$ : chr [1:2] "time" "status"
#> - attr(*, "type")= chr "right"str(example_snp_data)
#> int [1:200, 1:500] 2 2 2 2 2 2 2 2 2 2 ...
#> - attr(*, "dimnames")=List of 2
#> ..$ : chr [1:200] "1" "2" "3" "4" ...
#> ..$ : chr [1:500] "SNP_1683421" "SNP_226557" "SNP_1173455" "SNP_1785590" ...As we can see, in this synthetic GWAS, we had 200 subjects of which we documented 500 SNPs. These SNPs have a strong correlation structure, as we often see in covariates derived from the genome. One pair has a build in influence on the hazard rate according to a Cox proportional hazards model. Lets register a parallel back-end and use the package to analyse this dataset:
doParallel::registerDoParallel(2)
result <- twostagecoxph(example_survival_data, example_snp_data, multicore = TRUE,
control = twostagecoxph.control(progress = 0))print(result)
#>
#> Call:
#> twostagecoxph(survival.dataset = example_survival_data, covariate.matrix = example_snp_data,
#> multicore = TRUE, control = twostagecoxph.control(progress = 0))
#>
#> 191 covariates were marginally significant at level 0.05, resulting in 18145 second stage tests.
#>
#> These are the most significant interactions found with their respective p-values:
#> SNP_263757 x SNP_222191 SNP_412807 x SNP_120042
#> p-value epistasis 0.000684222 0.004159757
#> amount similar results 1 5
#> SNP_412807 x SNP_961485 SNP_348130 x SNP_478023
#> p-value epistasis 0.007463934 0.052284475
#> amount similar results 0 0
#> SNP_2224824 x SNP_478023
#> p-value epistasis 0.065584745
#> amount similar results 5As we can see, we have some significant results. Three pairs of covariates have a corresponding p-value that is lower than 0.05, so at this level of significance we find a significant interaction. Moreover, some other pairs have the same (significant) p-values as other pairs. All of these are listed in the object:
result$most.significant.results$duplicate.interactions
#> $`SNP_263757 x SNP_222191`
#> SNP_263757 x SNP_2128600
#> 0.0006842215
#>
#> $`SNP_412807 x SNP_120042`
#> SNP_412807 x SNP_651758 SNP_412807 x SNP_1487093 SNP_412807 x SNP_840515
#> 0.004159757 0.004159757 0.004159757
#> SNP_412807 x SNP_737709 SNP_412807 x SNP_892406
#> 0.004159757 0.004159757
#>
#> $`SNP_412807 x SNP_961485`
#> named numeric(0)
#>
#> $`SNP_348130 x SNP_478023`
#> named numeric(0)
#>
#> $`SNP_2224824 x SNP_478023`
#> SNP_595526 x SNP_478023 SNP_510866 x SNP_478023 SNP_105982 x SNP_478023
#> 0.06558475 0.06558475 0.06558475
#> SNP_2143785 x SNP_478023 SNP_341004 x SNP_478023
#> 0.06558475 0.06558475This means that the interaction of these pairs are so similar to the interactions displayed in , that the p-values differ too little. These are therefore reported in a separate part of the object.
This function is written with the intention that it would be used in a Genome Wide Association Study. In such a study, the amount of covariates can rise to the millions, which gives rise to problems concerning memory. To address this, a function has been implemented to read to covariates in batches from multiple files. This way, the user has control over how much memory is required for the operation of the function. The user has to provide the files themselves, and a way for the function to read them. The large variety of data structures does not allow for an easy way for the function the adapt the files to something useable, so the user is therefore required to do some work beforehand. The files should be structured in such a way, that the user-provided function results in a matrix where each column corresponds to a covariate. Additional requirements are described in the documentation. Lets make some files from the previous example to show the operations:
number.of.covs <- dim(example_snp_data)[2]
#Split the covariate matrix into various files.
number.of.files <- 6
#500 is not a multiple of 500, so the last file has less covariates than the other ones
temp.snpfile.paths <- tempfile(rep("snpfile", number.of.files),
tmpdir = tempdir(check = TRUE),
fileext = ".txt")
indices.matrix <- matrix(c(seq_len(number.of.covs),
rep(NA, ceiling(number.of.covs/number.of.files)*
number.of.files - number.of.covs)),
ncol = number.of.files)
for(file.num in 1:number.of.files){
indices <- indices.matrix[,file.num]
write.table(example_snp_data[,indices[!is.na(indices)]],
file = temp.snpfile.paths[file.num])
}
#Associated read function: function(x) as.matrix(read.table(x))As we can see, it is possible for the last file to have less covariates than the other files. This is however the only flexibility implemented; the function assumes all other files are structured similarly, e.g. same dimensions, same presence of headers, etc.
The result from the function using batches is then the same as the one using the basic function:
print(batched.result <-
batched.twostagecoxph(example_survival_data, temp.snpfile.paths,
number.of.covariates = number.of.covs))
#> First stage complete. Commencing second stage. Analysis completed. Runtime: 9.91
#> Call:
#> batched.twostagecoxph(survival.dataset = example_survival_data,
#> covariate.filepaths = temp.snpfile.paths, number.of.covariates = number.of.covs)
#>
#> 191 covariates were marginally significant at level 0.05, resulting in 18145 second stage tests.
#>
#> These are the most significant interactions found with their respective p-values:
#> SNP_263757 x SNP_222191 SNP_412807 x SNP_120042
#> p-value epistasis 0.000684222 0.004159757
#> amount similar results 1 5
#> SNP_412807 x SNP_961485 SNP_348130 x SNP_478023
#> p-value epistasis 0.007463934 0.052284475
#> amount similar results 0 0
#> SNP_2224824 x SNP_478023
#> p-value epistasis 0.065584745
#> amount similar results 5Note that we must provide the paths to the files in a vector, and that we let the function know how many covariates can be found in all files combined. The function does allow that to be done automatically, but this requires it to read through all files individually, which costs time. The user probably know this, so it is best if it is provided. Also note that this function performs its operations in parallel, just like the previous example with .
When using this function, a consideration must be made concerning the amount of covariates in each file. The function always has no more than the contents of two files in memory on each core, but, it also needs some additional memory for fitting the model, and saving the results, etc. It is therefore recommended that no more than one quarter of the available memory of each core would be used by one file. It is also recommended that the files are as large as possible. Reading a file requires making a connection, and since the computation is done on multiple cores, this introduces some overhead. To reduce this, the files must be as large as possible so to reduce the number of connections that must be made.
This package may be too limited in its implementation for all users. Since all research is different, we have decided to not try to implement too much variation, but only provide the most basic functioning package. If a user needs additional functionality, they are highly encouraged to implement their own adaptations. In order to aid this, a highly documented version of the source code is available on the GitHub repository of this package. This should help you in avoiding the many pitfalls lurking in implementing this two stage method.