ryan
/
resume
spogulis no https://github.com/DarwinAwardWinner/resume.git


			
				
					
						
						
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351
							<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=utf-8">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #408080; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #008000; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .ch { color: #408080; font-style: italic } /* Comment.Hashbang */
body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
body .cp { color: #BC7A00 } /* Comment.Preproc */
body .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */
body .c1 { color: #408080; font-style: italic } /* Comment.Single */
body .cs { color: #408080; font-style: italic } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #008000 } /* Keyword.Pseudo */
body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #B00040 } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BA2121 } /* Literal.String */
body .na { color: #7D9029 } /* Name.Attribute */
body .nb { color: #008000 } /* Name.Builtin */
body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #0000FF } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #19177C } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mb { color: #666666 } /* Literal.Number.Bin */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sa { color: #BA2121 } /* Literal.String.Affix */
body .sb { color: #BA2121 } /* Literal.String.Backtick */
body .sc { color: #BA2121 } /* Literal.String.Char */
body .dl { color: #BA2121 } /* Literal.String.Delimiter */
body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BA2121 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BA2121 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BA2121 } /* Literal.String.Single */
body .ss { color: #19177C } /* Literal.String.Symbol */
body .bp { color: #008000 } /* Name.Builtin.Pseudo */
body .fm { color: #0000FF } /* Name.Function.Magic */
body .vc { color: #19177C } /* Name.Variable.Class */
body .vg { color: #19177C } /* Name.Variable.Global */
body .vi { color: #19177C } /* Name.Variable.Instance */
body .vm { color: #19177C } /* Name.Variable.Magic */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span></span><span class="c1">#!/usr/bin/env Rscript</span>

<span class="kn">library</span><span class="p">(</span>xlsx<span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>frmaTools<span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>stringr<span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>magrittr<span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>plyr<span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>affy<span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>preprocessCore<span class="p">)</span>

training.data.dir <span class="o">&lt;-</span> <span class="s">&quot;Training Data&quot;</span>
datasets <span class="o">&lt;-</span> <span class="kt">data.frame</span><span class="p">(</span>Dataset<span class="o">=</span><span class="kp">list.files</span><span class="p">(</span>training.data.dir<span class="p">))</span>
<span class="kp">rownames</span><span class="p">(</span>datasets<span class="p">)</span> <span class="o">&lt;-</span> datasets<span class="o">$</span>Dataset
datasets<span class="o">$</span>Tissue <span class="o">&lt;-</span> <span class="kp">factor</span><span class="p">(</span>str_extract<span class="p">(</span>datasets<span class="o">$</span>Dataset<span class="p">,</span> <span class="s">&quot;\\b(PAX|BX)\\b&quot;</span><span class="p">))</span>

tsmsg <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span><span class="kc">...</span><span class="p">)</span> <span class="p">{</span>
  <span class="kp">message</span><span class="p">(</span><span class="kp">date</span><span class="p">(),</span> <span class="s">&quot;: &quot;</span><span class="p">,</span> <span class="kc">...</span><span class="p">)</span>
<span class="p">}</span>

<span class="c1">## Some Scan Dates are marked as identical for multiple batches, which</span>
<span class="c1">## is bad. But the dates embedded in the file names for these batches</span>
<span class="c1">## are different, so we use those dates instead.</span>
parse.date.from.filename <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>fname<span class="p">)</span> <span class="p">{</span>
    res1 <span class="o">&lt;-</span> str_match<span class="p">(</span>fname<span class="p">,</span> <span class="s">&quot;^(\\d\\d)(\\d\\d)(\\d\\d)&quot;</span><span class="p">)[,</span><span class="kt">c</span><span class="p">(</span><span class="m">4</span><span class="p">,</span><span class="m">2</span><span class="p">,</span><span class="m">3</span><span class="p">)]</span>
    res2 <span class="o">&lt;-</span> str_match<span class="p">(</span>fname<span class="p">,</span> <span class="s">&quot;^20(\\d\\d)_(\\d\\d)_(\\d\\d)&quot;</span><span class="p">)[,</span><span class="m">-1</span><span class="p">]</span>
    res1<span class="p">[</span><span class="kp">is.na</span><span class="p">(</span>res1<span class="p">)]</span> <span class="o">&lt;-</span> res2<span class="p">[</span><span class="kp">is.na</span><span class="p">(</span>res1<span class="p">)]</span>
    <span class="kp">colnames</span><span class="p">(</span>res1<span class="p">)</span> <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span><span class="s">&quot;year&quot;</span><span class="p">,</span> <span class="s">&quot;month&quot;</span><span class="p">,</span> <span class="s">&quot;day&quot;</span><span class="p">)</span>
    res1<span class="p">[,</span><span class="s">&quot;year&quot;</span><span class="p">]</span> <span class="o">%&lt;&gt;%</span> str_c<span class="p">(</span><span class="s">&quot;20&quot;</span><span class="p">,</span> <span class="m">.</span><span class="p">)</span>
    <span class="kp">as.Date</span><span class="p">(</span><span class="kp">do.call</span><span class="p">(</span><span class="kp">ISOdate</span><span class="p">,</span> <span class="kt">data.frame</span><span class="p">(</span>res1<span class="p">)))</span>
<span class="p">}</span>

makeVectorsAffyBatch <span class="o">&lt;-</span> <span class="kr">function</span> <span class="p">(</span>files<span class="p">,</span> batch.id<span class="p">,</span> background <span class="o">=</span> <span class="s">&quot;rma&quot;</span><span class="p">,</span> normalize <span class="o">=</span> <span class="s">&quot;quantile&quot;</span><span class="p">,</span>
    normVec <span class="o">=</span> <span class="kc">NULL</span><span class="p">,</span> cdfname <span class="o">=</span> <span class="kc">NULL</span><span class="p">,</span> file.dir <span class="o">=</span> <span class="s">&quot;.&quot;</span><span class="p">,</span> verbose <span class="o">=</span> <span class="kc">TRUE</span><span class="p">)</span>
<span class="p">{</span>
    wd <span class="o">&lt;-</span> <span class="kp">getwd</span><span class="p">()</span>
    <span class="kp">setwd</span><span class="p">(</span>file.dir<span class="p">)</span>
    object <span class="o">&lt;-</span> ReadAffy<span class="p">(</span>filenames <span class="o">=</span> files<span class="p">,</span> cdfname <span class="o">=</span> cdfname<span class="p">,</span>
        verbose <span class="o">=</span> verbose<span class="p">)</span>
    <span class="kp">setwd</span><span class="p">(</span>wd<span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Data loaded \n&quot;</span><span class="p">)</span>
    batch.size <span class="o">&lt;-</span> <span class="kp">table</span><span class="p">(</span>batch.id<span class="p">)[</span><span class="m">1</span><span class="p">]</span>
    <span class="kr">if</span> <span class="p">(</span><span class="o">!</span><span class="kp">all</span><span class="p">(</span><span class="kp">table</span><span class="p">(</span>batch.id<span class="p">)</span> <span class="o">==</span> batch.size<span class="p">))</span>
        <span class="kp">stop</span><span class="p">(</span><span class="s">&quot;Batches must be of the same size.&quot;</span><span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>background <span class="o">==</span> <span class="s">&quot;rma&quot;</span><span class="p">)</span> <span class="p">{</span>
        object <span class="o">&lt;-</span> bg.correct.rma<span class="p">(</span>object<span class="p">)</span>
        <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
            <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Background Corrected \n&quot;</span><span class="p">)</span>
        <span class="kp">gc</span><span class="p">()</span>
    <span class="p">}</span>
    pms <span class="o">&lt;-</span> pm<span class="p">(</span>object<span class="p">)</span>
    pns <span class="o">&lt;-</span> probeNames<span class="p">(</span>object<span class="p">)</span>
    pmi <span class="o">&lt;-</span> <span class="kp">unlist</span><span class="p">(</span>pmindex<span class="p">(</span>object<span class="p">))</span>
    <span class="kr">if</span> <span class="p">(</span><span class="o">!</span><span class="kp">all</span><span class="p">(</span><span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;%i&quot;</span><span class="p">,</span> pmi<span class="p">)</span> <span class="o">==</span> <span class="kp">rownames</span><span class="p">(</span>pms<span class="p">)))</span>
        <span class="kp">stop</span><span class="p">(</span><span class="s">&quot;Mismatch between pmindex and rownames of pms&quot;</span><span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>object<span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    <span class="kr">if</span> <span class="p">(</span>normalize <span class="o">==</span> <span class="s">&quot;quantile&quot;</span><span class="p">)</span> <span class="p">{</span>
        <span class="kr">if</span> <span class="p">(</span><span class="kp">is.null</span><span class="p">(</span>normVec<span class="p">))</span>
            normVec <span class="o">&lt;-</span> normalize.quantiles.determine.target<span class="p">(</span>pms<span class="p">)</span>
        pms <span class="o">&lt;-</span> normalize.quantiles.use.target<span class="p">(</span>pms<span class="p">,</span> normVec<span class="p">)</span>
        <span class="kp">names</span><span class="p">(</span>normVec<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">as.character</span><span class="p">(</span>pmi<span class="p">)</span>
        <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
            <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Normalized \n&quot;</span><span class="p">)</span>
    <span class="p">}</span>
    pms <span class="o">&lt;-</span> <span class="kp">log2</span><span class="p">(</span>pms<span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    N <span class="o">&lt;-</span> <span class="m">1</span><span class="o">:</span><span class="kp">dim</span><span class="p">(</span>pms<span class="p">)[</span><span class="m">1</span><span class="p">]</span>
    S <span class="o">&lt;-</span> <span class="kp">split</span><span class="p">(</span>N<span class="p">,</span> pns<span class="p">)</span>
    nc <span class="o">&lt;-</span> <span class="kp">ncol</span><span class="p">(</span>pms<span class="p">)</span>
    nr <span class="o">&lt;-</span> <span class="kp">nrow</span><span class="p">(</span>pms<span class="p">)</span>
    resids <span class="o">&lt;-</span> <span class="kt">matrix</span><span class="p">(</span>ncol <span class="o">=</span> nc<span class="p">,</span> nrow <span class="o">=</span> nr<span class="p">)</span>
    probeVec <span class="o">&lt;-</span> <span class="kt">vector</span><span class="p">(</span>length <span class="o">=</span> nr<span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Beginning Probe Effect Calculation ... \n&quot;</span><span class="p">)</span>
    <span class="kr">for</span> <span class="p">(</span>k <span class="kr">in</span> <span class="m">1</span><span class="o">:</span><span class="kp">length</span><span class="p">(</span>S<span class="p">))</span> <span class="p">{</span>
        fit <span class="o">&lt;-</span> rcModelPLM<span class="p">(</span>pms<span class="p">[</span>S<span class="p">[[</span>k<span class="p">]],</span> <span class="p">,</span> drop <span class="o">=</span> <span class="kc">FALSE</span><span class="p">])</span>
        resids<span class="p">[</span>S<span class="p">[[</span>k<span class="p">]],</span> <span class="p">]</span> <span class="o">&lt;-</span> fit<span class="o">$</span>Residuals
        probeVec<span class="p">[</span>S<span class="p">[[</span>k<span class="p">]]]</span> <span class="o">&lt;-</span> fit<span class="o">$</span>Estimates<span class="p">[(</span>nc <span class="o">+</span> <span class="m">1</span><span class="p">)</span><span class="o">:</span><span class="kp">length</span><span class="p">(</span>fit<span class="o">$</span>Estimates<span class="p">)]</span>
        <span class="kr">if</span> <span class="p">((</span>k<span class="o">%%</span><span class="m">1000</span><span class="p">)</span> <span class="o">==</span> <span class="m">0</span><span class="p">)</span> <span class="p">{</span>
            <span class="kp">message</span><span class="p">(</span><span class="kp">paste</span><span class="p">(</span><span class="s">&quot;Finished probeset:&quot;</span><span class="p">,</span> k<span class="p">,</span> <span class="s">&quot;\n&quot;</span><span class="p">))</span>
            <span class="kp">gc</span><span class="p">()</span>
        <span class="p">}</span>
    <span class="p">}</span>
    <span class="kp">names</span><span class="p">(</span>probeVec<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">as.character</span><span class="p">(</span>pmi<span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Probe Effects Calculated \n&quot;</span><span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    tmp <span class="o">&lt;-</span> <span class="kp">split</span><span class="p">(</span><span class="kp">t</span><span class="p">(</span>resids<span class="p">),</span> batch.id<span class="p">)</span>
    withinMean <span class="o">&lt;-</span> <span class="kp">lapply</span><span class="p">(</span>tmp<span class="p">,</span> frmaTools<span class="o">:::</span>getProbeMean<span class="p">,</span> batch.size<span class="p">)</span>
    withinVar <span class="o">&lt;-</span> <span class="kp">lapply</span><span class="p">(</span>tmp<span class="p">,</span> frmaTools<span class="o">:::</span>getProbeVar<span class="p">,</span> batch.size<span class="p">)</span>
    withinAvgVar <span class="o">&lt;-</span> <span class="kp">rowMeans</span><span class="p">(</span><span class="kt">matrix</span><span class="p">(</span><span class="kp">unlist</span><span class="p">(</span>withinVar<span class="p">),</span> ncol <span class="o">=</span> <span class="kp">length</span><span class="p">(</span>withinVar<span class="p">)))</span>
    btwVar <span class="o">&lt;-</span> <span class="kp">apply</span><span class="p">(</span><span class="kt">matrix</span><span class="p">(</span><span class="kp">unlist</span><span class="p">(</span>withinMean<span class="p">),</span> ncol <span class="o">=</span> <span class="kp">length</span><span class="p">(</span>withinMean<span class="p">)),</span>
        <span class="m">1</span><span class="p">,</span> var<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>tmp<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>withinMean<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>withinVar<span class="p">)</span>
    <span class="kp">names</span><span class="p">(</span>withinAvgVar<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">names</span><span class="p">(</span>btwVar<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">as.character</span><span class="p">(</span>pmi<span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Probe Variances Calculated \n&quot;</span><span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    tmp <span class="o">&lt;-</span> <span class="kp">split</span><span class="p">(</span>resids<span class="p">,</span> pns<span class="p">)</span>
    psetMAD <span class="o">&lt;-</span> <span class="kp">unlist</span><span class="p">(</span><span class="kp">lapply</span><span class="p">(</span>tmp<span class="p">,</span> frmaTools<span class="o">:::</span>getPsetMAD<span class="p">,</span> nc<span class="p">,</span> batch.id<span class="p">))</span>
    <span class="kp">names</span><span class="p">(</span>psetMAD<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">names</span><span class="p">(</span>tmp<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>tmp<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>resids<span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Probe Set SDs Calculated \n&quot;</span><span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    w <span class="o">&lt;-</span> <span class="m">1</span><span class="o">/</span><span class="p">(</span>withinAvgVar <span class="o">+</span> btwVar<span class="p">)</span>
    w<span class="p">[</span>w <span class="o">==</span> <span class="kc">Inf</span><span class="p">]</span> <span class="o">&lt;-</span> <span class="m">1</span>
    medianSE <span class="o">&lt;-</span> <span class="kt">vector</span><span class="p">(</span>length <span class="o">=</span> <span class="kp">length</span><span class="p">(</span>psetMAD<span class="p">))</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Beginning Median SE Calculation ... \n&quot;</span><span class="p">)</span>
    <span class="kr">for</span> <span class="p">(</span>k <span class="kr">in</span> <span class="m">1</span><span class="o">:</span><span class="kp">length</span><span class="p">(</span>S<span class="p">))</span> <span class="p">{</span>
        fit <span class="o">&lt;-</span> frmaTools<span class="o">:::</span>rwaFit2<span class="p">(</span>pms<span class="p">[</span>S<span class="p">[[</span>k<span class="p">]],</span> <span class="p">,</span> drop <span class="o">=</span> <span class="kc">FALSE</span><span class="p">],</span> w<span class="p">[</span>S<span class="p">[[</span>k<span class="p">]]],</span>
            probeVec<span class="p">[</span>S<span class="p">[[</span>k<span class="p">]]],</span> psetMAD<span class="p">[</span>k<span class="p">])</span>
        medianSE<span class="p">[</span>k<span class="p">]</span> <span class="o">&lt;-</span> median<span class="p">(</span>fit<span class="o">$</span>StdErrors<span class="p">)</span>
        <span class="kr">if</span> <span class="p">((</span>k<span class="o">%%</span><span class="m">1000</span><span class="p">)</span> <span class="o">==</span> <span class="m">0</span><span class="p">)</span> <span class="p">{</span>
            <span class="kp">message</span><span class="p">(</span><span class="kp">paste</span><span class="p">(</span><span class="s">&quot;Finished probeset:&quot;</span><span class="p">,</span> k<span class="p">,</span> <span class="s">&quot;\n&quot;</span><span class="p">))</span>
            <span class="kp">gc</span><span class="p">()</span>
        <span class="p">}</span>
    <span class="p">}</span>
    <span class="kp">names</span><span class="p">(</span>medianSE<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">names</span><span class="p">(</span>psetMAD<span class="p">)</span>
    <span class="kr">if</span> <span class="p">(</span>verbose<span class="p">)</span>
        <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Median SEs Calculated \n&quot;</span><span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    <span class="kp">rm</span><span class="p">(</span>w<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>pms<span class="p">)</span>
    <span class="kp">rm</span><span class="p">(</span>pns<span class="p">)</span>
    <span class="kp">gc</span><span class="p">()</span>
    <span class="kr">if</span> <span class="p">(</span><span class="kp">is.null</span><span class="p">(</span>cdfname<span class="p">))</span> <span class="p">{</span>
        vers <span class="o">&lt;-</span> <span class="s">&quot;&quot;</span>
    <span class="p">}</span> <span class="kr">else</span> <span class="p">{</span>
        vers <span class="o">&lt;-</span> <span class="kp">as.character</span><span class="p">(</span>packageVersion<span class="p">(</span>cdfname<span class="p">))</span>
    <span class="p">}</span>
    <span class="c1">## vers &lt;- ifelse(!is.null(cdfname), as.character(packageVersion(cdfname)),</span>
    <span class="c1">##     &quot;&quot;)</span>
    <span class="kr">return</span><span class="p">(</span><span class="kt">list</span><span class="p">(</span>normVec <span class="o">=</span> normVec<span class="p">,</span> probeVec <span class="o">=</span> probeVec<span class="p">,</span> probeVarWithin <span class="o">=</span> withinAvgVar<span class="p">,</span>
        probeVarBetween <span class="o">=</span> btwVar<span class="p">,</span> probesetSD <span class="o">=</span> psetMAD<span class="p">,</span> medianSE <span class="o">=</span> medianSE<span class="p">,</span>
        version <span class="o">=</span> vers<span class="p">))</span>
<span class="p">}</span>

<span class="c1">## This reads in the xlsx file for each of the 7 datasets and combines</span>
<span class="c1">## them into one big table of all samples. The Batch column contains</span>
<span class="c1">## the partitioning of samples into unique combinations of Dataset,</span>
<span class="c1">## Scan Date, and Phenotype. Finally, we split based on Tissue type to</span>
<span class="c1">## get one table for biopsies (BX), and one for blood (PAX).</span>
sample.tables <span class="o">&lt;-</span> ddply<span class="p">(</span>datasets<span class="p">,</span> <span class="m">.</span><span class="p">(</span>Dataset<span class="p">),</span> <span class="kr">function</span><span class="p">(</span>df<span class="p">)</span> <span class="p">{</span>
    df <span class="o">&lt;-</span> df<span class="p">[</span><span class="m">1</span><span class="p">,]</span>
    <span class="kp">rownames</span><span class="p">(</span>df<span class="p">)</span> <span class="o">&lt;-</span> <span class="kc">NULL</span>
    dset.dir <span class="o">&lt;-</span> <span class="kp">file.path</span><span class="p">(</span>training.data.dir<span class="p">,</span> df<span class="o">$</span>Dataset<span class="p">)</span>
    x <span class="o">&lt;-</span> read.xlsx<span class="p">(</span><span class="kp">list.files</span><span class="p">(</span>dset.dir<span class="p">,</span> pattern<span class="o">=</span>glob2rx<span class="p">(</span><span class="s">&quot;*.xlsx&quot;</span><span class="p">),</span> full.names<span class="o">=</span><span class="kc">TRUE</span><span class="p">)[</span><span class="m">1</span><span class="p">],</span> <span class="m">1</span><span class="p">)</span> <span class="o">%&gt;%</span>
        setNames<span class="p">(</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;Filename&quot;</span><span class="p">,</span> <span class="s">&quot;Phenotype&quot;</span><span class="p">,</span> <span class="s">&quot;ScanDate&quot;</span><span class="p">))</span>
    x<span class="o">$</span>Filename <span class="o">&lt;-</span> <span class="kp">as.character</span><span class="p">(</span>x<span class="o">$</span>Filename<span class="p">)</span>
    missing.CEL <span class="o">&lt;-</span> <span class="o">!</span>str_detect<span class="p">(</span>x<span class="o">$</span>Filename<span class="p">,</span> <span class="s">&quot;\\.CEL$&quot;</span><span class="p">)</span>
    x<span class="o">$</span>Filename<span class="p">[</span>missing.CEL<span class="p">]</span> <span class="o">&lt;-</span> str_c<span class="p">(</span>x<span class="o">$</span>Filename<span class="p">[</span>missing.CEL<span class="p">],</span> <span class="s">&quot;.CEL&quot;</span><span class="p">)</span>
    <span class="kp">stopifnot</span><span class="p">(</span><span class="kp">all</span><span class="p">(</span>str_detect<span class="p">(</span>x<span class="o">$</span>Filename<span class="p">,</span> <span class="s">&quot;\\.CEL$&quot;</span><span class="p">)))</span>
    parsed.date <span class="o">&lt;-</span> parse.date.from.filename<span class="p">(</span>x<span class="o">$</span>Filename<span class="p">)</span>
    x<span class="o">$</span>ScanDate<span class="p">[</span><span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>parsed.date<span class="p">)]</span> <span class="o">&lt;-</span> parsed.date<span class="p">[</span><span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>parsed.date<span class="p">)]</span>
    x <span class="o">%&gt;%</span> <span class="kp">cbind</span><span class="p">(</span>df<span class="p">)</span> <span class="o">%&gt;%</span>
        <span class="kp">transform</span><span class="p">(</span>Filename<span class="o">=</span><span class="kp">file.path</span><span class="p">(</span>dset.dir<span class="p">,</span> Filename<span class="p">),</span>
                  Batch<span class="o">=</span><span class="kp">droplevels</span><span class="p">(</span>Tissue<span class="o">:</span>Dataset<span class="o">:</span><span class="kp">factor</span><span class="p">(</span>ScanDate<span class="p">)</span><span class="o">:</span>Phenotype<span class="p">))</span> <span class="o">%&gt;%</span>
                      <span class="kp">subset</span><span class="p">(</span><span class="o">!</span> Filename <span class="o">%in%</span> blacklist<span class="p">)</span> <span class="o">%&gt;%</span>
                          <span class="kp">subset</span><span class="p">(</span><span class="o">!</span><span class="kp">duplicated</span><span class="p">(</span>Filename<span class="p">))</span>
<span class="p">})</span> <span class="o">%&gt;%</span>
    <span class="kp">split</span><span class="p">(</span><span class="m">.</span><span class="o">$</span>Tissue<span class="p">)</span> <span class="o">%&gt;%</span>
        <span class="kp">lapply</span><span class="p">(</span><span class="kp">droplevels</span><span class="p">)</span>

<span class="c1">## fRMA requires equal-sized batches, so for each batch size from 3 to</span>
<span class="c1">## 15, compute how many batches have at least that many samples.</span>
x <span class="o">&lt;-</span> <span class="kp">sapply</span><span class="p">(</span><span class="m">3</span><span class="o">:</span><span class="m">15</span><span class="p">,</span> <span class="kr">function</span><span class="p">(</span>i<span class="p">)</span> <span class="kp">sapply</span><span class="p">(</span>sample.tables<span class="p">,</span> <span class="m">.</span> <span class="o">%$%</span> Batch <span class="o">%&gt;%</span> table <span class="o">%&gt;%</span> as.vector <span class="o">%&gt;%</span> <span class="p">{</span><span class="kp">sum</span><span class="p">(</span><span class="m">.</span> <span class="o">&gt;=</span> i<span class="p">)}))</span>
<span class="kp">colnames</span><span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> <span class="m">3</span><span class="o">:</span><span class="m">15</span>

<span class="c1">## Based on the above and the recommendations in the frmaTools paper,</span>
<span class="c1">## I chose 5 as the optimal batch size. This could be optimized</span>
<span class="c1">## empirically, though.</span>
arrays.per.batch <span class="o">&lt;-</span> <span class="m">5</span>

<span class="c1">## For each tissue type, compute fRMA vectors.</span>
vectors <span class="o">&lt;-</span> <span class="kp">lapply</span><span class="p">(</span>sample.tables<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>stab<span class="p">)</span> <span class="p">{</span>
    <span class="kp">set.seed</span><span class="p">(</span><span class="m">1986</span><span class="p">)</span>

    tsmsg<span class="p">(</span><span class="s">&quot;Reading full dataset&quot;</span><span class="p">)</span>
    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>
    tsmsg<span class="p">(</span><span class="s">&quot;Getting reference normalization distribution from full dataset&quot;</span><span class="p">)</span>
    normVec <span class="o">&lt;-</span> normalize.quantiles.determine.target<span class="p">(</span>pm<span class="p">(</span>bg.correct.rma<span class="p">(</span>affy<span class="p">)))</span>
    <span class="kp">rm</span><span class="p">(</span>affy<span class="p">);</span> <span class="kp">gc</span><span class="p">()</span>

    tsmsg<span class="p">(</span><span class="s">&quot;Selecting batches&quot;</span><span class="p">)</span>
    <span class="c1">## Keep only arrays with enough samples</span>
    big.enough <span class="o">&lt;-</span> stab<span class="o">$</span>Batch <span class="o">%&gt;%</span> table <span class="o">%&gt;%</span> <span class="m">.</span><span class="p">[</span><span class="m">.</span><span class="o">&gt;=</span> arrays.per.batch<span class="p">]</span> <span class="o">%&gt;%</span> <span class="kp">names</span>
    stab <span class="o">&lt;-</span> stab<span class="p">[</span>stab<span class="o">$</span>Batch <span class="o">%in%</span> big.enough<span class="p">,]</span> <span class="o">%&gt;%</span> <span class="kp">droplevels</span>

    <span class="c1">## Sample an equal number of arrays from each batch</span>
    subtab <span class="o">&lt;-</span> ddply<span class="p">(</span>stab<span class="p">,</span> <span class="m">.</span><span class="p">(</span>Batch<span class="p">),</span> <span class="kr">function</span><span class="p">(</span>df<span class="p">)</span> <span class="p">{</span>
        df<span class="p">[</span><span class="kp">sample</span><span class="p">(</span><span class="kp">seq</span><span class="p">(</span><span class="kp">nrow</span><span class="p">(</span>df<span class="p">)),</span> size<span class="o">=</span>arrays.per.batch<span class="p">),]</span>
    <span class="p">})</span>

    tsmsg<span class="p">(</span><span class="s">&quot;Making fRMA vectors&quot;</span><span class="p">)</span>
    <span class="c1">## Make fRMA vectors, using normVec from full dataset</span>
    res <span class="o">&lt;-</span> makeVectorsAffyBatch<span class="p">(</span>subtab<span class="o">$</span>Filename<span class="p">,</span> subtab<span class="o">$</span>Batch<span class="p">,</span> normVec<span class="o">=</span>normVec<span class="p">)</span>

    tsmsg<span class="p">(</span><span class="s">&quot;Finished.&quot;</span><span class="p">)</span>
    res
<span class="p">})</span>

<span class="c1">## The code below here just takes the trained vectors and packages</span>
<span class="c1">## them up into installable R packages.</span>
makePackageFromVectors <span class="o">&lt;-</span>
    <span class="kr">function</span> <span class="p">(</span>vecs<span class="p">,</span> <span class="kp">version</span><span class="p">,</span> maintainer<span class="p">,</span> species<span class="p">,</span> annotation<span class="p">,</span>
              packageName<span class="p">,</span> file.dir <span class="o">=</span> <span class="s">&quot;.&quot;</span><span class="p">,</span>
              output.dir <span class="o">=</span> <span class="s">&quot;.&quot;</span><span class="p">,</span> unlink <span class="o">=</span> <span class="kc">TRUE</span><span class="p">)</span>
<span class="p">{</span>
    platform <span class="o">&lt;-</span> <span class="kp">gsub</span><span class="p">(</span><span class="s">&quot;cdf$&quot;</span><span class="p">,</span> <span class="s">&quot;&quot;</span><span class="p">,</span> annotation<span class="p">)</span>
    <span class="c1">## type &lt;- match.arg(type, c(&quot;AffyBatch&quot;, &quot;FeatureSet&quot;))</span>
    <span class="c1">## if (type == &quot;AffyBatch&quot;)</span>
    <span class="c1">##     platform &lt;- gsub(&quot;cdf&quot;, &quot;&quot;, annotation)</span>
    <span class="c1">## if (type == &quot;FeatureSet&quot;) {</span>
    <span class="c1">##     platform &lt;- annotation</span>
    <span class="c1">##     require(oligo)</span>
    <span class="c1">## }</span>
    thispkg <span class="o">&lt;-</span> <span class="s">&quot;frmaTools&quot;</span>
    desc <span class="o">&lt;-</span> packageDescription<span class="p">(</span>thispkg<span class="p">)</span>
    thispkgVers <span class="o">&lt;-</span> desc<span class="o">$</span>Version
    symbolValues <span class="o">&lt;-</span> <span class="kt">list</span><span class="p">(</span>ARRAYTYPE <span class="o">=</span> platform<span class="p">,</span> VERSION <span class="o">=</span> <span class="kp">version</span><span class="p">,</span>
        CREATOR <span class="o">=</span> <span class="kp">paste</span><span class="p">(</span><span class="s">&quot;package&quot;</span><span class="p">,</span> thispkg<span class="p">,</span> <span class="s">&quot;version&quot;</span><span class="p">,</span> thispkgVers<span class="p">),</span>
        FRMATOOLSVERSION <span class="o">=</span> thispkgVers<span class="p">,</span> MAINTAINER <span class="o">=</span> maintainer<span class="p">,</span>
        SPECIES <span class="o">=</span> species<span class="p">)</span>
    createdPkg <span class="o">&lt;-</span> createPackage<span class="p">(</span>packageName<span class="p">,</span> destinationDir <span class="o">=</span> output.dir<span class="p">,</span>
        originDir <span class="o">=</span> <span class="kp">system.file</span><span class="p">(</span><span class="s">&quot;VectorPkg-template&quot;</span><span class="p">,</span> package <span class="o">=</span> <span class="s">&quot;frmaTools&quot;</span><span class="p">),</span>
        symbolValues <span class="o">=</span> symbolValues<span class="p">,</span> unlink <span class="o">=</span> <span class="kp">unlink</span><span class="p">)</span>
    <span class="c1">## if (type == &quot;AffyBatch&quot;)</span>
    <span class="c1">##     vecs &lt;- makeVectorsAffyBatch(files, batch.id, background,</span>
    <span class="c1">##         normalize, normVec, annotation, file.dir, verbose)</span>
    <span class="c1">## if (type == &quot;FeatureSet&quot;)</span>
    <span class="c1">##     vecs &lt;- makeVectorsFeatureSet(files, batch.id, annotation,</span>
    <span class="c1">##         background, normalize, normVec, file.dir, verbose)</span>
    <span class="kp">assign</span><span class="p">(</span>packageName<span class="p">,</span> vecs<span class="p">)</span>
    <span class="kp">save</span><span class="p">(</span><span class="kt">list</span> <span class="o">=</span> <span class="kp">eval</span><span class="p">(</span>packageName<span class="p">),</span> file <span class="o">=</span> <span class="kp">file.path</span><span class="p">(</span>createdPkg<span class="o">$</span>pkgdir<span class="p">,</span>
        <span class="s">&quot;data&quot;</span><span class="p">,</span> <span class="kp">paste</span><span class="p">(</span>packageName<span class="p">,</span> <span class="s">&quot;.rda&quot;</span><span class="p">,</span> sep <span class="o">=</span> <span class="s">&quot;&quot;</span><span class="p">)),</span> compress <span class="o">=</span> <span class="kc">TRUE</span><span class="p">)</span>
<span class="p">}</span>

annotation <span class="o">&lt;-</span> cleancdfname<span class="p">(</span>affyio<span class="o">:::</span>read.celfile.header<span class="p">(</span>sample.tables<span class="p">[[</span><span class="m">1</span><span class="p">]]</span><span class="o">$</span>Filename<span class="p">[</span><span class="m">1</span><span class="p">])</span><span class="o">$</span>cdfName<span class="p">,</span> <span class="kc">FALSE</span><span class="p">)</span>

<span class="kp">dir.create</span><span class="p">(</span><span class="s">&quot;pkgs&quot;</span><span class="p">,</span> <span class="kc">FALSE</span><span class="p">,</span> <span class="kc">TRUE</span><span class="p">,</span> mode<span class="o">=</span><span class="s">&quot;755&quot;</span><span class="p">)</span>

<span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="kp">names</span><span class="p">(</span>vectors<span class="p">))</span> <span class="p">{</span>
    vecs <span class="o">&lt;-</span> vectors<span class="p">[[</span>i<span class="p">]]</span>
    pkgname <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;DSalomon.%s.%sfrmavecs&quot;</span><span class="p">,</span> i<span class="p">,</span> annotation<span class="p">)</span>
    <span class="kp">message</span><span class="p">(</span><span class="s">&quot;Making &quot;</span><span class="p">,</span> pkgname<span class="p">)</span>
    makePackageFromVectors<span class="p">(</span>
        vecs<span class="p">,</span>
        version<span class="o">=</span><span class="s">&quot;0.1&quot;</span><span class="p">,</span>
        maintainer<span class="o">=</span><span class="s">&quot;Ryan C. Thompson &lt;rcthomps@scripps.edu&gt;&quot;</span><span class="p">,</span>
        species<span class="o">=</span><span class="s">&quot;Homo_sapiens&quot;</span><span class="p">,</span>
        annotation<span class="o">=</span>annotation<span class="p">,</span>
        packageName<span class="o">=</span>pkgname<span class="p">,</span>
        output.dir <span class="o">=</span> <span class="s">&quot;pkgs&quot;</span><span class="p">)</span>
<span class="p">}</span>

<span class="kp">save.image</span><span class="p">(</span><span class="s">&quot;train-data.rda&quot;</span><span class="p">)</span>
</pre></div>
</body>
</html>