automation

docs/articles/vignette.html

@@ -183,7 +183,7 @@
 <a class="sourceLine" id="cb4-67" data-line-number="67"></a>
 <a class="sourceLine" id="cb4-68" data-line-number="68">}</a></code></pre></div>
 <p>On the virtual machine, first restart R and then execute this chunk to test for alternative splicing. (It seems that lme4::lmer in spliceQTL::adjust.covariates fails to release memory on Linux.)</p>
-<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb5-1" data-line-number="1"><span class="cf">for</span>(chr <span class="cf">in</span> <span class="dv">22</span><span class="op">:</span><span class="dv">1</span>){</a>
+<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb5-1" data-line-number="1"><span class="cf">for</span>(chr <span class="cf">in</span> <span class="kw">c</span>(<span class="dv">22</span>,<span class="dv">1</span>,<span class="dv">21</span><span class="op">:</span><span class="dv">2</span>)){</a>
 <a class="sourceLine" id="cb5-2" data-line-number="2">    <span class="cf">for</span>(data <span class="cf">in</span> <span class="kw">c</span>(<span class="st">"Geuvadis"</span>,<span class="st">"LLS"</span>)){</a>
 <a class="sourceLine" id="cb5-3" data-line-number="3">        <span class="kw">cat</span>(<span class="st">"Analysing"</span>,data,chr,<span class="st">":"</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="cb5-4" data-line-number="4">        </a>
@@ -195,7 +195,7 @@
 <a class="sourceLine" id="cb5-10" data-line-number="10">        </a>
 <a class="sourceLine" id="cb5-11" data-line-number="11">        <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">"data"</span>,<span class="st">"chr"</span>,<span class="st">"path"</span>))); <span class="kw">gc</span>(); <span class="kw">cat</span>(<span class="st">"."</span>)</a>
 <a class="sourceLine" id="cb5-12" data-line-number="12">        <span class="kw">load</span>(<span class="kw">file.path</span>(path,<span class="kw">paste0</span>(<span class="st">"temp."</span>,data,<span class="st">".chr"</span>,chr,<span class="st">".RData"</span>))); <span class="kw">cat</span>(<span class="st">"."</span>)</a>
-<a class="sourceLine" id="cb5-13" data-line-number="13">        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="dv">1</span>,<span class="dt">spec=</span><span class="dv">16</span>,<span class="dt">steps=</span><span class="dv">30</span>); <span class="kw">cat</span>(<span class="st">"."</span>)</a>
+<a class="sourceLine" id="cb5-13" data-line-number="13">        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="dv">1</span>),<span class="dt">spec=</span><span class="dv">16</span>); <span class="kw">cat</span>(<span class="st">"."</span>)</a>
 <a class="sourceLine" id="cb5-14" data-line-number="14">        <span class="kw">save</span>(<span class="dt">object=</span>pvalue,<span class="dt">file=</span><span class="kw">file.path</span>(path,<span class="kw">paste0</span>(<span class="st">"pval."</span>,data,<span class="st">".chr"</span>,chr,<span class="st">".RData"</span>))); <span class="kw">cat</span>(<span class="st">"</span><span class="ch">\n</span><span class="st">"</span>)</a>
 <a class="sourceLine" id="cb5-15" data-line-number="15">    }</a>
 <a class="sourceLine" id="cb5-16" data-line-number="16"><span class="co">#q()</span></a>
@@ -205,21 +205,22 @@
 <p>On the virtual machine, execute this chunk to compare the results between the Geuvadis and the BBMRI project.</p>
 <div class="sourceCode" id="cb6"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb6-1" data-line-number="1">cor &lt;-<span class="st"> </span>chisq &lt;-<span class="st"> </span><span class="kw">rep</span>(<span class="ot">NA</span>,<span class="dt">length=</span><span class="dv">22</span>)</a>
 <a class="sourceLine" id="cb6-2" data-line-number="2"><span class="cf">for</span>(chr <span class="cf">in</span> <span class="dv">22</span><span class="op">:</span><span class="dv">1</span>){</a>
-<a class="sourceLine" id="cb6-3" data-line-number="3">    <span class="kw">load</span>(<span class="kw">file.path</span>(path,<span class="kw">paste0</span>(<span class="st">"pval.Geuvadis.chr"</span>,chr,<span class="st">".RData"</span>)))</a>
-<a class="sourceLine" id="cb6-4" data-line-number="4">    a &lt;-<span class="st"> </span>pvalue; pvalue &lt;-<span class="st"> </span><span class="ot">NA</span></a>
-<a class="sourceLine" id="cb6-5" data-line-number="5">    <span class="kw">load</span>(<span class="kw">file.path</span>(path,<span class="kw">paste0</span>(<span class="st">"pval.LLS.chr"</span>,chr,<span class="st">".RData"</span>)))</a>
-<a class="sourceLine" id="cb6-6" data-line-number="6">    b &lt;-<span class="st"> </span>pvalue; pvalue &lt;-<span class="st"> </span><span class="ot">NA</span></a>
-<a class="sourceLine" id="cb6-7" data-line-number="7"></a>
-<a class="sourceLine" id="cb6-8" data-line-number="8">    names &lt;-<span class="st"> </span><span class="kw">intersect</span>(<span class="kw">rownames</span>(a),<span class="kw">rownames</span>(b))</a>
-<a class="sourceLine" id="cb6-9" data-line-number="9">    <span class="kw">plot</span>(<span class="kw">jitter</span>(<span class="op">-</span><span class="kw">log</span>(a[names,<span class="dv">1</span>])),<span class="kw">jitter</span>(<span class="op">-</span><span class="kw">log</span>(b[names,<span class="dv">1</span>])))</a>
-<a class="sourceLine" id="cb6-10" data-line-number="10">    cor[chr] &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/cor">cor</a></span>(a[names,<span class="dv">1</span>],b[names,<span class="dv">1</span>],<span class="dt">method=</span><span class="st">"spearman"</span>)</a>
-<a class="sourceLine" id="cb6-11" data-line-number="11"></a>
-<a class="sourceLine" id="cb6-12" data-line-number="12">    a &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/p.adjust">p.adjust</a></span>(a[names,<span class="dv">1</span>])<span class="op">&lt;</span><span class="fl">0.05</span></a>
-<a class="sourceLine" id="cb6-13" data-line-number="13">    b &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/p.adjust">p.adjust</a></span>(b[names,<span class="dv">1</span>])<span class="op">&lt;</span><span class="fl">0.05</span></a>
-<a class="sourceLine" id="cb6-14" data-line-number="14">    <span class="kw">print</span>(<span class="kw">paste0</span>(<span class="st">"chr"</span>,chr))</a>
-<a class="sourceLine" id="cb6-15" data-line-number="15">    <span class="kw">print</span>(<span class="kw">table</span>(a,b))</a>
-<a class="sourceLine" id="cb6-16" data-line-number="16">    chisq[chr] &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/chisq.test">chisq.test</a></span>(<span class="kw">table</span>(a,b))<span class="op">$</span>p.value</a>
-<a class="sourceLine" id="cb6-17" data-line-number="17">}</a></code></pre></div>
+<a class="sourceLine" id="cb6-3" data-line-number="3">    sel &lt;-<span class="st"> "rho=1"</span></a>
+<a class="sourceLine" id="cb6-4" data-line-number="4">    <span class="kw">load</span>(<span class="kw">file.path</span>(path,<span class="kw">paste0</span>(<span class="st">"pval.Geuvadis.chr"</span>,chr,<span class="st">".RData"</span>)))</a>
+<a class="sourceLine" id="cb6-5" data-line-number="5">    a &lt;-<span class="st"> </span>pvalue; pvalue &lt;-<span class="st"> </span><span class="ot">NA</span></a>
+<a class="sourceLine" id="cb6-6" data-line-number="6">    <span class="kw">load</span>(<span class="kw">file.path</span>(path,<span class="kw">paste0</span>(<span class="st">"pval.LLS.chr"</span>,chr,<span class="st">".RData"</span>)))</a>
+<a class="sourceLine" id="cb6-7" data-line-number="7">    b &lt;-<span class="st"> </span>pvalue; pvalue &lt;-<span class="st"> </span><span class="ot">NA</span></a>
+<a class="sourceLine" id="cb6-8" data-line-number="8"></a>
+<a class="sourceLine" id="cb6-9" data-line-number="9">    names &lt;-<span class="st"> </span><span class="kw">intersect</span>(<span class="kw">rownames</span>(a),<span class="kw">rownames</span>(b))</a>
+<a class="sourceLine" id="cb6-10" data-line-number="10">    <span class="kw">plot</span>(<span class="kw">jitter</span>(<span class="op">-</span><span class="kw">log</span>(a[names,sel])),<span class="kw">jitter</span>(<span class="op">-</span><span class="kw">log</span>(b[names,sel])))</a>
+<a class="sourceLine" id="cb6-11" data-line-number="11">    cor[chr] &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/cor">cor</a></span>(a[names,sel],b[names,sel],<span class="dt">method=</span><span class="st">"spearman"</span>)</a>
+<a class="sourceLine" id="cb6-12" data-line-number="12"></a>
+<a class="sourceLine" id="cb6-13" data-line-number="13">    a &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/p.adjust">p.adjust</a></span>(a[names,sel])<span class="op">&lt;</span><span class="fl">0.05</span></a>
+<a class="sourceLine" id="cb6-14" data-line-number="14">    b &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/p.adjust">p.adjust</a></span>(b[names,sel])<span class="op">&lt;</span><span class="fl">0.05</span></a>
+<a class="sourceLine" id="cb6-15" data-line-number="15">    <span class="kw">print</span>(<span class="kw">paste0</span>(<span class="st">"chr"</span>,chr))</a>
+<a class="sourceLine" id="cb6-16" data-line-number="16">    <span class="kw">print</span>(<span class="kw">table</span>(a,b))</a>
+<a class="sourceLine" id="cb6-17" data-line-number="17">    chisq[chr] &lt;-<span class="st"> </span>stats<span class="op">::</span><span class="kw"><a href="http://www.rdocumentation.org/packages/stats/topics/chisq.test">chisq.test</a></span>(<span class="kw">table</span>(a,b))<span class="op">$</span>p.value</a>
+<a class="sourceLine" id="cb6-18" data-line-number="18">}</a></code></pre></div>
 <p>memory and CPU usage Linux: htop</p>
   </div>
 

vignettes/vignette.Rmd

@@ -119,7 +119,7 @@ for(chr in 1:22){
 On the virtual machine, first restart R and then execute this chunk to test for alternative splicing. (It seems that lme4::lmer in spliceQTL::adjust.covariates fails to release memory on Linux.)
 
 ```{r test,eval=FALSE}
-for(chr in 22:1){
+for(chr in c(22,1,21:2)){
     for(data in c("Geuvadis","LLS")){
         cat("Analysing",data,chr,":",as.character(Sys.time()),"\n")
         
@@ -131,7 +131,7 @@ for(chr in 22:1){
         
         rm(list=setdiff(ls(),c("data","chr","path"))); gc(); cat(".")
         load(file.path(path,paste0("temp.",data,".chr",chr,".RData"))); cat(".")
-        pvalue <- spliceQTL::test.multiple(Y=exons,X=snps,map=map,rho=1,spec=16,steps=30); cat(".")
+        pvalue <- spliceQTL::test.multiple(Y=exons,X=snps,map=map,rho=c(0,1),spec=16); cat(".")
         save(object=pvalue,file=file.path(path,paste0("pval.",data,".chr",chr,".RData"))); cat("\n")
     }
 #q()
@@ -145,17 +145,18 @@ On the virtual machine, execute this chunk to compare the results between the Ge
 ```{r,eval=FALSE}
 cor <- chisq <- rep(NA,length=22)
 for(chr in 22:1){
+    sel <- "rho=1"
     load(file.path(path,paste0("pval.Geuvadis.chr",chr,".RData")))
     a <- pvalue; pvalue <- NA
     load(file.path(path,paste0("pval.LLS.chr",chr,".RData")))
     b <- pvalue; pvalue <- NA
 
     names <- intersect(rownames(a),rownames(b))
-    plot(jitter(-log(a[names,1])),jitter(-log(b[names,1])))
-    cor[chr] <- stats::cor(a[names,1],b[names,1],method="spearman")
+    plot(jitter(-log(a[names,sel])),jitter(-log(b[names,sel])))
+    cor[chr] <- stats::cor(a[names,sel],b[names,sel],method="spearman")
 
-    a <- stats::p.adjust(a[names,1])<0.05
-    b <- stats::p.adjust(b[names,1])<0.05
+    a <- stats::p.adjust(a[names,sel])<0.05
+    b <- stats::p.adjust(b[names,sel])<0.05
     print(paste0("chr",chr))
     print(table(a,b))
     chisq[chr] <- stats::chisq.test(table(a,b))$p.value