... | ... |
@@ -13,13 +13,7 @@ |
13 | 13 |
#' @param weight.function.int An integer used to assign family weights. Specifically, we use \code{weight.function.int} in a function that takes the weighted sum |
14 | 14 |
#' of the number of different SNPs and SNPs both equal to one as an argument, denoted as x, and |
15 | 15 |
#' returns a family weight equal to \code{weight.function.int}^x. Defaults to 2. |
16 |
-#' @param recessive.ref.prop The proportion to which the observed proportion of informative cases with the provisional risk genotype(s) will be compared |
|
17 |
-#' to determine whether to recode the SNP as recessive. Defaults to 0.75. |
|
18 |
-#' @param recode.test.stat For a given SNP, the minimum test statistic required to recode and recompute the fitness score using recessive coding. Defaults to 1.64. |
|
19 |
-#' See the GADGETS paper for specific details. |
|
20 |
-#' @param use.parents A integer indicating whether family level informativeness should be used alongside transmissions in computing GxE fitness scores. Defaults to 1, |
|
21 |
-#' indicating family level informativeness will be used. Specify 0 to only use transmission data. |
|
22 |
-#' @return A list of thee elements: |
|
16 |
+#' @return A list of three elements: |
|
23 | 17 |
#' \describe{ |
24 | 18 |
#' \item{pval}{The p-value of the test.} |
25 | 19 |
#' \item{obs.fitness.score}{The fitness score from the observed data} |
... | ... |
@@ -31,39 +25,39 @@ |
31 | 25 |
#' data(dad.gxe) |
32 | 26 |
#' data(mom.gxe) |
33 | 27 |
#' data(exposure) |
28 |
+#' exposure <- matrix(exposure - 1, ncol = 1, drop = FALSE) |
|
34 | 29 |
#' data(snp.annotations) |
35 | 30 |
#' pp.list <- preprocess.genetic.data(case.gxe, father.genetic.data = dad.gxe, |
36 | 31 |
#' mother.genetic.data = mom.gxe, |
37 | 32 |
#' ld.block.vec = rep(25, 4), |
38 |
-#' categorical.exposures = exposure) |
|
33 |
+#' continuous.exposures = exposure) |
|
34 |
+#' run.gadgets(pp.list, n.chromosomes = 5, chromosome.size = 3, |
|
35 |
+#' results.dir = "tmp_gxe", cluster.type = "interactive", |
|
36 |
+#' registryargs = list(file.dir = "tmp_reg_gxe", seed = 1300), |
|
37 |
+#' n.islands = 8, island.cluster.size = 4, |
|
38 |
+#' n.migrations = 1) |
|
39 | 39 |
#' |
40 |
-#' run.gadgets(pp.list, n.chromosomes = 5, chromosome.size = 3, |
|
41 |
-#' results.dir = "tmp_gxe", cluster.type = "interactive", |
|
42 |
-#' registryargs = list(file.dir = "tmp_reg_gxe", seed = 1300), |
|
43 |
-#' n.islands = 8, island.cluster.size = 4, |
|
44 |
-#' n.migrations = 1) |
|
45 | 40 |
#' |
46 | 41 |
#' combined.res <- combine.islands('tmp_gxe', snp.annotations, pp.list, 1) |
47 |
-#' null.info <- readRDS("tmp_gxe/null.mean.cov.info.rds") |
|
42 |
+#' null.info <- readRDS("tmp_gxe/null.mean.sd.info.rds") |
|
48 | 43 |
#' null.mean <- null.info[["null.mean"]] |
49 |
-#' null.cov <- null.info[["null.cov"]] |
|
44 |
+#' null.sd <- null.info[["null.sd"]] |
|
50 | 45 |
#' |
51 | 46 |
#' top.snps <- as.vector(t(combined.res[1, 1:3])) |
52 | 47 |
#' set.seed(10) |
53 |
-#' GxE.test.res <- GxE.test(top.snps, pp.list, null.mean, null.cov) |
|
48 |
+#' GxE.test.res <- GxE.test(top.snps, pp.list, null.mean, null.sd) |
|
54 | 49 |
#' |
55 | 50 |
#' unlink('tmp_gxe', recursive = TRUE) |
56 | 51 |
#' unlink('tmp_reg_gxe', recursive = TRUE) |
57 | 52 |
#' |
58 | 53 |
#' @export |
59 | 54 |
|
60 |
-GxE.test <- function(snp.cols, preprocessed.list, null.mean.vec, null.se.vec, |
|
55 |
+GxE.test <- function(snp.cols, preprocessed.list, null.mean.vec, null.sd.vec, |
|
61 | 56 |
n.permutes = 10000, n.different.snps.weight = 2, n.both.one.weight = 1, |
62 | 57 |
weight.function.int = 2) { |
63 | 58 |
|
64 | 59 |
# run the test via cpp |
65 |
- GxE_test(snp.cols, preprocessed.list, null.mean.vec, null.se.vec, n.permutes, |
|
66 |
- n.different.snps.weight, n.both.one.weight, weight.function.int, |
|
67 |
- recessive.ref.prop, recode.test.stat, use.parents) |
|
60 |
+ GxE_test(snp.cols, preprocessed.list, null.mean.vec, null.sd.vec, n.permutes, |
|
61 |
+ n.different.snps.weight, n.both.one.weight, weight.function.int) |
|
68 | 62 |
|
69 | 63 |
} |
... | ... |
@@ -181,11 +181,11 @@ epistasis_test <- function(snp_cols, preprocessed_list, n_permutes = 10000L, n_d |
181 | 181 |
.Call('_epistasisGAGE_epistasis_test', PACKAGE = 'epistasisGAGE', snp_cols, preprocessed_list, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, warn) |
182 | 182 |
} |
183 | 183 |
|
184 |
-GxE_test <- function(snp_cols, preprocessed_list, null_mean_vec, null_se_vec, n_permutes = 10000L, n_different_snps_weight = 2L, n_both_one_weight = 1L, weight_function_int = 2L, recessive_ref_prop = 0.75, recode_test_stat = 1.64, use_parents = 1L) { |
|
185 |
- .Call('_epistasisGAGE_GxE_test', PACKAGE = 'epistasisGAGE', snp_cols, preprocessed_list, null_mean_vec, null_se_vec, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, use_parents) |
|
184 |
+GxE_test <- function(target_snps, preprocessed_list, null_mean_vec, null_se_vec, n_permutes = 10000L, n_different_snps_weight = 2L, n_both_one_weight = 1L, weight_function_int = 2L) { |
|
185 |
+ .Call('_epistasisGAGE_GxE_test', PACKAGE = 'epistasisGAGE', target_snps, preprocessed_list, null_mean_vec, null_se_vec, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int) |
|
186 | 186 |
} |
187 | 187 |
|
188 |
-n2log_epistasis_pvals <- function(chromosome_list, in_data_list, null_mean_vec, null_se_vec, n_permutes = 10000L, n_different_snps_weight = 2L, n_both_one_weight = 1L, weight_function_int = 2L, recessive_ref_prop = 0.75, recode_test_stat = 1.64, GxE = FALSE, use_parents = 1L) { |
|
189 |
- .Call('_epistasisGAGE_n2log_epistasis_pvals', PACKAGE = 'epistasisGAGE', chromosome_list, in_data_list, null_mean_vec, null_se_vec, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, GxE, use_parents) |
|
188 |
+n2log_epistasis_pvals <- function(chromosome_list, preprocessed_list, n_permutes = 10000L, n_different_snps_weight = 2L, n_both_one_weight = 1L, weight_function_int = 2L, recessive_ref_prop = 0.75, recode_test_stat = 1.64, null_mean_vec_ = NULL, null_se_vec_ = NULL) { |
|
189 |
+ .Call('_epistasisGAGE_n2log_epistasis_pvals', PACKAGE = 'epistasisGAGE', chromosome_list, preprocessed_list, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, null_mean_vec_, null_se_vec_) |
|
190 | 190 |
} |
191 | 191 |
|
... | ... |
@@ -9,7 +9,7 @@ |
9 | 9 |
#' @param preprocessed.list The list output by \code{preprocess.genetic.data} run on the observed data. |
10 | 10 |
#' @param null.mean.vec A vector of null means used for comparison in the Mahalanobis distance based version of the |
11 | 11 |
#' GxE fitness score. Does not need to be specified otherwise, and can be left at its default. |
12 |
-#' @param null.se.vec A vector or estimated null standard errors for the components of the |
|
12 |
+#' @param null.sd.vec A vector or estimated null standard deviations for the components of the |
|
13 | 13 |
#' GxE fitness score. Does not need to be specified otherwise, and can be left at its default. |
14 | 14 |
#' @param score.type A character string specifying the method for aggregating SNP-pair scores across chromosome sizes. Options are |
15 | 15 |
#' 'max', 'sum', or 'logsum', defaulting to "logsum". For a given SNP-pair, it's graphical score will be the \code{score.type} of all |
... | ... |
@@ -93,7 +93,7 @@ |
93 | 93 |
compute.graphical.scores <- function(results.list, preprocessed.list, score.type = "logsum", pval.thresh = 0.05, |
94 | 94 |
n.permutes = 10000, n.different.snps.weight = 2, n.both.one.weight = 1, |
95 | 95 |
weight.function.int = 2, recessive.ref.prop = 0.75, recode.test.stat = 1.64, |
96 |
- bp.param = bpparam(), null.mean.vec = NULL, null.se.vec = NULL) { |
|
96 |
+ bp.param = bpparam(), null.mean.vec = NULL, null.sd.vec = NULL) { |
|
97 | 97 |
|
98 | 98 |
if (pval.thresh > 0.6){ |
99 | 99 |
|
... | ... |
@@ -102,9 +102,9 @@ compute.graphical.scores <- function(results.list, preprocessed.list, score.type |
102 | 102 |
} |
103 | 103 |
|
104 | 104 |
#need to specify both or neither of null mean and sd vec |
105 |
- if (is.null(null.mean.vec) & !is.null(null.se.vec) | !is.null(null.mean.vec) & is.null(null.se.vec)){ |
|
105 |
+ if (is.null(null.mean.vec) & !is.null(null.sd.vec) | !is.null(null.mean.vec) & is.null(null.sd.vec)){ |
|
106 | 106 |
|
107 |
- stop("null.mean.vec and null.se.vec must both be NULL or both not NULL") |
|
107 |
+ stop("null.mean.vec and null.sd.vec must both be NULL or both not NULL") |
|
108 | 108 |
|
109 | 109 |
} |
110 | 110 |
|
... | ... |
@@ -123,14 +123,14 @@ compute.graphical.scores <- function(results.list, preprocessed.list, score.type |
123 | 123 |
## compute graphical scores based on epistasis test |
124 | 124 |
GxE <- !is.null(preprocessed.list$exposure) |
125 | 125 |
n2log.epi.pvals <- bplapply(chrom.list, function(chrom.size.list, preprocessed.list, null.mean.vec, |
126 |
- null.se.vec, n.permutes, n.different.snps.weight, n.both.one.weight, |
|
126 |
+ null.sd.vec, n.permutes, n.different.snps.weight, n.both.one.weight, |
|
127 | 127 |
weight.function.int, recessive.ref.prop, recode.test.stat){ |
128 | 128 |
|
129 | 129 |
n2log_epistasis_pvals(chrom.size.list, preprocessed.list, n.permutes, |
130 | 130 |
n.different.snps.weight, n.both.one.weight, weight.function.int, |
131 |
- recessive.ref.prop, recode.test.stat, null.mean.vec, null.se.vec) |
|
131 |
+ recessive.ref.prop, recode.test.stat, null.mean.vec, null.sd.vec) |
|
132 | 132 |
|
133 |
- }, preprocessed.list = preprocessed.list, null.mean.vec = null.mean.vec, null.se.vec = null.se.vec, |
|
133 |
+ }, preprocessed.list = preprocessed.list, null.mean.vec = null.mean.vec, null.sd.vec = null.sd.vec, |
|
134 | 134 |
n.permutes = n.permutes, n.different.snps.weight = n.different.snps.weight, n.both.one.weight = n.both.one.weight, |
135 | 135 |
weight.function.int = weight.function.int, recessive.ref.prop = recessive.ref.prop, |
136 | 136 |
recode.test.stat = recode.test.stat, BPPARAM = bp.param) |
... | ... |
@@ -58,10 +58,10 @@ |
58 | 58 |
#' GxE fitness score. This does not need to be specified unless an analyst wants to replicate the results of a previous GADGETS |
59 | 59 |
#' GxE run, or if some of the islands of a run failed to complete, and the analyst forgot to set the seed prior to running this command. |
60 | 60 |
#' In that case, to use the same null mean vector in computing the fitness score, the analyst can find the |
61 |
-#' previously used null mean vector in the file "null.mean.se.info.rds" stored in the \code{results.dir} directory. |
|
62 |
-#' @param null.se.vec A vector of estimated null standard errors for the three components of the |
|
61 |
+#' previously used null mean vector in the file "null.mean.sd.info.rds" stored in the \code{results.dir} directory. |
|
62 |
+#' @param null.sd.vec A vector of estimated null standard errors for the three components of the |
|
63 | 63 |
#' GxE fitness score. See argument \code{null.mean.vec} for reasons this argument might be specified. For a given run, the |
64 |
-#' previously used vector can also be found in the file "null.mean.se.info.rds" stored in the \code{results.dir} directory. |
|
64 |
+#' previously used vector can also be found in the file "null.mean.sd.info.rds" stored in the \code{results.dir} directory. |
|
65 | 65 |
#' @return For each island, a list of two elements will be written to \code{results.dir}: |
66 | 66 |
#' \describe{ |
67 | 67 |
#' \item{top.chromosome.results}{A data.table of the final generation chromosomes, their fitness scores, their difference vectors, |
... | ... |
@@ -101,7 +101,7 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
101 | 101 |
generations = 500, gen.same.fitness = 50, initial.sample.duplicates = FALSE, |
102 | 102 |
snp.sampling.type = "chisq", crossover.prop = 0.8, n.islands = 1000, island.cluster.size = 4, migration.generations = 50, |
103 | 103 |
n.migrations = 20, recessive.ref.prop = 0.75, recode.test.stat = 1.64, n.random.chroms = 10000, null.mean.vec = NULL, |
104 |
- null.se.vec = NULL) { |
|
104 |
+ null.sd.vec = NULL) { |
|
105 | 105 |
|
106 | 106 |
### make sure if island clusters exist, the migration interval is set properly ### |
107 | 107 |
if (island.cluster.size > 1 & migration.generations >= generations & island.cluster.size != 1) { |
... | ... |
@@ -242,12 +242,9 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
242 | 242 |
|
243 | 243 |
} |
244 | 244 |
|
245 |
- ### if running GxE, compute the elements required for mahalanobis distance fitness score ### |
|
246 |
- #exposure <- data.list$exposure |
|
247 |
- # null.mean <- rep(0, 3) |
|
248 |
- # null.se <- rep(1, 3) |
|
245 |
+ ### if running GxE, compute the elements required for standardizing elements of the fitness score ### |
|
249 | 246 |
null.mean <- rep(0, 2) |
250 |
- null.se <- rep(1, 2) |
|
247 |
+ null.sd <- rep(1, 2) |
|
251 | 248 |
#if (!is.null(exposure)){ |
252 | 249 |
|
253 | 250 |
#storage.mode(exposure) <- "integer" |
... | ... |
@@ -255,10 +252,10 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
255 | 252 |
|
256 | 253 |
if (use.parents == 1){ |
257 | 254 |
|
258 |
- if (is.null(null.mean.vec) & is.null(null.se.vec)){ |
|
255 |
+ if (is.null(null.mean.vec) & is.null(null.sd.vec)){ |
|
259 | 256 |
|
260 | 257 |
# make sure we're not accidentally redoing this |
261 |
- out.file.name <- file.path(results.dir, "null.mean.se.info.rds") |
|
258 |
+ out.file.name <- file.path(results.dir, "null.mean.sd.info.rds") |
|
262 | 259 |
|
263 | 260 |
# # split genetic data by exposure status |
264 | 261 |
# case.genetic.data <- data.frame(data.list$case.genetic.data) |
... | ... |
@@ -292,7 +289,7 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
292 | 289 |
rep(0, 2), rep(1, 2), n.different.snps.weight, |
293 | 290 |
n.both.one.weight) |
294 | 291 |
null.mean <- colMeans(null.vec.mat) |
295 |
- null.se <- sqrt(diag(cov(null.vec.mat))) |
|
292 |
+ null.sd <- sqrt(diag(cov(null.vec.mat))) |
|
296 | 293 |
|
297 | 294 |
#save these if needed later |
298 | 295 |
if (file.exists(out.file.name)){ |
... | ... |
@@ -317,10 +314,10 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
317 | 314 |
} |
318 | 315 |
|
319 | 316 |
|
320 |
- } else if (!is.null(null.mean.vec) & !is.null(null.se.vec)){ |
|
317 |
+ } else if (!is.null(null.mean.vec) & !is.null(null.sd.vec)){ |
|
321 | 318 |
|
322 | 319 |
# make sure we're not accidentally redoing this |
323 |
- out.file.name <- file.path(results.dir, "null.mean.se.info.rds") |
|
320 |
+ out.file.name <- file.path(results.dir, "null.mean.sd.info.rds") |
|
324 | 321 |
if (file.exists(out.file.name)){ |
325 | 322 |
|
326 | 323 |
prev.res <- readRDS(out.file.name) |
... | ... |
@@ -329,18 +326,18 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
329 | 326 |
|
330 | 327 |
if (any(null.mean != prev.mean) | any(prev.se != null.se)){ |
331 | 328 |
|
332 |
- stop(paste("null mean and/or se vector do not match the previous values in", |
|
329 |
+ stop(paste("null mean and/or sd vector do not match the previous values in", |
|
333 | 330 |
out.file.name)) |
334 | 331 |
|
335 | 332 |
} |
336 | 333 |
null.mean <- null.mean.vec |
337 |
- null.se <- null.se.vec |
|
334 |
+ null.sd <- null.sd.vec |
|
338 | 335 |
|
339 | 336 |
} |
340 | 337 |
|
341 | 338 |
} else { |
342 | 339 |
|
343 |
- stop("null.mean.vec and null.se.vec must both be null or both not be null") |
|
340 |
+ stop("null.mean.vec and null.sd.vec must both be null or both not be null") |
|
344 | 341 |
|
345 | 342 |
} |
346 | 343 |
|
... | ... |
@@ -404,21 +401,12 @@ run.gadgets <- function(data.list, n.chromosomes, chromosome.size, results.dir, |
404 | 401 |
} |
405 | 402 |
|
406 | 403 |
# write jobs to registry |
407 |
- # ids <- batchMap(GADGETS, cluster.number = cluster.ids, more.args = list(results.dir = results.dir, n.migrations = n.migrations, |
|
408 |
- # case.genetic.data = data.list$case.genetic.data, complement.genetic.data = data.list$complement.genetic.data, |
|
409 |
- # ld.block.vec = data.list$ld.block.vec, n.chromosomes = n.chromosomes, chromosome.size = chromosome.size, snp.chisq = snp.chisq, |
|
410 |
- # weight.lookup = weight.lookup, null.mean.vec = null.mean, null.se.vec = null.se, island.cluster.size = island.cluster.size, |
|
411 |
- # n.different.snps.weight = n.different.snps.weight, n.both.one.weight = n.both.one.weight, migration.interval = migration.generations, |
|
412 |
- # gen.same.fitness = gen.same.fitness, max.generations = generations, initial.sample.duplicates = initial.sample.duplicates, crossover.prop = crossover.prop, |
|
413 |
- # recessive.ref.prop = recessive.ref.prop, recode.test.stat = recode.test.stat, exposure.levels = data.list$exposure.levels, |
|
414 |
- # exposure = exposure, use.parents = use.parents), |
|
415 |
- # reg = registry) |
|
416 | 404 |
ids <- batchMap(GADGETS, cluster.number = cluster.ids, |
417 | 405 |
more.args = list(results.dir = results.dir, n.migrations = n.migrations, |
418 | 406 |
case.genetic.data = case.genetic.data, complement.genetic.data = complement.genetic.data, case.genetic.data.n = case.genetic.data.n, |
419 | 407 |
complement.genetic.data.n = complement.genetic.data.n, exposure.mat = exposure.mat, weight.lookup.n = weight.lookup.n, |
420 | 408 |
ld.block.vec = data.list$ld.block.vec, n.chromosomes = n.chromosomes, chromosome.size = chromosome.size, snp.chisq = snp.chisq, |
421 |
- weight.lookup = weight.lookup, null.mean.vec = null.mean, null.se.vec = null.se, island.cluster.size = island.cluster.size, |
|
409 |
+ weight.lookup = weight.lookup, null.mean.vec = null.mean, null.se.vec = null.sd, island.cluster.size = island.cluster.size, |
|
422 | 410 |
n.different.snps.weight = n.different.snps.weight, n.both.one.weight = n.both.one.weight, migration.interval = migration.generations, |
423 | 411 |
gen.same.fitness = gen.same.fitness, max.generations = generations, initial.sample.duplicates = initial.sample.duplicates, |
424 | 412 |
crossover.prop = crossover.prop, recessive.ref.prop = recessive.ref.prop, recode.test.stat = recode.test.stat, |
... | ... |
@@ -8,14 +8,11 @@ GxE.test( |
8 | 8 |
snp.cols, |
9 | 9 |
preprocessed.list, |
10 | 10 |
null.mean.vec, |
11 |
- null.se.vec, |
|
11 |
+ null.sd.vec, |
|
12 | 12 |
n.permutes = 10000, |
13 | 13 |
n.different.snps.weight = 2, |
14 | 14 |
n.both.one.weight = 1, |
15 |
- weight.function.int = 2, |
|
16 |
- recessive.ref.prop = 0.75, |
|
17 |
- recode.test.stat = 1.64, |
|
18 |
- use.parents = 1 |
|
15 |
+ weight.function.int = 2 |
|
19 | 16 |
) |
20 | 17 |
} |
21 | 18 |
\arguments{ |
... | ... |
@@ -34,18 +31,9 @@ is multiplied in computing the family weights. Defaults to 1.} |
34 | 31 |
\item{weight.function.int}{An integer used to assign family weights. Specifically, we use \code{weight.function.int} in a function that takes the weighted sum |
35 | 32 |
of the number of different SNPs and SNPs both equal to one as an argument, denoted as x, and |
36 | 33 |
returns a family weight equal to \code{weight.function.int}^x. Defaults to 2.} |
37 |
- |
|
38 |
-\item{recessive.ref.prop}{The proportion to which the observed proportion of informative cases with the provisional risk genotype(s) will be compared |
|
39 |
-to determine whether to recode the SNP as recessive. Defaults to 0.75.} |
|
40 |
- |
|
41 |
-\item{recode.test.stat}{For a given SNP, the minimum test statistic required to recode and recompute the fitness score using recessive coding. Defaults to 1.64. |
|
42 |
-See the GADGETS paper for specific details.} |
|
43 |
- |
|
44 |
-\item{use.parents}{A integer indicating whether family level informativeness should be used alongside transmissions in computing GxE fitness scores. Defaults to 1, |
|
45 |
-indicating family level informativeness will be used. Specify 0 to only use transmission data.} |
|
46 | 34 |
} |
47 | 35 |
\value{ |
48 |
-A list of thee elements: |
|
36 |
+A list of three elements: |
|
49 | 37 |
\describe{ |
50 | 38 |
\item{pval}{The p-value of the test.} |
51 | 39 |
\item{obs.fitness.score}{The fitness score from the observed data} |
... | ... |
@@ -62,26 +50,27 @@ data(case.gxe) |
62 | 50 |
data(dad.gxe) |
63 | 51 |
data(mom.gxe) |
64 | 52 |
data(exposure) |
53 |
+exposure <- matrix(exposure - 1, ncol = 1, drop = FALSE) |
|
65 | 54 |
data(snp.annotations) |
66 | 55 |
pp.list <- preprocess.genetic.data(case.gxe, father.genetic.data = dad.gxe, |
67 | 56 |
mother.genetic.data = mom.gxe, |
68 | 57 |
ld.block.vec = rep(25, 4), |
69 |
- categorical.exposures = exposure) |
|
58 |
+ continuous.exposures = exposure) |
|
59 |
+ run.gadgets(pp.list, n.chromosomes = 5, chromosome.size = 3, |
|
60 |
+ results.dir = "tmp_gxe", cluster.type = "interactive", |
|
61 |
+ registryargs = list(file.dir = "tmp_reg_gxe", seed = 1300), |
|
62 |
+ n.islands = 8, island.cluster.size = 4, |
|
63 |
+ n.migrations = 1) |
|
70 | 64 |
|
71 |
-run.gadgets(pp.list, n.chromosomes = 5, chromosome.size = 3, |
|
72 |
- results.dir = "tmp_gxe", cluster.type = "interactive", |
|
73 |
- registryargs = list(file.dir = "tmp_reg_gxe", seed = 1300), |
|
74 |
- n.islands = 8, island.cluster.size = 4, |
|
75 |
- n.migrations = 1) |
|
76 | 65 |
|
77 | 66 |
combined.res <- combine.islands('tmp_gxe', snp.annotations, pp.list, 1) |
78 |
-null.info <- readRDS("tmp_gxe/null.mean.cov.info.rds") |
|
67 |
+null.info <- readRDS("tmp_gxe/null.mean.sd.info.rds") |
|
79 | 68 |
null.mean <- null.info[["null.mean"]] |
80 |
-null.cov <- null.info[["null.cov"]] |
|
69 |
+null.sd <- null.info[["null.sd"]] |
|
81 | 70 |
|
82 | 71 |
top.snps <- as.vector(t(combined.res[1, 1:3])) |
83 | 72 |
set.seed(10) |
84 |
-GxE.test.res <- GxE.test(top.snps, pp.list, null.mean, null.cov) |
|
73 |
+GxE.test.res <- GxE.test(top.snps, pp.list, null.mean, null.sd) |
|
85 | 74 |
|
86 | 75 |
unlink('tmp_gxe', recursive = TRUE) |
87 | 76 |
unlink('tmp_reg_gxe', recursive = TRUE) |
... | ... |
@@ -7,8 +7,6 @@ |
7 | 7 |
compute.graphical.scores( |
8 | 8 |
results.list, |
9 | 9 |
preprocessed.list, |
10 |
- null.mean.vec = rep(0, 3), |
|
11 |
- null.se.vec = rep(1, 3), |
|
12 | 10 |
score.type = "logsum", |
13 | 11 |
pval.thresh = 0.05, |
14 | 12 |
n.permutes = 10000, |
... | ... |
@@ -17,7 +15,9 @@ compute.graphical.scores( |
17 | 15 |
weight.function.int = 2, |
18 | 16 |
recessive.ref.prop = 0.75, |
19 | 17 |
recode.test.stat = 1.64, |
20 |
- bp.param = bpparam() |
|
18 |
+ bp.param = bpparam(), |
|
19 |
+ null.mean.vec = NULL, |
|
20 |
+ null.sd.vec = NULL |
|
21 | 21 |
) |
22 | 22 |
} |
23 | 23 |
\arguments{ |
... | ... |
@@ -28,12 +28,6 @@ otherwise an error will be returned.} |
28 | 28 |
|
29 | 29 |
\item{preprocessed.list}{The list output by \code{preprocess.genetic.data} run on the observed data.} |
30 | 30 |
|
31 |
-\item{null.mean.vec}{A vector of null means used for comparison in the Mahalanobis distance based version of the |
|
32 |
-GxE fitness score. Does not need to be specified otherwise, and can be left at its default.} |
|
33 |
- |
|
34 |
-\item{null.se.vec}{A vector or estimated null standard errors for the components of the |
|
35 |
-GxE fitness score. Does not need to be specified otherwise, and can be left at its default.} |
|
36 |
- |
|
37 | 31 |
\item{score.type}{A character string specifying the method for aggregating SNP-pair scores across chromosome sizes. Options are |
38 | 32 |
'max', 'sum', or 'logsum', defaulting to "logsum". For a given SNP-pair, it's graphical score will be the \code{score.type} of all |
39 | 33 |
graphical scores of chromosomes containing that pair across chromosome sizes. The choice of 'logsum' rather than 'sum' |
... | ... |
@@ -64,6 +58,12 @@ to determine whether to recode the SNP as recessive. Defaults to 0.75.} |
64 | 58 |
See the GADGETS paper for specific details.} |
65 | 59 |
|
66 | 60 |
\item{bp.param}{The BPPARAM argument to be passed to bplapply. See \code{BiocParallel::bplapply} for more details.} |
61 |
+ |
|
62 |
+\item{null.mean.vec}{A vector of null means used for comparison in the Mahalanobis distance based version of the |
|
63 |
+GxE fitness score. Does not need to be specified otherwise, and can be left at its default.} |
|
64 |
+ |
|
65 |
+\item{null.sd.vec}{A vector or estimated null standard deviations for the components of the |
|
66 |
+GxE fitness score. Does not need to be specified otherwise, and can be left at its default.} |
|
67 | 67 |
} |
68 | 68 |
\value{ |
69 | 69 |
A list of two elements: |
... | ... |
@@ -31,7 +31,7 @@ run.gadgets( |
31 | 31 |
recode.test.stat = 1.64, |
32 | 32 |
n.random.chroms = 10000, |
33 | 33 |
null.mean.vec = NULL, |
34 |
- null.se.vec = NULL |
|
34 |
+ null.sd.vec = NULL |
|
35 | 35 |
) |
36 | 36 |
} |
37 | 37 |
\arguments{ |
... | ... |
@@ -114,11 +114,11 @@ compute the GxE fitness score.} |
114 | 114 |
GxE fitness score. This does not need to be specified unless an analyst wants to replicate the results of a previous GADGETS |
115 | 115 |
GxE run, or if some of the islands of a run failed to complete, and the analyst forgot to set the seed prior to running this command. |
116 | 116 |
In that case, to use the same null mean vector in computing the fitness score, the analyst can find the |
117 |
-previously used null mean vector in the file "null.mean.se.info.rds" stored in the \code{results.dir} directory.} |
|
117 |
+previously used null mean vector in the file "null.mean.sd.info.rds" stored in the \code{results.dir} directory.} |
|
118 | 118 |
|
119 |
-\item{null.se.vec}{A vector of estimated null standard errors for the three components of the |
|
119 |
+\item{null.sd.vec}{A vector of estimated null standard errors for the three components of the |
|
120 | 120 |
GxE fitness score. See argument \code{null.mean.vec} for reasons this argument might be specified. For a given run, the |
121 |
-previously used vector can also be found in the file "null.mean.se.info.rds" stored in the \code{results.dir} directory.} |
|
121 |
+previously used vector can also be found in the file "null.mean.sd.info.rds" stored in the \code{results.dir} directory.} |
|
122 | 122 |
|
123 | 123 |
\item{use.parents}{A logical indicating whether parent data should be used in computing the fitness score. Defaults to FALSE. This should only be set to true |
124 | 124 |
if the population is homogenous with no exposure related population structure.} |
... | ... |
@@ -770,12 +770,12 @@ BEGIN_RCPP |
770 | 770 |
END_RCPP |
771 | 771 |
} |
772 | 772 |
// GxE_test |
773 |
-List GxE_test(IntegerVector snp_cols, List preprocessed_list, arma::vec null_mean_vec, arma::vec null_se_vec, int n_permutes, int n_different_snps_weight, int n_both_one_weight, int weight_function_int, double recessive_ref_prop, double recode_test_stat, int use_parents); |
|
774 |
-RcppExport SEXP _epistasisGAGE_GxE_test(SEXP snp_colsSEXP, SEXP preprocessed_listSEXP, SEXP null_mean_vecSEXP, SEXP null_se_vecSEXP, SEXP n_permutesSEXP, SEXP n_different_snps_weightSEXP, SEXP n_both_one_weightSEXP, SEXP weight_function_intSEXP, SEXP recessive_ref_propSEXP, SEXP recode_test_statSEXP, SEXP use_parentsSEXP) { |
|
773 |
+List GxE_test(arma::uvec target_snps, List preprocessed_list, arma::vec null_mean_vec, arma::vec null_se_vec, int n_permutes, int n_different_snps_weight, int n_both_one_weight, int weight_function_int); |
|
774 |
+RcppExport SEXP _epistasisGAGE_GxE_test(SEXP target_snpsSEXP, SEXP preprocessed_listSEXP, SEXP null_mean_vecSEXP, SEXP null_se_vecSEXP, SEXP n_permutesSEXP, SEXP n_different_snps_weightSEXP, SEXP n_both_one_weightSEXP, SEXP weight_function_intSEXP) { |
|
775 | 775 |
BEGIN_RCPP |
776 | 776 |
Rcpp::RObject rcpp_result_gen; |
777 | 777 |
Rcpp::RNGScope rcpp_rngScope_gen; |
778 |
- Rcpp::traits::input_parameter< IntegerVector >::type snp_cols(snp_colsSEXP); |
|
778 |
+ Rcpp::traits::input_parameter< arma::uvec >::type target_snps(target_snpsSEXP); |
|
779 | 779 |
Rcpp::traits::input_parameter< List >::type preprocessed_list(preprocessed_listSEXP); |
780 | 780 |
Rcpp::traits::input_parameter< arma::vec >::type null_mean_vec(null_mean_vecSEXP); |
781 | 781 |
Rcpp::traits::input_parameter< arma::vec >::type null_se_vec(null_se_vecSEXP); |
... | ... |
@@ -783,32 +783,27 @@ BEGIN_RCPP |
783 | 783 |
Rcpp::traits::input_parameter< int >::type n_different_snps_weight(n_different_snps_weightSEXP); |
784 | 784 |
Rcpp::traits::input_parameter< int >::type n_both_one_weight(n_both_one_weightSEXP); |
785 | 785 |
Rcpp::traits::input_parameter< int >::type weight_function_int(weight_function_intSEXP); |
786 |
- Rcpp::traits::input_parameter< double >::type recessive_ref_prop(recessive_ref_propSEXP); |
|
787 |
- Rcpp::traits::input_parameter< double >::type recode_test_stat(recode_test_statSEXP); |
|
788 |
- Rcpp::traits::input_parameter< int >::type use_parents(use_parentsSEXP); |
|
789 |
- rcpp_result_gen = Rcpp::wrap(GxE_test(snp_cols, preprocessed_list, null_mean_vec, null_se_vec, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, use_parents)); |
|
786 |
+ rcpp_result_gen = Rcpp::wrap(GxE_test(target_snps, preprocessed_list, null_mean_vec, null_se_vec, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int)); |
|
790 | 787 |
return rcpp_result_gen; |
791 | 788 |
END_RCPP |
792 | 789 |
} |
793 | 790 |
// n2log_epistasis_pvals |
794 |
-NumericVector n2log_epistasis_pvals(ListOf<IntegerVector> chromosome_list, List in_data_list, arma::vec null_mean_vec, arma::vec null_se_vec, int n_permutes, int n_different_snps_weight, int n_both_one_weight, int weight_function_int, double recessive_ref_prop, double recode_test_stat, bool GxE, int use_parents); |
|
795 |
-RcppExport SEXP _epistasisGAGE_n2log_epistasis_pvals(SEXP chromosome_listSEXP, SEXP in_data_listSEXP, SEXP null_mean_vecSEXP, SEXP null_se_vecSEXP, SEXP n_permutesSEXP, SEXP n_different_snps_weightSEXP, SEXP n_both_one_weightSEXP, SEXP weight_function_intSEXP, SEXP recessive_ref_propSEXP, SEXP recode_test_statSEXP, SEXP GxESEXP, SEXP use_parentsSEXP) { |
|
791 |
+NumericVector n2log_epistasis_pvals(ListOf<IntegerVector> chromosome_list, List preprocessed_list, int n_permutes, int n_different_snps_weight, int n_both_one_weight, int weight_function_int, double recessive_ref_prop, double recode_test_stat, Nullable<NumericVector> null_mean_vec_, Nullable<NumericVector> null_se_vec_); |
|
792 |
+RcppExport SEXP _epistasisGAGE_n2log_epistasis_pvals(SEXP chromosome_listSEXP, SEXP preprocessed_listSEXP, SEXP n_permutesSEXP, SEXP n_different_snps_weightSEXP, SEXP n_both_one_weightSEXP, SEXP weight_function_intSEXP, SEXP recessive_ref_propSEXP, SEXP recode_test_statSEXP, SEXP null_mean_vec_SEXP, SEXP null_se_vec_SEXP) { |
|
796 | 793 |
BEGIN_RCPP |
797 | 794 |
Rcpp::RObject rcpp_result_gen; |
798 | 795 |
Rcpp::RNGScope rcpp_rngScope_gen; |
799 | 796 |
Rcpp::traits::input_parameter< ListOf<IntegerVector> >::type chromosome_list(chromosome_listSEXP); |
800 |
- Rcpp::traits::input_parameter< List >::type in_data_list(in_data_listSEXP); |
|
801 |
- Rcpp::traits::input_parameter< arma::vec >::type null_mean_vec(null_mean_vecSEXP); |
|
802 |
- Rcpp::traits::input_parameter< arma::vec >::type null_se_vec(null_se_vecSEXP); |
|
797 |
+ Rcpp::traits::input_parameter< List >::type preprocessed_list(preprocessed_listSEXP); |
|
803 | 798 |
Rcpp::traits::input_parameter< int >::type n_permutes(n_permutesSEXP); |
804 | 799 |
Rcpp::traits::input_parameter< int >::type n_different_snps_weight(n_different_snps_weightSEXP); |
805 | 800 |
Rcpp::traits::input_parameter< int >::type n_both_one_weight(n_both_one_weightSEXP); |
806 | 801 |
Rcpp::traits::input_parameter< int >::type weight_function_int(weight_function_intSEXP); |
807 | 802 |
Rcpp::traits::input_parameter< double >::type recessive_ref_prop(recessive_ref_propSEXP); |
808 | 803 |
Rcpp::traits::input_parameter< double >::type recode_test_stat(recode_test_statSEXP); |
809 |
- Rcpp::traits::input_parameter< bool >::type GxE(GxESEXP); |
|
810 |
- Rcpp::traits::input_parameter< int >::type use_parents(use_parentsSEXP); |
|
811 |
- rcpp_result_gen = Rcpp::wrap(n2log_epistasis_pvals(chromosome_list, in_data_list, null_mean_vec, null_se_vec, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, GxE, use_parents)); |
|
804 |
+ Rcpp::traits::input_parameter< Nullable<NumericVector> >::type null_mean_vec_(null_mean_vec_SEXP); |
|
805 |
+ Rcpp::traits::input_parameter< Nullable<NumericVector> >::type null_se_vec_(null_se_vec_SEXP); |
|
806 |
+ rcpp_result_gen = Rcpp::wrap(n2log_epistasis_pvals(chromosome_list, preprocessed_list, n_permutes, n_different_snps_weight, n_both_one_weight, weight_function_int, recessive_ref_prop, recode_test_stat, null_mean_vec_, null_se_vec_)); |
|
812 | 807 |
return rcpp_result_gen; |
813 | 808 |
END_RCPP |
814 | 809 |
} |
... | ... |
@@ -859,8 +854,8 @@ static const R_CallMethodDef CallEntries[] = { |
859 | 854 |
{"_epistasisGAGE_epistasis_test_permute", (DL_FUNC) &_epistasisGAGE_epistasis_test_permute, 10}, |
860 | 855 |
{"_epistasisGAGE_epistasis_test_null_scores", (DL_FUNC) &_epistasisGAGE_epistasis_test_null_scores, 11}, |
861 | 856 |
{"_epistasisGAGE_epistasis_test", (DL_FUNC) &_epistasisGAGE_epistasis_test, 9}, |
862 |
- {"_epistasisGAGE_GxE_test", (DL_FUNC) &_epistasisGAGE_GxE_test, 11}, |
|
863 |
- {"_epistasisGAGE_n2log_epistasis_pvals", (DL_FUNC) &_epistasisGAGE_n2log_epistasis_pvals, 12}, |
|
857 |
+ {"_epistasisGAGE_GxE_test", (DL_FUNC) &_epistasisGAGE_GxE_test, 8}, |
|
858 |
+ {"_epistasisGAGE_n2log_epistasis_pvals", (DL_FUNC) &_epistasisGAGE_n2log_epistasis_pvals, 10}, |
|
864 | 859 |
{NULL, NULL, 0} |
865 | 860 |
}; |
866 | 861 |
|