consistency-evaluate.R.html 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163
  1. <!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
  2. "http://www.w3.org/TR/html4/strict.dtd">
  3. <html>
  4. <head>
  5. <title></title>
  6. <meta http-equiv="content-type" content="text/html; charset=utf-8">
  7. <style type="text/css">
  8. td.linenos { background-color: #f0f0f0; padding-right: 10px; }
  9. span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
  10. pre { line-height: 125%; }
  11. body .hll { background-color: #ffffcc }
  12. body { background: #f8f8f8; }
  13. body .c { color: #408080; font-style: italic } /* Comment */
  14. body .err { border: 1px solid #FF0000 } /* Error */
  15. body .k { color: #008000; font-weight: bold } /* Keyword */
  16. body .o { color: #666666 } /* Operator */
  17. body .ch { color: #408080; font-style: italic } /* Comment.Hashbang */
  18. body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
  19. body .cp { color: #BC7A00 } /* Comment.Preproc */
  20. body .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */
  21. body .c1 { color: #408080; font-style: italic } /* Comment.Single */
  22. body .cs { color: #408080; font-style: italic } /* Comment.Special */
  23. body .gd { color: #A00000 } /* Generic.Deleted */
  24. body .ge { font-style: italic } /* Generic.Emph */
  25. body .gr { color: #FF0000 } /* Generic.Error */
  26. body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
  27. body .gi { color: #00A000 } /* Generic.Inserted */
  28. body .go { color: #888888 } /* Generic.Output */
  29. body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
  30. body .gs { font-weight: bold } /* Generic.Strong */
  31. body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
  32. body .gt { color: #0044DD } /* Generic.Traceback */
  33. body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
  34. body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
  35. body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
  36. body .kp { color: #008000 } /* Keyword.Pseudo */
  37. body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
  38. body .kt { color: #B00040 } /* Keyword.Type */
  39. body .m { color: #666666 } /* Literal.Number */
  40. body .s { color: #BA2121 } /* Literal.String */
  41. body .na { color: #7D9029 } /* Name.Attribute */
  42. body .nb { color: #008000 } /* Name.Builtin */
  43. body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
  44. body .no { color: #880000 } /* Name.Constant */
  45. body .nd { color: #AA22FF } /* Name.Decorator */
  46. body .ni { color: #999999; font-weight: bold } /* Name.Entity */
  47. body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
  48. body .nf { color: #0000FF } /* Name.Function */
  49. body .nl { color: #A0A000 } /* Name.Label */
  50. body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
  51. body .nt { color: #008000; font-weight: bold } /* Name.Tag */
  52. body .nv { color: #19177C } /* Name.Variable */
  53. body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
  54. body .w { color: #bbbbbb } /* Text.Whitespace */
  55. body .mb { color: #666666 } /* Literal.Number.Bin */
  56. body .mf { color: #666666 } /* Literal.Number.Float */
  57. body .mh { color: #666666 } /* Literal.Number.Hex */
  58. body .mi { color: #666666 } /* Literal.Number.Integer */
  59. body .mo { color: #666666 } /* Literal.Number.Oct */
  60. body .sa { color: #BA2121 } /* Literal.String.Affix */
  61. body .sb { color: #BA2121 } /* Literal.String.Backtick */
  62. body .sc { color: #BA2121 } /* Literal.String.Char */
  63. body .dl { color: #BA2121 } /* Literal.String.Delimiter */
  64. body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
  65. body .s2 { color: #BA2121 } /* Literal.String.Double */
  66. body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
  67. body .sh { color: #BA2121 } /* Literal.String.Heredoc */
  68. body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
  69. body .sx { color: #008000 } /* Literal.String.Other */
  70. body .sr { color: #BB6688 } /* Literal.String.Regex */
  71. body .s1 { color: #BA2121 } /* Literal.String.Single */
  72. body .ss { color: #19177C } /* Literal.String.Symbol */
  73. body .bp { color: #008000 } /* Name.Builtin.Pseudo */
  74. body .fm { color: #0000FF } /* Name.Function.Magic */
  75. body .vc { color: #19177C } /* Name.Variable.Class */
  76. body .vg { color: #19177C } /* Name.Variable.Global */
  77. body .vi { color: #19177C } /* Name.Variable.Instance */
  78. body .vm { color: #19177C } /* Name.Variable.Magic */
  79. body .il { color: #666666 } /* Literal.Number.Integer.Long */
  80. </style>
  81. </head>
  82. <body>
  83. <h2></h2>
  84. <div class="highlight"><pre><span></span><span class="c1">#!/usr/bin/env Rscript</span>
  85. <span class="c1"># Script to check reproducibility of fRMA training process</span>
  86. <span class="kn">library</span><span class="p">(</span>xlsx<span class="p">)</span>
  87. <span class="kn">library</span><span class="p">(</span>frma<span class="p">)</span>
  88. <span class="kn">library</span><span class="p">(</span>frmaTools<span class="p">)</span>
  89. <span class="kn">library</span><span class="p">(</span>stringr<span class="p">)</span>
  90. <span class="kn">library</span><span class="p">(</span>magrittr<span class="p">)</span>
  91. <span class="kn">library</span><span class="p">(</span>plyr<span class="p">)</span>
  92. <span class="kn">library</span><span class="p">(</span>affy<span class="p">)</span>
  93. <span class="kn">library</span><span class="p">(</span>preprocessCore<span class="p">)</span>
  94. <span class="kn">library</span><span class="p">(</span>ggplot2<span class="p">)</span>
  95. <span class="kn">library</span><span class="p">(</span>proto<span class="p">)</span>
  96. <span class="kn">library</span><span class="p">(</span>dplyr<span class="p">)</span>
  97. <span class="kp">load</span><span class="p">(</span><span class="s">&quot;consistency.rda&quot;</span><span class="p">)</span>
  98. <span class="c1">## Select a random subset of 20 arrays from each tissue (Would rather</span>
  99. <span class="c1">## use entire dataset but ENOMEM)</span>
  100. <span class="kp">set.seed</span><span class="p">(</span><span class="m">1986</span><span class="p">)</span>
  101. norm.exprs <span class="o">&lt;-</span> <span class="kp">lapply</span><span class="p">(</span><span class="kp">names</span><span class="p">(</span>vectors<span class="p">),</span> <span class="kr">function</span><span class="p">(</span>ttype<span class="p">)</span> <span class="p">{</span>
  102. stab <span class="o">&lt;-</span> sample.tables<span class="p">[[</span>ttype<span class="p">]]</span> <span class="o">%&gt;%</span> sample_n<span class="p">(</span><span class="m">20</span><span class="p">)</span>
  103. tsmsg<span class="p">(</span><span class="s">&quot;Reading 20 random arrays for &quot;</span><span class="p">,</span> ttype<span class="p">)</span>
  104. affy <span class="o">&lt;-</span> ReadAffy<span class="p">(</span>filenames<span class="o">=</span>stab<span class="o">$</span>Filename<span class="p">,</span> sampleNames<span class="o">=</span><span class="kp">rownames</span><span class="p">(</span>stab<span class="p">))</span>
  105. tsmsg<span class="p">(</span><span class="s">&quot;Normalizing with RMA for comparison&quot;</span><span class="p">)</span>
  106. eset.rma <span class="o">&lt;-</span> rma<span class="p">(</span>affy<span class="p">)</span>
  107. rma.exprs <span class="o">&lt;-</span> eset.rma <span class="o">%&gt;%</span> exprs <span class="o">%&gt;%</span> <span class="kp">as.vector</span>
  108. tsmsg<span class="p">(</span><span class="s">&quot;Normalizing with 5 trains fRMA vector sets&quot;</span><span class="p">)</span>
  109. esets.frma <span class="o">&lt;-</span> <span class="kp">lapply</span><span class="p">(</span>vectors<span class="p">[[</span>ttype<span class="p">]],</span> <span class="m">.</span> <span class="o">%&gt;%</span> frma<span class="p">(</span>affy<span class="p">,</span> input.vecs<span class="o">=</span><span class="m">.</span><span class="p">))</span>
  110. frma.exprs <span class="o">&lt;-</span> <span class="kp">sapply</span><span class="p">(</span>esets.frma<span class="p">,</span> <span class="m">.</span> <span class="o">%&gt;%</span> exprs <span class="o">%&gt;%</span> <span class="kp">as.vector</span><span class="p">)</span>
  111. <span class="kt">data.frame</span><span class="p">(</span>RMA<span class="o">=</span>rma.exprs<span class="p">,</span> fRMA<span class="o">=</span>frma.exprs<span class="p">)</span>
  112. <span class="p">})</span> <span class="o">%&gt;%</span> setNames<span class="p">(</span><span class="kp">names</span><span class="p">(</span>vectors<span class="p">))</span>
  113. <span class="c1">## Save because the above takes a while</span>
  114. <span class="kp">save.image</span><span class="p">(</span><span class="s">&quot;consistency.rda&quot;</span><span class="p">)</span>
  115. <span class="kp">dir.create</span><span class="p">(</span><span class="s">&quot;fRMA_consistency_results&quot;</span><span class="p">,</span> <span class="kc">FALSE</span><span class="p">)</span>
  116. <span class="c1">## Compute M/A data for all pairwise comparisons</span>
  117. ma.data <span class="o">&lt;-</span> <span class="kp">lapply</span><span class="p">(</span>norm.exprs<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>normexprs<span class="p">)</span> <span class="p">{</span>
  118. normexprs <span class="o">%&gt;%</span> names <span class="o">%&gt;%</span> combn<span class="p">(</span><span class="m">2</span><span class="p">,</span> simplify<span class="o">=</span><span class="kc">FALSE</span><span class="p">)</span> <span class="o">%&gt;%</span> ldply<span class="p">(</span><span class="m">.</span> <span class="o">%&gt;%</span> <span class="p">{</span>
  119. col1 <span class="o">&lt;-</span> <span class="m">.</span><span class="p">[</span><span class="m">1</span><span class="p">]</span>
  120. col2 <span class="o">&lt;-</span> <span class="m">.</span><span class="p">[</span><span class="m">2</span><span class="p">]</span>
  121. x1 <span class="o">&lt;-</span> normexprs<span class="p">[[</span>col1<span class="p">]]</span>
  122. x2 <span class="o">&lt;-</span> normexprs<span class="p">[[</span>col2<span class="p">]]</span>
  123. <span class="kt">data.frame</span><span class="p">(</span>Comparison <span class="o">=</span> str_c<span class="p">(</span>col1<span class="p">,</span><span class="s">&quot;.vs.&quot;</span><span class="p">,</span>col2<span class="p">),</span>
  124. x1<span class="o">=</span>x1<span class="p">,</span> x2<span class="o">=</span>x2<span class="p">,</span>
  125. M <span class="o">=</span> x2 <span class="o">-</span> x1<span class="p">,</span>
  126. A <span class="o">=</span> <span class="p">(</span>x1<span class="o">+</span>x2<span class="p">)</span><span class="o">/</span><span class="m">2</span><span class="p">)</span>
  127. <span class="p">})</span>
  128. <span class="p">})</span>
  129. <span class="kr">for</span> <span class="p">(</span>ttype <span class="kr">in</span> <span class="kp">names</span><span class="p">(</span>norm.exprs<span class="p">))</span> <span class="p">{</span>
  130. madata <span class="o">&lt;-</span> ma.data<span class="p">[[</span>ttype<span class="p">]]</span>
  131. pdf<span class="p">(</span>str_c<span class="p">(</span><span class="s">&quot;fRMA_consistency_results/MA Plots for &quot;</span><span class="p">,</span> ttype<span class="p">,</span> <span class="s">&quot;.pdf&quot;</span><span class="p">))</span>
  132. <span class="c1">## MA plot for every pair of normalizations</span>
  133. ddply<span class="p">(</span>madata<span class="p">,</span> <span class="m">.</span><span class="p">(</span>Comparison<span class="p">),</span> <span class="m">.</span> <span class="o">%$%</span> <span class="p">{</span>
  134. smoothScatter<span class="p">(</span>x<span class="o">=</span>A<span class="p">,</span> y<span class="o">=</span>M<span class="p">,</span> nbin<span class="o">=</span><span class="m">512</span><span class="p">,</span> main<span class="o">=</span>Comparison<span class="p">[</span><span class="m">1</span><span class="p">],</span> xlab<span class="o">=</span><span class="s">&quot;A&quot;</span><span class="p">,</span> ylab<span class="o">=</span><span class="s">&quot;M&quot;</span><span class="p">)</span>
  135. <span class="p">})</span>
  136. dev.off<span class="p">()</span>
  137. <span class="c1">## M boxplots &amp; violin plots</span>
  138. pdf<span class="p">(</span>str_c<span class="p">(</span><span class="s">&quot;fRMA_consistency_results/M Boxplots for &quot;</span><span class="p">,</span> ttype<span class="p">,</span> <span class="s">&quot;.pdf&quot;</span><span class="p">),</span>
  139. width<span class="o">=</span><span class="m">8</span><span class="p">,</span> height<span class="o">=</span><span class="m">12</span><span class="p">)</span>
  140. p <span class="o">&lt;-</span> ggplot<span class="p">(</span>madata<span class="p">)</span> <span class="o">+</span> aes<span class="p">(</span>x<span class="o">=</span>Comparison<span class="p">,</span> y<span class="o">=</span>M<span class="p">)</span> <span class="o">+</span>
  141. scale_x_discrete<span class="p">(</span>limits <span class="o">=</span> <span class="kp">rev</span><span class="p">(</span><span class="kp">levels</span><span class="p">(</span>madata<span class="o">$</span>Comparison<span class="p">)))</span> <span class="o">+</span>
  142. coord_flip<span class="p">()</span>
  143. <span class="kp">print</span><span class="p">(</span>p <span class="o">+</span> geom_boxplot<span class="p">(</span>notch<span class="o">=</span><span class="kc">TRUE</span><span class="p">,</span> outlier.shape <span class="o">=</span> <span class="kc">NA</span><span class="p">)</span> <span class="o">+</span>
  144. ggtitle<span class="p">(</span>str_c<span class="p">(</span><span class="s">&quot;Boxplots of M value distributions for &quot;</span><span class="p">,</span> ttype<span class="p">)))</span>
  145. <span class="kp">print</span><span class="p">(</span>p <span class="o">+</span> geom_violin<span class="p">(</span>scale<span class="o">=</span><span class="s">&quot;width&quot;</span><span class="p">)</span> <span class="o">+</span>
  146. ggtitle<span class="p">(</span>str_c<span class="p">(</span><span class="s">&quot;Violin plots of M value distributions for &quot;</span><span class="p">,</span> ttype<span class="p">)))</span>
  147. dev.off<span class="p">()</span>
  148. <span class="p">}</span>
  149. </pre></div>
  150. </body>
  151. </html>