Commit a294a69f authored by Rauschenberger's avatar Rauschenberger
Browse files

automation

parent dd327fff
......@@ -813,6 +813,40 @@ test.single <- function(Y,X,map,i,limit=NULL,steps=NULL,rho=c(0,0.5,1)){
return(pvalue)
}
# test.trial <- function(y,x,limit=NULL,steps=NULL,rho=c(0,0.5,1)){
#
# if(is.null(limit)){limit <- 5}
# if(is.null(steps)){steps <- c(10,20,20,50)}
#
# # check input
# 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(steps)|!is.vector(steps)){
# stop("Argument \"steps\" is no numeric vector.",call.=FALSE)
# }
# if(sum(steps)<2){
# stop("Too few permutations \"sum(steps)\".",call.=FALSE)
# }
#
# # test effects
# pvalue <- rep(x=NA,times=length(rho))
# for(j in seq_along(rho)){
# tstat <- spliceQTL:::G2.multin(
# dep.data=y,indep.data=x,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[j])$Sg[-1])
# if(sum(tstat >= tstat[1]) >= limit){break}
# }
# pvalue[j] <- mean(tstat >= tstat[1],na.rm=TRUE)
# }
#
# return(pvalue)
# }
#' @export
#' @title
......@@ -848,7 +882,7 @@ test.single <- function(Y,X,map,i,limit=NULL,steps=NULL,rho=c(0,0.5,1)){
#' @examples
#' NA
#'
test.multiple <- function(Y,X,map,rho=c(0,0.5,1),spec=4){
test.multiple <- function(Y,X,map,rho=c(0,0.5,1),spec=1){
p <- nrow(map$genes)
......@@ -873,7 +907,7 @@ test.multiple <- function(Y,X,map,rho=c(0,0.5,1),spec=4){
if(max != sum(steps)){stop("Invalid combination?",call.=FALSE)}
# parallel computation
## parallel computation
#type <- ifelse(test=.Platform$OS.type=="windows",yes="PSOCK",no="FORK")
#cluster <- parallel::makeCluster(spec=spec,type=type)
#parallel::clusterSetRNGStream(cl=cluster,iseed=1)
......@@ -881,14 +915,17 @@ test.multiple <- function(Y,X,map,rho=c(0,0.5,1),spec=4){
#start <- Sys.time()
#parallel::clusterEvalQ(cl=cluster,library(spliceQTL))
#pvalue <- parallel::parLapply(cl=cluster,X=seq_len(p),fun=function(i) test.single(Y=Y,X=X,map=map,i=i,limit=limit,steps=steps,rho=rho))
#pvalue <- parallel::parLapply(cl=cluster,X=seq_len(p),fun=function(i) test.trial(y=Y[,map$exons[[i]],drop=FALSE],x=X[,seq(from=map$snps$from[i],to=map$snps$to[i],by=1),drop=FALSE],limit=limit,steps=steps,rho=rho))
#end <- Sys.time()
#parallel::stopCluster(cluster)
#rm(cluster)
## trial
#pvalue <- parallel::mclapply(X=seq_len(p),FUN=function(i) test.single(Y=Y,X=X,map=map,i=i,limit=limit,steps=steps,rho=rho))
pvalue <- lapply(X=seq_len(p),FUN=function(i) test.single(Y=Y,X=X,map=map,i=i,limit=limit,steps=steps,rho=rho))
## sequential computation
pvalue <- lapply(X=seq_len(p),FUN=function(i) spliceQTL::test.single(Y=Y,X=X,map=map,i=i,limit=limit,steps=steps,rho=rho))
# tyding up
pvalue <- do.call(what=rbind,args=pvalue)
colnames(pvalue) <- paste0("rho=",rho)
......
......@@ -127,13 +127,13 @@
<a class="sourceLine" id="cb4-17" data-line-number="17"> exons &lt;-<span class="st"> </span>list<span class="op">$</span>exons; snps &lt;-<span class="st"> </span>list<span class="op">$</span>snps; <span class="kw">rm</span>(list)</a>
<a class="sourceLine" id="cb4-18" data-line-number="18"> </a>
<a class="sourceLine" id="cb4-19" data-line-number="19"> <span class="kw">cat</span>(<span class="st">"Adjusting samples:"</span>,<span class="st">"</span><span class="ch">\n</span><span class="st">"</span>)</a>
<a class="sourceLine" id="cb4-20" data-line-number="20"> <span class="co">#exons &lt;- spliceQTL::adjust.samples(x=exons) # slow!</span></a>
<a class="sourceLine" id="cb4-20" data-line-number="20"> exons &lt;-<span class="st"> </span>spliceQTL<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/spliceQTL/topics/adjust.samples">adjust.samples</a></span>(<span class="dt">x=</span>exons) <span class="co"># slow!</span></a>
<a class="sourceLine" id="cb4-21" data-line-number="21"> exons &lt;-<span class="st"> </span><span class="kw">asinh</span>(<span class="dt">x=</span>exons)</a>
<a class="sourceLine" id="cb4-22" data-line-number="22"> </a>
<a class="sourceLine" id="cb4-23" data-line-number="23"> <span class="kw">cat</span>(<span class="st">"Adjusting covariates:"</span>,<span class="st">"</span><span class="ch">\n</span><span class="st">"</span>)</a>
<a class="sourceLine" id="cb4-24" data-line-number="24"> names &lt;-<span class="st"> </span><span class="kw">strsplit</span>(<span class="dt">x=</span><span class="kw">colnames</span>(exons),<span class="dt">split=</span><span class="st">"_"</span>) <span class="co"># exon names</span></a>
<a class="sourceLine" id="cb4-25" data-line-number="25"> length &lt;-<span class="st"> </span><span class="kw">sapply</span>(names,<span class="cf">function</span>(x) <span class="kw">as.integer</span>(x[[<span class="dv">3</span>]])<span class="op">-</span><span class="kw">as.integer</span>(x[[<span class="dv">2</span>]])) <span class="co"># exon length</span></a>
<a class="sourceLine" id="cb4-26" data-line-number="26"> <span class="co">#exons &lt;- spliceQTL::adjust.covariates(x=exons,group=gene_id,offset=length) # slow!</span></a>
<a class="sourceLine" id="cb4-26" data-line-number="26"> exons &lt;-<span class="st"> </span>spliceQTL<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/spliceQTL/topics/adjust.covariates">adjust.covariates</a></span>(<span class="dt">x=</span>exons,<span class="dt">group=</span>gene_id,<span class="dt">offset=</span>length) <span class="co"># slow!</span></a>
<a class="sourceLine" id="cb4-27" data-line-number="27"> </a>
<a class="sourceLine" id="cb4-28" data-line-number="28"> <span class="kw">cat</span>(<span class="st">"Mapping exons:"</span>,<span class="st">"</span><span class="ch">\n</span><span class="st">"</span>)</a>
<a class="sourceLine" id="cb4-29" data-line-number="29"> map &lt;-<span class="st"> </span><span class="kw">list</span>()</a>
......@@ -154,7 +154,7 @@
<a class="sourceLine" id="cb4-44" data-line-number="44"> </a>
<a class="sourceLine" id="cb4-45" data-line-number="45"> <span class="kw">cat</span>(<span class="st">"Testing:"</span>,<span class="kw">as.character</span>(<span class="kw">Sys.time</span>()),<span class="st">" -&gt; "</span>)</a>
<a class="sourceLine" id="cb4-46" data-line-number="46"> <span class="kw">rm</span>(<span class="dt">list=</span><span class="kw">setdiff</span>(<span class="kw">ls</span>(),<span class="kw">c</span>(<span class="st">"exons"</span>,<span class="st">"snps"</span>,<span class="st">"map"</span>,<span class="st">"data"</span>,<span class="st">"chr"</span>))); <span class="kw">gc</span>()</a>
<a class="sourceLine" id="cb4-47" data-line-number="47"> pvalue &lt;-<span class="st"> </span><span class="kw"><a href="../reference/test.multiple.html">test.multiple</a></span>(<span class="dt">Y=</span>exons,<span class="dt">X=</span>snps,<span class="dt">map=</span>map,<span class="dt">rho=</span><span class="kw">c</span>(<span class="dv">0</span>,<span class="fl">0.5</span>))</a>
<a class="sourceLine" id="cb4-47" data-line-number="47"> <span class="co">#pvalue &lt;- test.multiple(Y=exons,X=snps,map=map,rho=c(0,0.5))</span></a>
<a class="sourceLine" id="cb4-48" data-line-number="48"> pvalue &lt;-<span class="st"> </span>spliceQTL<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/spliceQTL/topics/test.multiple">test.multiple</a></span>(<span class="dt">Y=</span>exons,<span class="dt">X=</span>snps,<span class="dt">map=</span>map,<span class="dt">rho=</span><span class="kw">c</span>(<span class="dv">0</span>,<span class="fl">0.5</span>),<span class="dt">spec=</span><span class="dv">1</span>)</a>
<a class="sourceLine" id="cb4-49" data-line-number="49"> <span class="kw">cat</span>(<span class="kw">as.character</span>(<span class="kw">Sys.time</span>()),<span class="st">"</span><span class="ch">\n</span><span class="st">"</span>)</a>
<a class="sourceLine" id="cb4-50" data-line-number="50"> </a>
......
......@@ -120,7 +120,7 @@
</div>
<pre class="usage"><span class='fu'>test.multiple</span>(<span class='no'>Y</span>, <span class='no'>X</span>, <span class='no'>map</span>, <span class='kw'>rho</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='fl'>0</span>, <span class='fl'>0.5</span>, <span class='fl'>1</span>), <span class='kw'>spec</span> <span class='kw'>=</span> <span class='fl'>4</span>)</pre>
<pre class="usage"><span class='fu'>test.multiple</span>(<span class='no'>Y</span>, <span class='no'>X</span>, <span class='no'>map</span>, <span class='kw'>rho</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='fl'>0</span>, <span class='fl'>0.5</span>, <span class='fl'>1</span>), <span class='kw'>spec</span> <span class='kw'>=</span> <span class='fl'>1</span>)</pre>
<h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a>Arguments</h2>
<table class="ref-arguments">
......
......@@ -4,7 +4,7 @@
\alias{test.multiple}
\title{Conduct multiple tests}
\usage{
test.multiple(Y, X, map, rho = c(0, 0.5, 1), spec = 4)
test.multiple(Y, X, map, rho = c(0, 0.5, 1), spec = 1)
}
\arguments{
\item{Y}{exon expression\strong{:}
......
......@@ -59,13 +59,13 @@ for(chr in 1:22){
exons <- list$exons; snps <- list$snps; rm(list)
cat("Adjusting samples:","\n")
#exons <- spliceQTL::adjust.samples(x=exons) # slow!
exons <- spliceQTL::adjust.samples(x=exons) # slow!
exons <- asinh(x=exons)
cat("Adjusting covariates:","\n")
names <- strsplit(x=colnames(exons),split="_") # exon names
length <- sapply(names,function(x) as.integer(x[[3]])-as.integer(x[[2]])) # exon length
#exons <- spliceQTL::adjust.covariates(x=exons,group=gene_id,offset=length) # slow!
exons <- spliceQTL::adjust.covariates(x=exons,group=gene_id,offset=length) # slow!
cat("Mapping exons:","\n")
map <- list()
......@@ -86,7 +86,7 @@ for(chr in 1:22){
cat("Testing:",as.character(Sys.time())," -> ")
rm(list=setdiff(ls(),c("exons","snps","map","data","chr"))); gc()
pvalue <- test.multiple(Y=exons,X=snps,map=map,rho=c(0,0.5))
#pvalue <- test.multiple(Y=exons,X=snps,map=map,rho=c(0,0.5))
pvalue <- spliceQTL::test.multiple(Y=exons,X=snps,map=map,rho=c(0,0.5),spec=1)
cat(as.character(Sys.time()),"\n")
......
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