Commit 74a4b5b9 authored by Rauschenberger's avatar Rauschenberger
Browse files

automation

parent 1fc69fad
......@@ -133,8 +133,11 @@ adjust.covariates <- function(x,offset,group){
if(!is.numeric(x)|!is.matrix(x)){
stop("Argument \"x\" is no numeric matrix.",call.=FALSE)
}
if(!is.numeric(offset)|!is.vector(offset)|any(offset<0)){
stop("Argument \"offset\" is no (positive) numeric vector.",call.=FALSE)
if(!is.numeric(offset)|!is.vector(offset)){
stop("Argument \"offset\" is no numeric vector.",call.=FALSE)
}
if(any(offset<0)){
stop("Argument \"offset\" takes negative values",call.=FALSE)
}
if(!is.character(group)|!is.vector(group)){
stop("Argument \"group\" is no character vector.",call.=FALSE)
......@@ -385,43 +388,57 @@ drop.trivial.genes <- function(map){
#'
test.single <- function(Y,X,map,i,limit,steps){
### extraction ###
if(!is.numeric(limit)){
stop("Argument \"limit\" is not numeric.",call.=FALSE)
}
if(limit<1){
stop("Argument \"limit\" is below one.",call.=FALSE)
}
if(!is.numeric(step)|!is.vector(steps)){
stop("Argument \"steps\" is no numeric vector.",call.=FALSE)
}
if(sum(steps)<2){
stop("Too few permutations \"sum(steps)\".",call.=FALSE)
}
ys <- map$exons[[i]]
y <- Y[,ys,drop=FALSE]
xs <- seq(from=map$snps$from[i],to=map$snps$to[i],by=1)
x <- X[,xs,drop=FALSE]
pvalue <- list()
### level 1: spliceQTL test ###
for(rho in c(0,0.2,0.5,1)){
rho <- c(0,0.2,0.5,1)
pvalue <- rep(x=NA,times=length(rho))
### spliceQTL test ###
for(j in seq_along(rho)){
tstat <- spliceQTL:::G2.multin(dep.data=y,indep.data=x,
nperm=steps[1]-1,rho=rho)$Sg
nperm=steps[1]-1,rho=rho[j])$Sg
for(nperm in steps[-1]){
tstat <- c(tstat,spliceQTL:::G2.multin(dep.data=y,indep.data=x,
nperm=nperm,rho=rho)$Sg[-1])
nperm=nperm,rho=rho[j])$Sg[-1])
if(sum(tstat >= tstat[1]) >= limit){break}
}
pvalue[[paste("splice",rho,sep="")]] <- mean(tstat >= tstat[1],na.rm=TRUE)
pvalue[j] <- mean(tstat >= tstat[1],na.rm=TRUE)
}
#### level 2: global test ###
pvalue$global <- rep(NA,times=length(ys))
names(pvalue$global) <- colnames(y)
for(j in seq_along(ys)){
gt <- globaltest::gt(response=y[,j],alternative=x)
pvalue$global[j] <- globaltest::p.value(gt)
}
#pvalue$global <- rep(NA,times=length(ys))
#names(pvalue$global) <- colnames(y)
#for(j in seq_along(ys)){
# gt <- globaltest::gt(response=y[,j],alternative=x)
# pvalue$global[j] <- globaltest::p.value(gt)
#}
### level 3: pairwise test ###
pvalue$pair <- matrix(NA,nrow=length(ys),ncol=length(xs))
rownames(pvalue$pair) <- colnames(y)
colnames(pvalue$pair) <- colnames(x)
for(j in seq_along(ys)){
for(k in seq_along(xs)){
lm <- stats::lm(y[,j] ~ x[,k])
pvalue$pair[j,k] <- summary(lm)$coef["x[, k]","Pr(>|t|)"]
}
}
#pvalue$pair <- matrix(NA,nrow=length(ys),ncol=length(xs))
#rownames(pvalue$pair) <- colnames(y)
#colnames(pvalue$pair) <- colnames(x)
#for(j in seq_along(ys)){
# for(k in seq_along(xs)){
# lm <- stats::lm(y[,j] ~ x[,k])
# pvalue$pair[j,k] <- summary(lm)$coef["x[,k]","Pr(>|t|)"]
# }
#}
return(pvalue)
}
......@@ -445,9 +462,9 @@ test.multiple <- function(Y,X,map){
cluster <- parallel::makeCluster(spec=8,type=type)
parallel::clusterSetRNGStream(cl=cluster,iseed=1)
test.single <- spliceQTL::test.single
parallel::clusterExport(cl=cluster,varlist=c("test.single","Y","X","map","limit","steps"),envir=environmnet())
parallel::clusterExport(cl=cluster,varlist=c("test.single","Y","X","map","limit","steps"),envir=environment())
parallel::clusterEvalQ(cl=cluster,expr=library(spliceQTL))
parallel::clusterEvalQ(cl=cluster,expr=library(globaltest))
#parallel::clusterEvalQ(cl=cluster,expr=library(globaltest))
start <- Sys.time()
pvalue <- parallel::parLapply(cl=cluster,X=1:p,fun=function(i) test.single(Y=Y,X=X,map=map,i=i,limit=limit,steps=steps))
end <- Sys.time()
......@@ -455,10 +472,11 @@ test.multiple <- function(Y,X,map){
rm(cluster)
# tyding up
names <- names(pvalue[[1]])
pvalue <- lapply(names,function(x) lapply(pvalue,function(y) y[[x]]))
names(pvalue) <- names
attributes(pvalue)$time <- end - start
#names <- names(pvalue[[1]])
#pvalue <- lapply(names,function(x) lapply(pvalue,function(y) y[[x]]))
#names(pvalue) <- names
#attributes(pvalue)$time <- end - start
return(pvalue)
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment