#' Pre-compute QR decompositions from genome cache for efficient fixed effect model mediation scans.
#'
#' This function calculates the QR decomposition for all null model and alternative models (per chromatin mediators)
#' for a genome chromatin mediation scan. This is only possible with a model that does not include any random effects.
#' Currently designed for chromatin mediation (measured by ATAC-seq), though could be generalized.
#'
#' @param genomecache The path to the genome cache directory. The genome cache is a particularly structured
#' directory that stores the haplotype probabilities/dosages at each locus. It has an additive model
#' subdirectory and a full model subdirectory. Each contains subdirectories for each chromosome, which then
#' store .RData files for the probabilities/dosages of each locus.
#' @param id DEFAULT: "SUBJECT.NAME". The is the individual-level ID that is associated with data points
#' in the phenotype data. This should be unique for each data point.
#' @param pos.is.bp DEFAULT: TRUE. If the pos variable in data is bp, use the default. If it is Mb, then
#' set to FALSE.
#' @param locus The marker or interval at which the QTL is detected. For the mediation scan, the detected
#' QTL is fit in the alternative model for each mediator.
#' @param data A data frame with mediators and potential covariates. Should also have IDs
#' that link to IDs in the genome cache, often with the individual-level ID named "SUBJECT.NAME", though others
#' can be specified with pheno.id. Currently designed to filter data to the specified covariates and the
#' chromatin variables.
#' @param formula The right side of an lm-style formula with functions of covariates contained in data frame. First
#' symbol should be "~". Functions of the outcome will be specified in mediation.scan.qr().
#' @param model DEFAULT: "additive". Specifies how to model the founder haplotype probabilities. The additive options specifies
#' use of haplotype dosages, and is most commonly used. The full option regresses the phenotype on the actual
#' diplotype probabilities.
#' @param condition.loci DEFAULT: NULL. If loci are specified, they will be included in the null and alternative models.
#' The loci must be present in the genome cache. Alternatively, conditional scans can be accomplished by regressing out
#' the loci effects with condition.out.locus.for.scan(), which does not require a new qr.object.
#' @param chr DEFAULT: "all". Specifies which chromosomes to scan.
#' @param use.progress.bar DEFAULT: TRUE. Results in a progress bar being displayed.
#' @export extract.mediation.qr
#' @examples extract.mediation.qr()
extract.mediation.qr <- function(genomecache,
id="SUBJECT.NAME",
pos.is.bp = TRUE,
locus,
data,
formula,
model=c("additive", "full"),
condition.loci=NULL,
chr="all",
use.progress.bar=TRUE){
K <- NULL
model <- model[1]
## Extracting chromatin phenotypes
chromatin <- grep(names(data), pattern="^chr", perl=TRUE, value=TRUE)
## Extracting chr for chromatin phenotypes
chromatin.chr <- gsub(chromatin, pattern="\\.[0-9]+\\.[0-9]+$", replacement="", perl=TRUE)
chromatin.chr <- gsub(chromatin.chr, pattern="^chr", replacement="", perl=TRUE)
## Extracting chromatin position
pos <- gsub(chromatin, pattern="^chr[0-9X]+\\.", replacement="", perl=TRUE)
if (pos.is.bp) {
pos <- as.numeric(gsub(pos, pattern="\\.[0-9]+\\.", replacement="", perl=TRUE))/1e6
}
h <- DiploprobReader$new(genomecache)
founders <- h$getFounders()
num.founders <- length(founders)
locus.chr <- h$getChromOfLocus(locus)
if(chr[1] != "all"){
chromatin <- chromatin[chromatin.chr %in% chr]
chromatin.chr <- chromatin.chr[chromatin.chr %in% chr]
}
subjects <- as.character(data[,id])
X.0 <- model.matrix(formula, data=data)
X.locus <- h$getLocusMatrix(locus, model=model, subjects=subjects)
keep.col <- 1:ncol(X.locus)
max.column <- which.max(colSums(X.locus, na.rm=TRUE))[1]
keep.col <- keep.col[keep.col != max.column]
if(!is.null(condition.loci)){
for(i in 1:length(condition.loci)){
X.condition <- h$getLocusMatrix(condition.loci[i], model=model, subjects=subjects)
keep.col <- 1:ncol(X.condition)
max.column <- which.max(colSums(X.condition, na.rm=TRUE))[1]
keep.col <- keep.col[keep.col != max.column]
X.0 <- cbind(X.0, X.condition[,keep.col])
}
}
if(use.progress.bar){ pb <- txtProgressBar(min=0, max=length(chromatin), style=3) }
qr.0.list <- qr.alt.list <- list()
intercept.allele <- rep(founders[max.column], length(chromatin)) # For allele effects
for(i in 1:length(chromatin)){
X <- rint(data[,chromatin[i]])
this.X.0 <- cbind(X.0, X)
this.X <- cbind(X.0, X.locus, X)
qr.0.list[[i]] <- qr(this.X.0)
qr.alt.list[[i]] <- qr(this.X)
if(use.progress.bar){ setTxtProgressBar(pb, i) }
}
names(qr.alt.list) <- names(qr.0.list) <- chromatin
qr.object <- list(qr.alt.list=qr.alt.list,
qr.0.list=qr.0.list,
X.0=X.0,
X.locus=X.locus,
intercept.allele=intercept.allele,
condition.loci=condition.loci,
chr=chromatin.chr,
locus=locus,
locus.chr=locus.chr,
pos=list(cM=NA,
Mb=pos),
model=model,
founders=founders,
subjects=subjects,
formula=Reduce(paste, deparse(formula)))
}
#' Runs genome-wide mediation scan from a QR decomposition object based on a specified QTL and phenotype data file.
#'
#' This function runs the genome scan based on a QR decomposition object and phenotype data file. This functionality only
#' works for fixed effect models. Currently written for chromatin mediation of local-eQTL, though could easily be
#' generalized to any pairing of mediator and outcome.
#'
#' @param mediation.qr.object Output object from extract.mediation.qr()
#' @param data A data frame with outcome. Should also have IDs
#' that link to IDs in the genome cache, often with the individual-level ID named "SUBJECT.NAME", though others
#' can be specified with pheno.id.
#' @param phenotype The column name (or function of column name) of a variable in data. This will become the outcome
#' of the mediation scan.
#' @param chr DEFAULT: "all". Specifies which chromosomes to scan.
#' @param id DEFAULT: "SUBJECT.NAME". The is the individual-level ID that is associated with data points
#' in the phenotype data. This should be unique for each data point.
#' @param debug.single.fit DEFAULT: FALSE. If TRUE, a browser() call is activated after the first locus is fit. This option
#' allows developers to more easily debug while still using the actual R package.
#' @param use.progress.bar DEFAULT: TRUE. Results in a progress bar being displayed.
#' @export mediation.scan.qr
#' @examples mediation.scan.qr()
mediation.scan.qr <- function(mediation.qr.object,
data,
phenotype,
chr="all",
id="SUBJECT.NAME",
debug.single.fit=FALSE,
use.progress.bar=TRUE,
...){
model <- mediation.qr.object$model
subjects <- mediation.qr.object$subjects
founders <- mediation.qr.object$founders
num.founders <- length(founders)
mediator <- names(mediation.qr.object$qr.alt.list)
mediator.chr <- mediation.qr.object$chr
mediator.pos <- mediation.qr.object$pos$Mb
rh.formula <- mediation.qr.object$formula
## Matching the subject order in the data with the qr object
reorder <- match(subjects, data[,id])
data <- data[reorder,]
n <- nrow(data)
if (all.vars(formula(paste0(phenotype, "~1")))[1] %in% mediator) {
remove.outcome <- which(all.vars(formula(paste0(phenotype, "~1")))[1] == mediator)
mediator <- mediator[-remove.outcome]
mediator.chr <- mediator.chr[-remove.outcome]
mediator.pos <- mediator.pos[-remove.outcome]
}
if (chr[1] != "all") {
mediator <- mediator[mediator.chr %in% chr]
mediator.chr <- mediator.chr[mediator.chr %in% chr]
mediator.pos <- mediator.pos[mediator.chr %in% chr]
}
formula.string <- paste(phenotype, rh.formula)
formula <- formula(formula.string)
allele.effects <- NULL
p.vec <- rep(NA, length(mediator))
y <- as.vector(model.frame(formula, data=data)[,1])
names(y) <- subjects
# Progress bar
if(use.progress.bar){ pb <- txtProgressBar(min=0, max=length(mediator), style=3) }
mediator.index <- unlist(sapply(1:length(mediator), function(x) which(mediator[x] == names(mediation.qr.object$qr.alt.list))))
for(i in 1:length(mediator)){
p.vec[i] <- get.f.stat.p.val(qr.alt=mediation.qr.object$qr.alt.list[[mediator.index[i]]],
qr.null=mediation.qr.object$qr.0.list[[mediator.index[i]]],
y=y)
if(is.nan(p.vec[i])){ p.vec[i] <- 1 }
if(debug.single.fit){ browser() }
# Update progress bar
if(use.progress.bar){ setTxtProgressBar(pb, i) }
}
names(p.vec) <- mediator
if(!is.null(mediation.qr.object$condition.loci)){
formula.string <- paste(Reduce(paste, deparse(formula)), paste(mediation.qr.object$condition.loci, collapse="+"), sep="+")
}
output <- list(LOD=NULL,
p.value=p.vec,
df=NULL,
pos=list(Mb=mediator.pos,
cM=mediation.qr.object$pos$cM),
loci=mediator,
chr=mediator.chr,
y=y,
formula=formula.string,
model.type=model,
p.value.method="ANOVA",
locus.effect.type="fixed",
locus=mediation.qr.object$locus,
n=length(y))
return(output)
}
#' Runs mediation permutation scans for significance thresholds based on an index permutation matrix object, the phenotype data, and
#' pre-computed QR decompositions for all the models.
#'
#' This function runs permutation scans in order to calculate significance thresholds for mediation analysis. Its results are most valid when
#' done in a fixed effect model setting and observations are exchangeable.
#'
#' @param perm.ind.matrix Output object from generate.qr.permutation.index.matrix().
#' @param mediation.qr.object Output object from extract.mediation.qr().
#' @param genomecache The path to the genome cache directory. The genome cache is a particularly structured
#' directory that stores the haplotype probabilities/dosages at each locus. It has an additive model
#' subdirectory and a full model subdirectory. Each contains subdirectories for each chromosome, which then
#' store .RData files for the probabilities/dosages of each locus.
#' @param phenotype The column name (or function of column name) of a variable in data. This will become the outcome
#' of the genome scan.
#' @param data A data frame with outcome and potential covariates. Should also have IDs
#' that link to IDs in the genome cache, often with the individual-level ID named "SUBJECT.NAME", though others
#' can be specified with id.
#' @param keep.full.scans DEFAULT: FALSE. If TRUE, all p-values are kept from all loci. If FALSE, only the minimum p-value
#' is kept from each scan, greatly reducing size of output.
#' @param scan.index DEFAULT: NULL. If NULL, all permutations are run. Integer vector can be specified to run just a
#' subset of the permutations.
#' @param id DEFAULT: "SUBJECT.NAME". This is the individual-level ID that is associated with data points
#' in the phenotype data. This should be unique for each data point.
#' @param chr DEFAULT: "all". Specifies which chromosomes to scan.
#' @param use.progress.bar DEFAULT: FALSE. Results in a progress bar while code runs.
#' @param pos.is.bp DEFAULT: TRUE. If the pos variable in data is bp, use the default. If it is Mb, then
#' set to FALSE.
#' @export run.qr.permutation.threshold.mediation.scans
#' @examples run.qr.permutation.threshold.mediation.scans()
run.qr.permutation.threshold.mediation.scans <- function(perm.ind.matrix,
mediation.qr.object,
genomecache,
phenotype,
data,
keep.full.scans=FALSE,
scan.index=NULL,
id="SUBJECT.NAME",
chr="all",
use.progress.bar=TRUE,
pos.is.bp = TRUE,
...){
if(is.null(scan.index)){ scan.index <- 1:ncol(perm.ind.matrix) }
model <- mediation.qr.object$model
condition.loci <- mediation.qr.object$condition.loci
locus <- mediation.qr.object$locus
chromatin <- names(mediation.qr.object$qr.0.list)
rh.formula <- mediation.qr.object$formula
chromatin.chr <- mediation.qr.object$chr
if(chr[1] != "all"){
chromatin <- chromatin[chromatin.chr %in% chr]
}
formula.string <- paste(phenotype, rh.formula, "+", id)
formula <- formula(formula.string)
if(chr == "all"){ chr.levels <- unique(chromatin.chr) }
else{ chr.levels <- chr }
full.p <- these.pos <- NULL
if(keep.full.scans){
full.p <- matrix(NA, nrow=length(scan.index), ncol=length(chromatin))
colnames(full.p) <- chromatin
these.pos <- list(Mb=mediation.qr.object$pos$Mb[chromatin.chr %in% chr.levels],
cM=mediation.qr.object$pos$cM[chromatin.chr %in% chr.levels])
}
min.p <- max.p <- rep(NA, length(scan.index))
y <- model.frame(formula, data=data)[,1]
data <- data[,!grepl(pattern="^gene_", perl=TRUE, x=colnames(data))]
permute.var <- !(colnames(data) %in% all.vars(formula)[-1])
for(i in 1:length(scan.index)){
## Permuting all variables but covariates
this.data <- data.frame(y=y, data[perm.ind.matrix[,scan.index[i]], permute.var], data[, !permute.var])
this.mediation.qr.object <- extract.mediation.qr(genomecache=genomecache, id=id,
data=this.data, formula=as.formula(rh.formula),
model=model, condition.loci=condition.loci,
chr=chr, locus=locus, use.progress.bar=FALSE,
pos.is.bp = pos.is.bp)
this.scan <- mediation.scan.qr(mediation.qr.object=this.mediation.qr.object, data=this.data,
phenotype="y", id=id, chr=chr,
return.allele.effects=FALSE, use.progress.bar=use.progress.bar,
...)
if(keep.full.scans){
full.p[i,] <- this.scan$p.value
}
min.p[i] <- min(this.scan$p.value)
max.p[i] <- max(this.scan$p.value)
cat("\n", "Threshold scan: index", scan.index[i], "complete ---------- final index of this run:", scan.index[length(scan.index)], "\n")
}
return(list(full.results=list(LOD=NULL,
p.value=full.p,
chr=chromatin.chr[chromatin.chr %in% chr.levels],
pos=these.pos),
max.statistics=list(LOD=NULL,
p.value=list(min=min.p,
max=max.p))))
}
#' Extracts chromosome-level statistics from the full output of genome-wide permutation scans, allowing chromosome-wide and genome-wide significance
#' to be determined from the same output.
#'
#' This function extracts the associations for a specified chromosome from genome-wide permutation scans. This allows both genome-wide and chromosome-wide
#' thresholds to be determined from the same output.
#'
#' @param full.perm.scans Output object from generate.qr.permutation.index.matrix().
#' @param chr The chromosome that will have its association scores extracted from the full.perm.scans, essentially
#' creating permutation scans specific to the chromosome.
#' @param use.lod DEFAULT: FALSE. If TRUE, LOD scores are recorded. If FALSE, p-values are recorded, either based on the likelihood ratio test or ANOVA.
#' @export extract.chr.max.statistics.from.genomewide
#' @examples extract.chr.max.statistics.from.genomewide()
extract.chr.max.statistics.from.genomewide <- function(full.perm.scans,
chr,
use.lod=FALSE){
statistic.type <- ifelse(use.lod, "LOD", "p.value")
full.results <- full.perm.scans$full.results[[statistic.type]][,full.perm.scans$full.results$chr == chr]
return(list(full.results=NULL,
max.statistics=list(LOD=NULL,
p.value=list(min=apply(full.results, 1, function(x) min(x, na.rm=TRUE)),
max=apply(full.results, 1, function(x) max(x, na.rm=TRUE))))))
}
#' Pre-compute QR decompositions from genome cache for efficient fixed effect model mediation scans.
#'
#' This function calculates the QR decomposition for all null model and alternative models (per gene expression mediators)
#' for a genome expression mediation scan. This is only possible with a model that does not include any random effects.
#' Currently designed for gene expression mediation (measured by RNA-seq), though could be generalized, essentially
#' consolidating extract.mediation.expression.qr() and extract.mediation.qr().
#'
#' @param genomecache The path to the genome cache directory. The genome cache is a particularly structured
#' directory that stores the haplotype probabilities/dosages at each locus. It has an additive model
#' subdirectory and a full model subdirectory. Each contains subdirectories for each chromosome, which then
#' store .RData files for the probabilities/dosages of each locus.
#' @param id DEFAULT: "SUBJECT.NAME". The is the individual-level ID that is associated with data points
#' in the phenotype data. This should be unique for each data point.
#' @param pos.is.bp DEFAULT: TRUE. If the pos variable in gene.data is bp, use the default. If it is Mb, then
#' set to FALSE.
#' @param locus The marker or interval at which the QTL is detected. For the mediation scan, the detected
#' QTL is fit in the alternative model for each mediator.
#' @param data A data frame with mediators and potential covariates. Should also have IDs
#' that link to IDs in the genome cache, often with the individual-level ID named "SUBJECT.NAME", though others
#' can be specified with pheno.id. Currently designed to filter data to the specified covariates and the
#' gene expression variables.
#' @param gene.data Gene annotation data for the gene expression variables in data. Specifically, the gene transcription
#' start site is needed, named gene.start in gene.data.
#' @param formula The right side of an lm-style formula with functions of covariates contained in data frame. First
#' symbol should be "~". Functions of the outcome will be specified in mediation.scan.qr().
#' @param model DEFAULT: "additive". Specifies how to model the founder haplotype probabilities. The additive options specifies
#' use of haplotype dosages, and is most commonly used. The full option regresses the phenotype on the actual
#' diplotype probabilities.
#' @param condition.loci DEFAULT: NULL. If loci are specified, they will be included in the null and alternative models.
#' The loci must be present in the genome cache. Alternatively, conditional scans can be accomplished by regressing out
#' the loci effects with condition.out.locus.for.scan(), which does not require a new qr.object.
#' @param chr DEFAULT: "all". Specifies which chromosomes to scan.
#' @param use.progress.bar DEFAULT: TRUE. Results in a progress bar being displayed.
#' @export extract.mediation.expression.qr
#' @examples extract.mediation.expression.qr()
extract.mediation.expression.qr <- function(genomecache,
id="SUBJECT.NAME",
pos.is.bp = TRUE,
locus,
data,
gene.data,
formula,
model=c("additive", "full"),
condition.loci=NULL,
chr="all",
use.progress.bar=TRUE){
K <- NULL
model <- model[1]
## Extracting gene expressions
genes <- grep(names(data), pattern="^gene", perl=TRUE, value=TRUE)
## Extracting chr for expression phenotypes
ref.genes <- paste0("gene_", gene.data[,4])
ref.genes <- chartr(old="()-", new="...", ref.genes)
gene.index <- unlist(sapply(1:length(genes), function(x) which(genes[x] == ref.genes)))
gene.chr <- gene.data[,1][gene.index]
gene.chr <- gsub(gene.chr, pattern="^chr", replacement="", perl=TRUE)
## Extracting gene position
pos <- gene.data[,2][gene.index]
if (pos.is.bp) {
pos <- pos/1e6
}
h <- DiploprobReader$new(genomecache)
founders <- h$getFounders()
num.founders <- length(founders)
locus.chr <- h$getChromOfLocus(locus)
if(chr[1] != "all"){
genes <- genes[gene.chr %in% chr]
gene.chr <- gene.chr[gene.chr %in% chr]
pos <- pos[gene.chr %in% chr]
}
subjects <- as.character(data[,id])
X.0 <- model.matrix(formula, data=data)
X.locus <- h$getLocusMatrix(locus, model=model, subjects=subjects)
keep.col <- 1:ncol(X.locus)
max.column <- which.max(colSums(X.locus, na.rm=TRUE))[1]
keep.col <- keep.col[keep.col != max.column]
if(!is.null(condition.loci)){
for(i in 1:length(condition.loci)){
X.condition <- h$getLocusMatrix(condition.loci[i], model=model, subjects=subjects)
keep.col <- 1:ncol(X.condition)
max.column <- which.max(colSums(X.condition, na.rm=TRUE))[1]
keep.col <- keep.col[keep.col != max.column]
X.0 <- cbind(X.0, X.condition[,keep.col])
}
}
if(use.progress.bar){ pb <- txtProgressBar(min=0, max=length(genes), style=3) }
qr.0.list <- qr.alt.list <- list()
intercept.allele <- rep(founders[max.column], length(genes)) # For allele effects
for(i in 1:length(genes)){
X <- rint(data[,genes[i]])
this.X.0 <- cbind(X.0, X)
this.X <- cbind(X.0, X.locus, X)
qr.0.list[[i]] <- qr(this.X.0)
qr.alt.list[[i]] <- qr(this.X)
if(use.progress.bar){ setTxtProgressBar(pb, i) }
}
names(qr.alt.list) <- names(qr.0.list) <- genes
qr.object <- list(qr.alt.list=qr.alt.list,
qr.0.list=qr.0.list,
X.0=X.0,
X.locus=X.locus,
intercept.allele=intercept.allele,
condition.loci=condition.loci,
chr=gene.chr,
locus=locus,
locus.chr=locus.chr,
pos=list(cM=NA,
Mb=pos),
model=model,
founders=founders,
subjects=subjects,
formula=Reduce(paste, deparse(formula)))
}
#' Runs mediation permutation scans for significance thresholds based on an index permutation matrix object, the phenotype data, and
#' pre-computed QR decompositions for all the models. Currently designed for mediation through gene expression.
#'
#' This function runs permutation scans in order to calculate significance thresholds for mediation analysis. Its results are most valid when
#' done in a fixed effect model setting and observations are exchangeable.
#'
#' @param perm.ind.matrix Output object from generate.qr.permutation.index.matrix().
#' @param mediation.qr.object Output object from extract.mediation.expression.qr().
#' @param genomecache The path to the genome cache directory. The genome cache is a particularly structured
#' directory that stores the haplotype probabilities/dosages at each locus. It has an additive model
#' subdirectory and a full model subdirectory. Each contains subdirectories for each chromosome, which then
#' store .RData files for the probabilities/dosages of each locus.
#' @param phenotype The column name (or function of column name) of a variable in data. This will become the outcome
#' of the genome scan.
#' @param data A data frame with outcome and potential covariates. Should also have IDs
#' that link to IDs in the genome cache, often with the individual-level ID named "SUBJECT.NAME", though others
#' can be specified with id.
#' @param keep.full.scans DEFAULT: FALSE. If TRUE, all p-values are kept from all loci. If FALSE, only the minimum p-value
#' is kept from each scan, greatly reducing size of output.
#' @param scan.index DEFAULT: NULL. If NULL, all permutations are run. Integer vector can be specified to run just a
#' subset of the permutations.
#' @param id DEFAULT: "SUBJECT.NAME". This is the individual-level ID that is associated with data points
#' in the phenotype data. This should be unique for each data point.
#' @param chr DEFAULT: "all". Specifies which chromosomes to scan.
#' @param use.progress.bar DEFAULT: FALSE. Results in a progress bar while code runs.
#' @param pos.is.bp DEFAULT: TRUE. If the pos variable in data is bp, use the default. If it is Mb, then
#' set to FALSE.
#' @export run.qr.permutation.threshold.mediation.expression.scans
#' @examples run.qr.permutation.threshold.mediation.expression.scans()
run.qr.permutation.threshold.mediation.expression.scans <- function(perm.ind.matrix,
mediation.qr.object,
genomecache,
phenotype,
data,
gene.data,
keep.full.scans=FALSE,
scan.index=NULL,
id="SUBJECT.NAME",
chr="all",
pos.is.bp=TRUE,
use.progress.bar=TRUE,
...){
if(is.null(scan.index)){ scan.index <- 1:ncol(perm.ind.matrix) }
model <- mediation.qr.object$model
condition.loci <- mediation.qr.object$condition.loci
locus <- mediation.qr.object$locus
genes <- names(mediation.qr.object$qr.0.list)
rh.formula <- mediation.qr.object$formula
genes.chr <- mediation.qr.object$chr
if(chr[1] != "all"){
genes <- genes[genes.chr %in% chr]
}
formula.string <- paste(phenotype, rh.formula, "+", id)
formula <- formula(formula.string)
if(chr == "all"){ chr.levels <- unique(genes.chr) }
else{ chr.levels <- chr }
full.p <- these.pos <- NULL
if(keep.full.scans){
full.p <- matrix(NA, nrow=length(scan.index), ncol=length(genes))
colnames(full.p) <- genes
these.pos <- list(Mb=mediation.qr.object$pos$Mb[genes.chr %in% chr.levels],
cM=mediation.qr.object$pos$cM[genes.chr %in% chr.levels])
}
min.p <- max.p <- rep(NA, length(scan.index))
y <- model.frame(formula, data=data)[,1]
permute.var <- !(colnames(data) %in% all.vars(formula)[-1])
for(i in 1:length(scan.index)){
## Permuting all variables but covariates
this.data <- data.frame(y=y, data[perm.ind.matrix[,scan.index[i]], permute.var], data[, !permute.var])
this.mediation.qr.object <- extract.mediation.expression.qr(genomecache=genomecache, id=id,
data=this.data,
gene.data=gene.data,
formula=as.formula(rh.formula),
pos.is.bp=pos.is.bp,
model=model, condition.loci=condition.loci,
chr=chr, locus=locus, use.progress.bar=FALSE)
this.scan <- mediation.scan.qr(mediation.qr.object=this.mediation.qr.object, data=this.data,
phenotype="y", id=id, chr=chr,
return.allele.effects=FALSE, use.progress.bar=use.progress.bar,
...)
if(keep.full.scans){
full.p[i,] <- this.scan$p.value
}
min.p[i] <- min(this.scan$p.value)
max.p[i] <- max(this.scan$p.value)
cat("\n", "Threshold scan: index", scan.index[i], "complete ---------- final index of this run:", scan.index[length(scan.index)], "\n")
}
return(list(full.results=list(LOD=NULL,
p.value=full.p,
chr=genes.chr[genes.chr %in% chr.levels],
pos=these.pos),
max.statistics=list(LOD=NULL,
p.value=list(min=min.p,
max=max.p))))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.