automation

R/functions.R

@@ -908,24 +908,26 @@ test.multiple <- function(Y,X,map,rho=c(0,0.5,1),spec=1){
     
     if(max != sum(steps)){stop("Invalid combination?",call.=FALSE)}
     
-    ## 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)
-    ##parallel::clusterExport(cl=cluster,varlist=c("Y","X","map","limit","steps","rho"),envir=environment())
-    #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)
+    if(spec==1){
+        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))
+    } else {
+        
+        ## 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)
+        parallel::clusterExport(cl=cluster,varlist=c("Y","X","map","limit","steps","rho"),envir=environment())
+        parallel::clusterEvalQ(cl=cluster,library(spliceQTL,lib.loc="/virdir/Scratch/arauschenberger/library"))
+        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))
+        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))
     
-    ## 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)

vignettes/vignette.Rmd

@@ -52,8 +52,6 @@ for(chr in 1:22){
       cond <- sapply(strsplit(x=rownames(exons),split=":"),function(x) x[[1]]==data)
       exons <- exons[cond,]
     }
-    cond <- sapply(strsplit(x=colnames(exons),split="_"),function(x) x[[1]]==chr)
-    exons <- exons[,cond]
     load(file.path(path,paste0(data,".chr",chr,".RData")),verbose=TRUE)
     
     cat("Matching samples:","\n")
@@ -61,13 +59,18 @@ 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!
+    
+    # subset chromosome
+    cond <- sapply(strsplit(x=colnames(exons),split="_"),function(x) x[[1]]==chr)
+    exons <- exons[,cond]
+    gene_id = gene_id[cond] # correct?
     
     cat("Mapping exons:","\n")
     map <- list()
@@ -116,4 +119,10 @@ for(chr in 1:22){
 ```
 
 
-memory usage Linux: ps hax -o rss,user | awk '{a[$2]+=$1;}END{for(i in a)print i" "int(a[i]/1024+0.5);}' | sort -rnk2
+memory usage Linux:
+
+ps hax -o rss,user | awk '{a[$2]+=$1;}END{for(i in a)print i" "int(a[i]/1024+0.5);}' | sort -rnk2
+
+ps -U root --no-headers -o rss | awk '{ sum+=$1} END {print int(sum/1024) "MB"}'
+
+