R/GADGETS.R
In mnodzenski/epistasisGA: An R package to identify multi-snp effects in nuclear family studies using the GADGETS method
Defines functions
GADGETS
Documented in
GADGETS
#' A function to run the GADGETS method
#'
#' This function runs the GADGETS method on a given cluster of islands.
#' It is a wrapper for
#' the underlying Rcpp function run_GADGETS.
#'
#' @param cluster.number An integer indicating the cluster number
#' (used for labeling the output file).
#' @param results.dir The directory to which island results will be saved.
#' @param case.genetic.data The genetic data of the disease affected children
#' from case-parent trios or disease-discordant sibling pairs. If searching for
#' maternal SNPs that are related to risk of disease in the child, some of the
#' columns in \code{case.genetic.data} may contain maternal SNP genotypes
#' (See argument \code{mother.snps} for how to indicate which SNPs columns
#' correspond to maternal genotypes). Columns are SNP allele counts, and rows
#' are individuals.This object should be of class 'matrix'. The ordering of the
#' columns must be consistent with the LD structure specified in
#' \code{ld.block.vec}. The genotypes cannot be dosages imputed with
#' uncertainty. If any data are missing for a particular family member for a
#' particular SNP, that SNP's genotype should be coded as -9 for each member of
#' the entire family: (\code{case.genetic.data} and
#' \code{father.genetic.data}/\code{mother.genetic.data},
#' or \code{case.genetic.data} and \code{complement.genetic.data}). If running
#' the experimental E-GADGETS method, this argument should be set to a 1x1
#' matrix whose only value is 0.0, and \code{case.genetic.data.n} will be used
#' to specify the case genotypes.
#' @param complement.genetic.data A genetic dataset for the controls
#' corresponding to the genotypes in \code{case.genetic.data}.For SNPs that
#' correspond to the affected child in \code{case.genetic.data}, the
#' corresponding column in \code{complement.genetic.data} should be set equal to
#' mother allele count + father allele count - case allele count. If using
#' disease-discordant siblings this argument should be the genotypes for the
#' unaffected siblings. For SNPs in \code{case.genetic.data} that represent
#' maternal genotypes (if any) the corresponding column in
#' \code{complement.genetic.data} should be the paternal genotypes for that SNP.
#' This object should be of class 'matrix'. Columns are SNP allele counts,
#' rows are families. If not specified, \code{father.genetic.data} and
#' \code{mother.genetic.data} must be specified. The genotypes cannot be dosages
#' imputed with uncertainty. If any data are missing for a particular family for
#' a particular SNP, that SNP's genotype should be coded as -9 for the entire
#' family (\code{case.genetic.data} and \code{complement.genetic.data}) for that
#' SNP. If running the experimental E-GADGETS method, this argument should be
#' set to a 1x1 matrix whose only value is 0.0, and
#' \code{complement.genetic.data.n} will be used to specify the complement
#' genotypes.
#' @param case.genetic.data.n (experimental) A matrix, to be used in the
#' experimental E-GADGETS method, containing the same data as described above
#' for \code{case.genetic.data}, but the genotypes here are stored as floating
#' point values, as opposed to integer values in \code{case.genetic.data}. If
#' not running E-GADGETS, this should be specified as a 1x1 matrix whose only
#' value is 0.0.
#' @param mother.genetic.data.n (experimental) A matrix, to be used in the
#' experimental E-GADGETS method, containing the genotypes for the mothers of
#' \code{case.genetic.data}, where the genotypes are stored as
#' floating point values, as opposed to integer values. If not running
#' E-GADGETS, this should be specified as a 1x1 matrix whose only value is 0.0.
#' @param father.genetic.data.n (experimental) A matrix, to be used in the
#' experimental E-GADGETS method, containing the genotypes for the fathers of
#' \code{case.genetic.data}, where the genotypes are stored as
#' floating point values, as opposed to integer values. If not running
#' E-GADGETS, this should be specified as a 1x1 matrix whose only value is 0.0.
#' @param exposure.mat (experimental) A matrix of the input categorical
#' and continuous exposures, if specified, to be used in the experimental
#' E-GADGETS method. If not running E-GADGETS, this should be a 1x1 matrix whose
#' only entry is 0.0.
#' @param weight.lookup.n (experimental) A vector that maps a family weight to
#' the weighted sum of the number of different SNPs and SNPs both equal to one,
#' to be used by the experimetnal E-GADGETS method. The vector should store
#' values as floating point values, not integers. If not running E-GADGETS, this
#' argument should be specified as 0.0, and will not be used in the GA. Instead,
#' for GADGETS, computation of the family weights will be based on argument
#' \code{weight.lookup}, which is computed in the same way, except stores values
#' as integers.
#' @param ld.block.vec An integer vector specifying the linkage blocks of the
#' input SNPs. As an example, for 100 candidate SNPs, suppose we specify
#' \code{ld.block.vec <- c(25, 75, 100)}. This vector indicates that the input
#' genetic data has 3 distinct linkage blocks, with SNPs 1-25 in the first
#' linkage block, 26-75 in the second block, and 76-100 in the third block. Note
#' that this means the ordering of the columns (SNPs) in
#' \code{case.genetic.data} must be consistent with the LD blocks specified in
#' \code{ld.block.vec}. In the absence of outside information, a reasonable
#' default is to consider SNPs to be in LD if they are located on the same
#' biological chromosome. If \code{case.genetic.data} includes both maternal
#' and child SNP genotypes, we recommend considering any maternal SNP and any
#' child SNP located on the same nominal biological chromosome as 'in linkage'.
#' E.g., we recommend considering any maternal SNPs located on chromosome 1
#' as being 'linked' to any child SNPs located on chromosome 1, even though,
#' strictly speaking, the maternal and child SNPs are located on separate pieces
#' of DNA. If running E-GADGETS, this argument should be specified as 0, and
#' will not be used.
#' @param n.chromosomes An integer specifying the number of chromosomes to use
#' in the GA.
#' @param chromosome.size An integer specifying the number of SNPs on each
#' chromosome.
#' @param snp.chisq A vector of statistics to be used in sampling SNPs for
#' mutation. By default, these are the square roots of
#' the chi-square marginal SNP-disease association statistics for each column in
#' \code{case.genetic.data}, but can also be manually specified or uniformly 1
#' (corresponding to totally random sampling).
#' @param weight.lookup A vector that maps a family weight to the weighted sum
#' of the number of different SNPs and SNPs both equal to one. This should
#' store values as integers.
#' @param null.mean.vec (experimental) A vector of estimated null means for each
#' component of the E-GADGETS (GxGxE) fitness score. For all other uses, this
#' should be specified as rep(0, 2) and will not be used.
#' @param null.se.vec (experimental) A vector of estimated null standard
#' deviations for each component of the E-GADGETS (GxGxE) fitness score. For all
#' other uses, this should be specified as rep(0, 2) and will not be used.
#' @param island.cluster.size An integer specifying the number of islands in
#' the cluster. See code{run.gadgets} for additional details.
#' @param n.migrations The number of chromosomes that migrate among islands.
#' This value must be less than \code{n.chromosomes} and greater than 0,
#' defaulting to 20.
#' @param n.different.snps.weight The number by which the number of different
#' SNPs between a case and complement is multiplied in computing the family
#' weights. Defaults to 2.
#' @param n.both.one.weight The number by which the number of SNPs equal to 1 in
#' both the case and complement is multiplied in computing the family weights.
#' Defaults to 1.
#' @param migration.interval The interval of generations for which GADGETS will
#' run prior to migration of top chromosomes among islands in a cluster.
#' Defaults to 50. In other words, top chromosomes will migrate among cluster
#' islands every \code{migration.interval} generations. We also check for
#' convergence at each of these intervals.
#' @param gen.same.fitness The number of consecutive generations with the same
#' fitness score required for algorithm termination. Defaults to 50.
#' @param max.generations The maximum number of generations for which GADGETS
#' will run. Defaults to 500.
#' @param initial.sample.duplicates A logical indicating whether the same SNP
#' can appear in more than one chromosome in the initial sample of chromosomes
#' (the same SNP may appear in more than one chromosome thereafter, regardless).
#' Defaults to FALSE.
#' @param crossover.prop A numeric between 0 and 1 indicating the proportion of
#' chromosomes to be subjected to cross over. The remaining proportion will be
#' mutated. Defaults to 0.8.
#' @param recessive.ref.prop The proportion to which the observed proportion of
#' informative cases with the provisional risk genotype(s) will be compared
#' to determine whether to recode the SNP as recessive. Defaults to 0.75.
#' @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. See the GADGETS paper for specific details.
#' @param E_GADGETS (experimental) A boolean indicating whether to run the
#' experimental 'E_GADGETS' method.
#' @return For each island in the cluster, an rds object containing a list with
#' the following elements will be written to \code{results.dir}.
#' \describe{
#' \item{top.chromosome.results}{A data.table of the final generation
#' chromosomes, their fitness scores, and, for GADGETS, additional information
#' pertaining to nominated risk-related genotypes. See the package vignette for
#' an example and for additional details.}
#' \item{n.generations}{The total number of generations run.}
#' }
#'
#' @examples
#'
#' set.seed(10)
#' data(case)
#' case <- as.matrix(case)
#' data(dad)
#' dad <- as.matrix(dad)
#' data(mom)
#' mom <- as.matrix(mom)
#' data.list <- preprocess.genetic.data(case[, 1:10],
#' father.genetic.data = dad[ , 1:10],
#' mother.genetic.data = mom[ , 1:10],
#' ld.block.vec = c(10))
#'
#' chisq.stats <- sqrt(data.list$chisq.stats)
#' ld.block.vec <- data.list$ld.block.vec
#' case.genetic.data <- data.list$case.genetic.data
#' complement.genetic.data <- data.list$complement.genetic.data
#'
#' #required inputs but not actually used in function below
#' case.genetic.data.n <- matrix(0.0, 1, 1)
#' mother.genetic.data.n <- matrix(0.0, 1, 1)
#' father.genetic.data.n <- matrix(0.0, 1, 1)
#' exposure.mat <- data.list$exposure.mat + 0.0
#'
#' weight.lookup <- vapply(seq_len(6), function(x) 2^x, 1)
#' dir.create('tmp')
#' GADGETS(cluster.number = 1, results.dir = 'tmp',
#' case.genetic.data = case.genetic.data,
#' complement.genetic.data = complement.genetic.data,
#' case.genetic.data.n = case.genetic.data.n,
#' mother.genetic.data.n = mother.genetic.data.n,
#' father.genetic.data.n = father.genetic.data.n,
#' exposure.mat = exposure.mat,
#' weight.lookup.n = weight.lookup + 0.0,
#' ld.block.vec = ld.block.vec,
#' n.chromosomes = 10, chromosome.size = 3, snp.chisq = chisq.stats,
#' weight.lookup = weight.lookup, n.migrations = 2,
#' migration.interval = 5,
#' gen.same.fitness = 10, max.generations = 10)
#'
#' @importFrom data.table as.data.table setorder setDT rbindlist transpose
#' @useDynLib epistasisGA
#' @export
GADGETS <- function(cluster.number, results.dir, case.genetic.data,
complement.genetic.data, case.genetic.data.n,
mother.genetic.data.n, father.genetic.data.n, exposure.mat,
weight.lookup.n, ld.block.vec, n.chromosomes,
chromosome.size, snp.chisq,
weight.lookup, null.mean.vec = c(0, 0),
null.se.vec = c(1, 1), island.cluster.size = 4,
n.migrations = 20, n.different.snps.weight = 2,
n.both.one.weight = 1, migration.interval = 50,
gen.same.fitness = 50, max.generations = 500,
initial.sample.duplicates = FALSE, crossover.prop = 0.8,
recessive.ref.prop = 0.75, recode.test.stat = 1.64,
E_GADGETS = FALSE) {
### run rcpp version of GADGETS ##
rcpp.res <- run_GADGETS(island.cluster.size, n.migrations, ld.block.vec,
n.chromosomes, chromosome.size, weight.lookup,
snp.chisq, case.genetic.data,
complement.genetic.data, case.genetic.data.n,
mother.genetic.data.n, father.genetic.data.n,
exposure.mat, weight.lookup.n, null.mean.vec,
null.se.vec, n.different.snps.weight,
n.both.one.weight, migration.interval,
gen.same.fitness, max.generations,
initial.sample.duplicates, crossover.prop,
recessive.ref.prop, recode.test.stat,
E_GADGETS)
### clean up and output results
lapply(seq_along(rcpp.res), function(island.number) {
#pick out the pieces from rcpp output
rcpp.res.length <- length(rcpp.res[[island.number]])
n.generations <- rcpp.res[[island.number]][["generation"]]
final.population.list <- rcpp.res[[island.number]][["current_fitness"]]
chromosome.list <- final.population.list[["gen_original_cols"]]
chromosome.dt <- as.data.table(do.call(rbind, chromosome.list))
colnames(chromosome.dt) <- paste0("snp", seq_len(chromosome.size))
fitness.score.dt <- data.table(fitness.score =
final.population.list[["fitness_scores"]])
if (!E_GADGETS){
dif.vec.list <- final.population.list[["sum_dif_vecs"]]
dif.vec.dt <- as.data.table(do.call(rbind, dif.vec.list))
colnames(dif.vec.dt) <- paste0("snp",
seq_len(chromosome.size), ".diff.vec")
risk.allele.vec.list <- final.population.list[["risk_allele_vecs"]]
risk.allele.vec.dt <- as.data.table(do.call(rbind,
risk.allele.vec.list))
colnames(risk.allele.vec.dt) <- paste0("snp",
seq_len(chromosome.size), ".allele.copies")
n.case.risk.geno.dt <- data.table(n.cases.risk.geno =
final.population.list[["n_case_risk_geno_vec"]])
n.comp.risk.geno.dt <- data.table(n.comps.risk.geno =
final.population.list[["n_comp_risk_geno_vec"]])
final.result <- cbind(chromosome.dt, dif.vec.dt,
risk.allele.vec.dt, fitness.score.dt,
n.case.risk.geno.dt, n.comp.risk.geno.dt)
setorder(final.result, -fitness.score)
} else {
# for E-GADGETS, the dif vecs to be used for risk allele nominations
# is stored in the "risk_allele_vecs" part of the list
dif.vec.list <- final.population.list[["risk_allele_vecs"]]
dif.vec.dt <- as.data.table(do.call(rbind, dif.vec.list))
colnames(dif.vec.dt) <- paste0("snp",
seq_len(chromosome.size),
".diff.vec")
risk.allele.vec.dt <- as.data.table(matrix("1+",
nrow(dif.vec.dt),
ncol(dif.vec.dt)))
colnames(risk.allele.vec.dt) <- paste0("snp",
seq_len(chromosome.size),
".allele.copies")
# also get the betas from the lm of prob disease geno
# based on exposure
beta.list <- final.population.list[["beta_exposure_prob_disease_vecs"]]
beta.dt <- as.data.table(do.call(rbind, beta.list))
colnames(beta.dt) <- paste0("risk.exp.beta",
seq_len(ncol(beta.dt)))
final.result <- cbind(chromosome.dt, dif.vec.dt, risk.allele.vec.dt,
beta.dt, fitness.score.dt)
setorder(final.result, -fitness.score)
}
#output list
final.list <- list(top.chromosome.results = final.result,
n.generations = n.generations)
#write to file
out.file <- file.path(results.dir, paste0("cluster",
cluster.number,
".island",
island.number,".rds"))
saveRDS(final.list, out.file)
})
}
mnodzenski/epistasisGA documentation built on Jan. 17, 2023, 7:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions GADGETS
Documented in GADGETS
#' A function to run the GADGETS method
#'
#' This function runs the GADGETS method on a given cluster of islands.
#' It is a wrapper for
#' the underlying Rcpp function run_GADGETS.
#'
#' @param cluster.number An integer indicating the cluster number
#' (used for labeling the output file).
#' @param results.dir The directory to which island results will be saved.
#' @param case.genetic.data The genetic data of the disease affected children
#' from case-parent trios or disease-discordant sibling pairs. If searching for
#' maternal SNPs that are related to risk of disease in the child, some of the
#' columns in \code{case.genetic.data} may contain maternal SNP genotypes
#' (See argument \code{mother.snps} for how to indicate which SNPs columns
#' correspond to maternal genotypes). Columns are SNP allele counts, and rows
#' are individuals.This object should be of class 'matrix'. The ordering of the
#' columns must be consistent with the LD structure specified in
#' \code{ld.block.vec}. The genotypes cannot be dosages imputed with
#' uncertainty. If any data are missing for a particular family member for a
#' particular SNP, that SNP's genotype should be coded as -9 for each member of
#' the entire family: (\code{case.genetic.data} and
#' \code{father.genetic.data}/\code{mother.genetic.data},
#' or \code{case.genetic.data} and \code{complement.genetic.data}). If running
#' the experimental E-GADGETS method, this argument should be set to a 1x1
#' matrix whose only value is 0.0, and \code{case.genetic.data.n} will be used
#' to specify the case genotypes.
#' @param complement.genetic.data A genetic dataset for the controls
#' corresponding to the genotypes in \code{case.genetic.data}.For SNPs that
#' correspond to the affected child in \code{case.genetic.data}, the
#' corresponding column in \code{complement.genetic.data} should be set equal to
#' mother allele count + father allele count - case allele count. If using
#' disease-discordant siblings this argument should be the genotypes for the
#' unaffected siblings. For SNPs in \code{case.genetic.data} that represent
#' maternal genotypes (if any) the corresponding column in
#' \code{complement.genetic.data} should be the paternal genotypes for that SNP.
#' This object should be of class 'matrix'. Columns are SNP allele counts,
#' rows are families. If not specified, \code{father.genetic.data} and
#' \code{mother.genetic.data} must be specified. The genotypes cannot be dosages
#' imputed with uncertainty. If any data are missing for a particular family for
#' a particular SNP, that SNP's genotype should be coded as -9 for the entire
#' family (\code{case.genetic.data} and \code{complement.genetic.data}) for that
#' SNP. If running the experimental E-GADGETS method, this argument should be
#' set to a 1x1 matrix whose only value is 0.0, and
#' \code{complement.genetic.data.n} will be used to specify the complement
#' genotypes.
#' @param case.genetic.data.n (experimental) A matrix, to be used in the
#' experimental E-GADGETS method, containing the same data as described above
#' for \code{case.genetic.data}, but the genotypes here are stored as floating
#' point values, as opposed to integer values in \code{case.genetic.data}. If
#' not running E-GADGETS, this should be specified as a 1x1 matrix whose only
#' value is 0.0.
#' @param mother.genetic.data.n (experimental) A matrix, to be used in the
#' experimental E-GADGETS method, containing the genotypes for the mothers of
#' \code{case.genetic.data}, where the genotypes are stored as
#' floating point values, as opposed to integer values. If not running
#' E-GADGETS, this should be specified as a 1x1 matrix whose only value is 0.0.
#' @param father.genetic.data.n (experimental) A matrix, to be used in the
#' experimental E-GADGETS method, containing the genotypes for the fathers of
#' \code{case.genetic.data}, where the genotypes are stored as
#' floating point values, as opposed to integer values. If not running
#' E-GADGETS, this should be specified as a 1x1 matrix whose only value is 0.0.
#' @param exposure.mat (experimental) A matrix of the input categorical
#' and continuous exposures, if specified, to be used in the experimental
#' E-GADGETS method. If not running E-GADGETS, this should be a 1x1 matrix whose
#' only entry is 0.0.
#' @param weight.lookup.n (experimental) A vector that maps a family weight to
#' the weighted sum of the number of different SNPs and SNPs both equal to one,
#' to be used by the experimetnal E-GADGETS method. The vector should store
#' values as floating point values, not integers. If not running E-GADGETS, this
#' argument should be specified as 0.0, and will not be used in the GA. Instead,
#' for GADGETS, computation of the family weights will be based on argument
#' \code{weight.lookup}, which is computed in the same way, except stores values
#' as integers.
#' @param ld.block.vec An integer vector specifying the linkage blocks of the
#' input SNPs. As an example, for 100 candidate SNPs, suppose we specify
#' \code{ld.block.vec <- c(25, 75, 100)}. This vector indicates that the input
#' genetic data has 3 distinct linkage blocks, with SNPs 1-25 in the first
#' linkage block, 26-75 in the second block, and 76-100 in the third block. Note
#' that this means the ordering of the columns (SNPs) in
#' \code{case.genetic.data} must be consistent with the LD blocks specified in
#' \code{ld.block.vec}. In the absence of outside information, a reasonable
#' default is to consider SNPs to be in LD if they are located on the same
#' biological chromosome. If \code{case.genetic.data} includes both maternal
#' and child SNP genotypes, we recommend considering any maternal SNP and any
#' child SNP located on the same nominal biological chromosome as 'in linkage'.
#' E.g., we recommend considering any maternal SNPs located on chromosome 1
#' as being 'linked' to any child SNPs located on chromosome 1, even though,
#' strictly speaking, the maternal and child SNPs are located on separate pieces
#' of DNA. If running E-GADGETS, this argument should be specified as 0, and
#' will not be used.
#' @param n.chromosomes An integer specifying the number of chromosomes to use
#' in the GA.
#' @param chromosome.size An integer specifying the number of SNPs on each
#' chromosome.
#' @param snp.chisq A vector of statistics to be used in sampling SNPs for
#' mutation. By default, these are the square roots of
#' the chi-square marginal SNP-disease association statistics for each column in
#' \code{case.genetic.data}, but can also be manually specified or uniformly 1
#' (corresponding to totally random sampling).
#' @param weight.lookup A vector that maps a family weight to the weighted sum
#' of the number of different SNPs and SNPs both equal to one. This should
#' store values as integers.
#' @param null.mean.vec (experimental) A vector of estimated null means for each
#' component of the E-GADGETS (GxGxE) fitness score. For all other uses, this
#' should be specified as rep(0, 2) and will not be used.
#' @param null.se.vec (experimental) A vector of estimated null standard
#' deviations for each component of the E-GADGETS (GxGxE) fitness score. For all
#' other uses, this should be specified as rep(0, 2) and will not be used.
#' @param island.cluster.size An integer specifying the number of islands in
#' the cluster. See code{run.gadgets} for additional details.
#' @param n.migrations The number of chromosomes that migrate among islands.
#' This value must be less than \code{n.chromosomes} and greater than 0,
#' defaulting to 20.
#' @param n.different.snps.weight The number by which the number of different
#' SNPs between a case and complement is multiplied in computing the family
#' weights. Defaults to 2.
#' @param n.both.one.weight The number by which the number of SNPs equal to 1 in
#' both the case and complement is multiplied in computing the family weights.
#' Defaults to 1.
#' @param migration.interval The interval of generations for which GADGETS will
#' run prior to migration of top chromosomes among islands in a cluster.
#' Defaults to 50. In other words, top chromosomes will migrate among cluster
#' islands every \code{migration.interval} generations. We also check for
#' convergence at each of these intervals.
#' @param gen.same.fitness The number of consecutive generations with the same
#' fitness score required for algorithm termination. Defaults to 50.
#' @param max.generations The maximum number of generations for which GADGETS
#' will run. Defaults to 500.
#' @param initial.sample.duplicates A logical indicating whether the same SNP
#' can appear in more than one chromosome in the initial sample of chromosomes
#' (the same SNP may appear in more than one chromosome thereafter, regardless).
#' Defaults to FALSE.
#' @param crossover.prop A numeric between 0 and 1 indicating the proportion of
#' chromosomes to be subjected to cross over. The remaining proportion will be
#' mutated. Defaults to 0.8.
#' @param recessive.ref.prop The proportion to which the observed proportion of
#' informative cases with the provisional risk genotype(s) will be compared
#' to determine whether to recode the SNP as recessive. Defaults to 0.75.
#' @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. See the GADGETS paper for specific details.
#' @param E_GADGETS (experimental) A boolean indicating whether to run the
#' experimental 'E_GADGETS' method.
#' @return For each island in the cluster, an rds object containing a list with
#' the following elements will be written to \code{results.dir}.
#' \describe{
#' \item{top.chromosome.results}{A data.table of the final generation
#' chromosomes, their fitness scores, and, for GADGETS, additional information
#' pertaining to nominated risk-related genotypes. See the package vignette for
#' an example and for additional details.}
#' \item{n.generations}{The total number of generations run.}
#' }
#'
#' @examples
#'
#' set.seed(10)
#' data(case)
#' case <- as.matrix(case)
#' data(dad)
#' dad <- as.matrix(dad)
#' data(mom)
#' mom <- as.matrix(mom)
#' data.list <- preprocess.genetic.data(case[, 1:10],
#' father.genetic.data = dad[ , 1:10],
#' mother.genetic.data = mom[ , 1:10],
#' ld.block.vec = c(10))
#'
#' chisq.stats <- sqrt(data.list$chisq.stats)
#' ld.block.vec <- data.list$ld.block.vec
#' case.genetic.data <- data.list$case.genetic.data
#' complement.genetic.data <- data.list$complement.genetic.data
#'
#' #required inputs but not actually used in function below
#' case.genetic.data.n <- matrix(0.0, 1, 1)
#' mother.genetic.data.n <- matrix(0.0, 1, 1)
#' father.genetic.data.n <- matrix(0.0, 1, 1)
#' exposure.mat <- data.list$exposure.mat + 0.0
#'
#' weight.lookup <- vapply(seq_len(6), function(x) 2^x, 1)
#' dir.create('tmp')
#' GADGETS(cluster.number = 1, results.dir = 'tmp',
#' case.genetic.data = case.genetic.data,
#' complement.genetic.data = complement.genetic.data,
#' case.genetic.data.n = case.genetic.data.n,
#' mother.genetic.data.n = mother.genetic.data.n,
#' father.genetic.data.n = father.genetic.data.n,
#' exposure.mat = exposure.mat,
#' weight.lookup.n = weight.lookup + 0.0,
#' ld.block.vec = ld.block.vec,
#' n.chromosomes = 10, chromosome.size = 3, snp.chisq = chisq.stats,
#' weight.lookup = weight.lookup, n.migrations = 2,
#' migration.interval = 5,
#' gen.same.fitness = 10, max.generations = 10)
#'
#' @importFrom data.table as.data.table setorder setDT rbindlist transpose
#' @useDynLib epistasisGA
#' @export
GADGETS <- function(cluster.number, results.dir, case.genetic.data,
complement.genetic.data, case.genetic.data.n,
mother.genetic.data.n, father.genetic.data.n, exposure.mat,
weight.lookup.n, ld.block.vec, n.chromosomes,
chromosome.size, snp.chisq,
weight.lookup, null.mean.vec = c(0, 0),
null.se.vec = c(1, 1), island.cluster.size = 4,
n.migrations = 20, n.different.snps.weight = 2,
n.both.one.weight = 1, migration.interval = 50,
gen.same.fitness = 50, max.generations = 500,
initial.sample.duplicates = FALSE, crossover.prop = 0.8,
recessive.ref.prop = 0.75, recode.test.stat = 1.64,
E_GADGETS = FALSE) {
### run rcpp version of GADGETS ##
rcpp.res <- run_GADGETS(island.cluster.size, n.migrations, ld.block.vec,
n.chromosomes, chromosome.size, weight.lookup,
snp.chisq, case.genetic.data,
complement.genetic.data, case.genetic.data.n,
mother.genetic.data.n, father.genetic.data.n,
exposure.mat, weight.lookup.n, null.mean.vec,
null.se.vec, n.different.snps.weight,
n.both.one.weight, migration.interval,
gen.same.fitness, max.generations,
initial.sample.duplicates, crossover.prop,
recessive.ref.prop, recode.test.stat,
E_GADGETS)
### clean up and output results
lapply(seq_along(rcpp.res), function(island.number) {
#pick out the pieces from rcpp output
rcpp.res.length <- length(rcpp.res[[island.number]])
n.generations <- rcpp.res[[island.number]][["generation"]]
final.population.list <- rcpp.res[[island.number]][["current_fitness"]]
chromosome.list <- final.population.list[["gen_original_cols"]]
chromosome.dt <- as.data.table(do.call(rbind, chromosome.list))
colnames(chromosome.dt) <- paste0("snp", seq_len(chromosome.size))
fitness.score.dt <- data.table(fitness.score =
final.population.list[["fitness_scores"]])
if (!E_GADGETS){
dif.vec.list <- final.population.list[["sum_dif_vecs"]]
dif.vec.dt <- as.data.table(do.call(rbind, dif.vec.list))
colnames(dif.vec.dt) <- paste0("snp",
seq_len(chromosome.size), ".diff.vec")
risk.allele.vec.list <- final.population.list[["risk_allele_vecs"]]
risk.allele.vec.dt <- as.data.table(do.call(rbind,
risk.allele.vec.list))
colnames(risk.allele.vec.dt) <- paste0("snp",
seq_len(chromosome.size), ".allele.copies")
n.case.risk.geno.dt <- data.table(n.cases.risk.geno =
final.population.list[["n_case_risk_geno_vec"]])
n.comp.risk.geno.dt <- data.table(n.comps.risk.geno =
final.population.list[["n_comp_risk_geno_vec"]])
final.result <- cbind(chromosome.dt, dif.vec.dt,
risk.allele.vec.dt, fitness.score.dt,
n.case.risk.geno.dt, n.comp.risk.geno.dt)
setorder(final.result, -fitness.score)
} else {
# for E-GADGETS, the dif vecs to be used for risk allele nominations
# is stored in the "risk_allele_vecs" part of the list
dif.vec.list <- final.population.list[["risk_allele_vecs"]]
dif.vec.dt <- as.data.table(do.call(rbind, dif.vec.list))
colnames(dif.vec.dt) <- paste0("snp",
seq_len(chromosome.size),
".diff.vec")
risk.allele.vec.dt <- as.data.table(matrix("1+",
nrow(dif.vec.dt),
ncol(dif.vec.dt)))
colnames(risk.allele.vec.dt) <- paste0("snp",
seq_len(chromosome.size),
".allele.copies")
# also get the betas from the lm of prob disease geno
# based on exposure
beta.list <- final.population.list[["beta_exposure_prob_disease_vecs"]]
beta.dt <- as.data.table(do.call(rbind, beta.list))
colnames(beta.dt) <- paste0("risk.exp.beta",
seq_len(ncol(beta.dt)))
final.result <- cbind(chromosome.dt, dif.vec.dt, risk.allele.vec.dt,
beta.dt, fitness.score.dt)
setorder(final.result, -fitness.score)
}
#output list
final.list <- list(top.chromosome.results = final.result,
n.generations = n.generations)
#write to file
out.file <- file.path(results.dir, paste0("cluster",
cluster.number,
".island",
island.number,".rds"))
saveRDS(final.list, out.file)
})
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.