ryan
/
resume
mirror of https://github.com/DarwinAwardWinner/resume.git


			
							123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306
							<!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">source</span><span class="p">(</span><span class="s">&quot;common.R&quot;</span><span class="p">)</span>
<span class="kn">library</span><span class="p">(</span>BSgenome.Hsapiens.UCSC.hg19<span class="p">)</span>

textConnectionFromLines <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>lines<span class="p">,</span> linesep<span class="o">=</span><span class="s">&quot;\n&quot;</span><span class="p">)</span> <span class="p">{</span>
  <span class="kp">textConnection</span><span class="p">(</span>str_c<span class="p">(</span><span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;%s%s&quot;</span><span class="p">,</span> lines<span class="p">,</span> linesep<span class="p">),</span> collapse<span class="o">=</span><span class="s">&quot;&quot;</span><span class="p">))</span>
<span class="p">}</span>

<span class="c1">## Takes a GRangesList</span>
calculate.block.entropy <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>grl<span class="p">,</span> expr.column<span class="o">=</span><span class="s">&quot;tagExpression&quot;</span><span class="p">)</span> <span class="p">{</span>
  groupfac <span class="o">&lt;-</span> <span class="kp">factor</span><span class="p">(</span><span class="kp">rep</span><span class="p">(</span><span class="kp">names</span><span class="p">(</span>grl<span class="p">),</span> elementLengths<span class="p">(</span>grl<span class="p">)))</span>
  exprs <span class="o">&lt;-</span> <span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span><span class="kp">unlist</span><span class="p">(</span>grl<span class="p">))[[</span>expr.column<span class="p">]])</span>
  total.exprs <span class="o">&lt;-</span> aggregate<span class="p">(</span>exprs<span class="p">,</span> by<span class="o">=</span><span class="kt">list</span><span class="p">(</span>groupfac<span class="p">),</span> FUN<span class="o">=</span><span class="kp">sum</span><span class="p">)</span><span class="o">$</span>x
  qi <span class="o">&lt;-</span> exprs <span class="o">/</span> <span class="kp">rep</span><span class="p">(</span>total.exprs<span class="p">,</span> elementLengths<span class="p">(</span>grl<span class="p">))</span>
  qi.times.log <span class="o">&lt;-</span> qi <span class="o">*</span> <span class="kp">log2</span><span class="p">(</span>qi<span class="p">)</span>
  results <span class="o">&lt;-</span> <span class="o">-</span>aggregate<span class="p">(</span>qi.times.log<span class="p">,</span> by<span class="o">=</span><span class="kt">list</span><span class="p">(</span>groupfac<span class="p">),</span> FUN<span class="o">=</span><span class="kp">sum</span><span class="p">)</span><span class="o">$</span>x
  <span class="kp">names</span><span class="p">(</span>results<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">names</span><span class="p">(</span>grl<span class="p">)</span>
  results
<span class="p">}</span>

<span class="c1">## http://hoffmann.bioinf.uni-leipzig.de/LIFE/blockbuster.html</span>
read.blockbuster.tag.output <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>bbout<span class="p">)</span> <span class="p">{</span>
  x <span class="o">&lt;-</span> <span class="kp">readLines</span><span class="p">(</span>bbout<span class="p">)</span>
  cluster.line.numbers <span class="o">&lt;-</span> <span class="kp">which</span><span class="p">(</span>str_sub<span class="p">(</span>x<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">==</span> <span class="s">&quot;&gt;&quot;</span><span class="p">)</span>
  cluster.lines <span class="o">&lt;-</span> str_sub<span class="p">(</span>x<span class="p">[</span>cluster.line.numbers<span class="p">],</span> <span class="m">2</span><span class="p">)</span>
  cluster.table <span class="o">&lt;-</span> read.table<span class="p">(</span>textConnectionFromLines<span class="p">(</span>cluster.lines<span class="p">),</span>
                              col.names<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;clusterID&quot;</span><span class="p">,</span> <span class="s">&quot;chrom&quot;</span><span class="p">,</span> <span class="s">&quot;clusterStart&quot;</span><span class="p">,</span> <span class="s">&quot;clusterEnd&quot;</span><span class="p">,</span> <span class="s">&quot;strand&quot;</span><span class="p">,</span> <span class="s">&quot;ClusterExpression&quot;</span><span class="p">,</span> <span class="s">&quot;tagCount&quot;</span><span class="p">,</span> <span class="s">&quot;blockCount&quot;</span><span class="p">),</span>
                              colClasses<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;character&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;numeric&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">))</span>
  tag.lines <span class="o">&lt;-</span> x<span class="p">[</span><span class="o">-</span>cluster.line.numbers<span class="p">]</span>
  cluster.tag.line.counts <span class="o">&lt;-</span> cluster.line.numbers<span class="p">[</span><span class="m">-1</span><span class="p">]</span> <span class="o">-</span> cluster.line.numbers<span class="p">[</span><span class="o">-</span><span class="kp">length</span><span class="p">(</span>cluster.line.numbers<span class="p">)]</span> <span class="o">-</span> <span class="m">1</span>
  cluster.tag.line.counts <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span>cluster.tag.line.counts<span class="p">,</span> <span class="kp">length</span><span class="p">(</span>tag.lines<span class="p">)</span> <span class="o">-</span> <span class="kp">sum</span><span class="p">(</span>cluster.tag.line.counts<span class="p">))</span>
  tag.table <span class="o">&lt;-</span> read.table<span class="p">(</span>textConnectionFromLines<span class="p">(</span>tag.lines<span class="p">),</span>
                          col.names<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;tagChrom&quot;</span><span class="p">,</span> <span class="s">&quot;tagStart&quot;</span><span class="p">,</span> <span class="s">&quot;tagEnd&quot;</span><span class="p">,</span> <span class="s">&quot;tagID&quot;</span><span class="p">,</span> <span class="s">&quot;tagExpression&quot;</span><span class="p">,</span> <span class="s">&quot;tagStrand&quot;</span><span class="p">,</span> <span class="s">&quot;blockNb&quot;</span><span class="p">),</span>
                          colClasses<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;character&quot;</span><span class="p">,</span> <span class="s">&quot;numeric&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">))</span>
  tag.table<span class="o">$</span>clusterID <span class="o">&lt;-</span> <span class="kp">rep</span><span class="p">(</span>cluster.table<span class="o">$</span>clusterID<span class="p">,</span> cluster.tag.line.counts<span class="p">)</span>
  tag.table<span class="o">$</span>blockID <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;%s/B%s&quot;</span><span class="p">,</span> tag.table<span class="o">$</span>clusterID<span class="p">,</span> tag.table<span class="o">$</span>blockNb<span class="p">)</span>
  tags <span class="o">&lt;-</span> table.to.granges<span class="p">(</span>tag.table<span class="p">,</span> seqnames.column<span class="o">=</span><span class="s">&quot;tagChrom&quot;</span><span class="p">,</span> start.column<span class="o">=</span><span class="s">&quot;tagStart&quot;</span><span class="p">,</span> end.column<span class="o">=</span><span class="s">&quot;tagEnd&quot;</span><span class="p">,</span> strand.column<span class="o">=</span><span class="s">&quot;tagStrand&quot;</span><span class="p">,</span> seqlengths<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">)</span>
  tags <span class="o">&lt;-</span> <span class="kp">split</span><span class="p">(</span>tags<span class="p">,</span> <span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>tags<span class="p">)</span><span class="o">$</span>blockID<span class="p">))</span>
  <span class="kr">return</span><span class="p">(</span>tags<span class="p">)</span>
<span class="p">}</span>

read.blockbuster.output <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>bbout<span class="p">,</span> bbout.tags<span class="o">=</span><span class="kc">NULL</span><span class="p">)</span> <span class="p">{</span>
  x <span class="o">&lt;-</span> <span class="kp">readLines</span><span class="p">(</span>bbout<span class="p">)</span>
  cluster.line.numbers <span class="o">&lt;-</span> <span class="kp">which</span><span class="p">(</span>str_sub<span class="p">(</span>x<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">==</span> <span class="s">&quot;&gt;&quot;</span><span class="p">)</span>
  cluster.lines <span class="o">&lt;-</span> str_sub<span class="p">(</span>x<span class="p">[</span>cluster.line.numbers<span class="p">],</span> <span class="m">2</span><span class="p">)</span>
  block.lines <span class="o">&lt;-</span> x<span class="p">[</span><span class="o">-</span>cluster.line.numbers<span class="p">]</span>
  cluster.table <span class="o">&lt;-</span> read.table<span class="p">(</span>textConnectionFromLines<span class="p">(</span>cluster.lines<span class="p">),</span>
                              col.names<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;clusterID&quot;</span><span class="p">,</span> <span class="s">&quot;chrom&quot;</span><span class="p">,</span> <span class="s">&quot;clusterStart&quot;</span><span class="p">,</span> <span class="s">&quot;clusterEnd&quot;</span><span class="p">,</span> <span class="s">&quot;strand&quot;</span><span class="p">,</span> <span class="s">&quot;ClusterExpression&quot;</span><span class="p">,</span> <span class="s">&quot;tagCount&quot;</span><span class="p">,</span> <span class="s">&quot;blockCount&quot;</span><span class="p">),</span>
                              colClasses<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;character&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;numeric&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">))</span>
  clusters <span class="o">&lt;-</span> table.to.granges<span class="p">(</span>cluster.table<span class="p">,</span> seqnames.column<span class="o">=</span><span class="s">&quot;chrom&quot;</span><span class="p">,</span> start.column<span class="o">=</span><span class="s">&quot;clusterStart&quot;</span><span class="p">,</span> end.column<span class="o">=</span><span class="s">&quot;clusterEnd&quot;</span><span class="p">,</span> seqlengths<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">)</span>
  block.table <span class="o">&lt;-</span> read.table<span class="p">(</span>textConnectionFromLines<span class="p">(</span>block.lines<span class="p">),</span>
                            col.names<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;blockNb&quot;</span><span class="p">,</span> <span class="s">&quot;blockChrom&quot;</span><span class="p">,</span> <span class="s">&quot;blockStart&quot;</span><span class="p">,</span> <span class="s">&quot;blockEnd&quot;</span><span class="p">,</span> <span class="s">&quot;blockStrand&quot;</span><span class="p">,</span> <span class="s">&quot;blockExpression&quot;</span><span class="p">,</span> <span class="s">&quot;readCount&quot;</span><span class="p">),</span>
                            colClasses<span class="o">=</span><span class="kt">c</span><span class="p">(</span><span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">,</span> <span class="s">&quot;factor&quot;</span><span class="p">,</span> <span class="s">&quot;numeric&quot;</span><span class="p">,</span> <span class="s">&quot;integer&quot;</span><span class="p">))</span>
  block.table<span class="o">$</span>clusterID <span class="o">&lt;-</span> <span class="kp">rep</span><span class="p">(</span>cluster.table<span class="o">$</span>clusterID<span class="p">,</span> cluster.table<span class="o">$</span>blockCount<span class="p">)</span>
  block.table<span class="o">$</span>blockID <span class="o">&lt;-</span> <span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;%s/B%s&quot;</span><span class="p">,</span> block.table<span class="o">$</span>clusterID<span class="p">,</span> block.table<span class="o">$</span>blockNb<span class="p">)</span>
  blocks. <span class="o">&lt;-</span> table.to.granges<span class="p">(</span>block.table<span class="p">,</span> seqnames.column<span class="o">=</span><span class="s">&quot;blockChrom&quot;</span><span class="p">,</span> start.column<span class="o">=</span><span class="s">&quot;blockStart&quot;</span><span class="p">,</span> end.column<span class="o">=</span><span class="s">&quot;blockEnd&quot;</span><span class="p">,</span> strand.column<span class="o">=</span><span class="s">&quot;blockStrand&quot;</span><span class="p">,</span> seqlengths<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">)</span>
  <span class="c1">## blocks. &lt;- split(blocks., as.vector(elementMetadata(blocks.)$clusterID))</span>
  retval <span class="o">&lt;-</span> <span class="kt">list</span><span class="p">(</span>clusters<span class="o">=</span>clusters<span class="p">,</span> blocks<span class="o">=</span>blocks.<span class="p">)</span>
  <span class="kr">if</span> <span class="p">(</span><span class="o">!</span><span class="kp">is.null</span><span class="p">(</span>bbout.tags<span class="p">))</span> <span class="p">{</span>
    retval<span class="o">$</span>tags <span class="o">&lt;-</span> read.blockbuster.tag.output<span class="p">(</span>bbout.tags<span class="p">)</span>
  <span class="p">}</span>
  <span class="kr">return</span><span class="p">(</span>retval<span class="p">)</span>
<span class="p">}</span>


blockbuster.path <span class="o">&lt;-</span> <span class="s">&quot;/home/ryan/bin/blockbuster&quot;</span>

run.blockbuster <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>gr<span class="p">,</span> <span class="kc">...</span><span class="p">)</span> <span class="p">{</span>
  temp.bed.file <span class="o">&lt;-</span> <span class="kp">tempfile</span><span class="p">(</span>fileext<span class="o">=</span><span class="s">&quot;.bed&quot;</span><span class="p">)</span>
  temp.bbout.file <span class="o">&lt;-</span> <span class="kp">tempfile</span><span class="p">(</span>fileext<span class="o">=</span><span class="s">&quot;.bbout&quot;</span><span class="p">)</span>
  temp.bbout.tag.file <span class="o">&lt;-</span> <span class="kp">tempfile</span><span class="p">(</span>fileext<span class="o">=</span><span class="s">&quot;.tag.bbout&quot;</span><span class="p">)</span>
  <span class="kp">tryCatch</span><span class="p">({</span>
    export<span class="p">(</span><span class="kp">sort</span><span class="p">(</span>gr<span class="p">),</span> temp.bed.file<span class="p">,</span> format<span class="o">=</span><span class="s">&quot;BED&quot;</span><span class="p">)</span>
    extra.args <span class="o">&lt;-</span> <span class="kt">list</span><span class="p">(</span><span class="kc">...</span><span class="p">)</span>
    bbargs <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span><span class="kp">rbind</span><span class="p">(</span><span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;-%s&quot;</span><span class="p">,</span> <span class="kp">names</span><span class="p">(</span>extra.args<span class="p">)),</span> extra.args<span class="p">),</span> <span class="s">&quot;-print&quot;</span><span class="p">,</span> <span class="s">&quot;1&quot;</span><span class="p">,</span> temp.bed.file<span class="p">)</span>
    bbargs.tags <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span><span class="kp">rbind</span><span class="p">(</span><span class="kp">sprintf</span><span class="p">(</span><span class="s">&quot;-%s&quot;</span><span class="p">,</span> <span class="kp">names</span><span class="p">(</span>extra.args<span class="p">)),</span> extra.args<span class="p">),</span> <span class="s">&quot;-print&quot;</span><span class="p">,</span> <span class="s">&quot;2&quot;</span><span class="p">,</span> temp.bed.file<span class="p">)</span>
    <span class="kp">system2</span><span class="p">(</span>blockbuster.path<span class="p">,</span> args<span class="o">=</span>bbargs<span class="p">,</span>
            stdout<span class="o">=</span>temp.bbout.file<span class="p">)</span>
    <span class="kp">system2</span><span class="p">(</span>blockbuster.path<span class="p">,</span> args<span class="o">=</span>bbargs.tags<span class="p">,</span>
            stdout<span class="o">=</span>temp.bbout.tag.file<span class="p">)</span>
    x <span class="o">&lt;-</span> read.blockbuster.output<span class="p">(</span>temp.bbout.file<span class="p">,</span> temp.bbout.tag.file<span class="p">)</span>
    x
  <span class="p">},</span> finally<span class="o">=</span><span class="p">{</span>
    <span class="kp">unlink</span><span class="p">(</span><span class="kt">c</span><span class="p">(</span>temp.bed.file<span class="p">,</span> temp.bbout.file<span class="p">,</span> temp.bbout.tag.file<span class="p">))</span>
  <span class="p">})</span>
<span class="p">}</span>

<span class="c1">## Load the reads from the bam file</span>
infile <span class="o">&lt;-</span> <span class="s">&quot;./all-results.bam&quot;</span>
read.ranges <span class="o">&lt;-</span> <span class="p">{</span>
  x <span class="o">&lt;-</span> readAlignedRanges<span class="p">(</span>infile<span class="p">,</span> include.reads<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>
  <span class="c1">## Throw away unneeded columns</span>
  elementMetadata<span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">subset</span><span class="p">(</span>elementMetadata<span class="p">(</span>x<span class="p">),</span> select<span class="o">=-</span><span class="kt">c</span><span class="p">(</span>id<span class="p">,</span>qual<span class="p">,</span>flag<span class="p">))</span>
  read.multi.map.counts <span class="o">&lt;-</span> <span class="kp">table</span><span class="p">(</span><span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span><span class="kp">seq</span><span class="p">))</span>
  elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span>multimap <span class="o">&lt;-</span> Rle<span class="p">(</span><span class="kp">as.vector</span><span class="p">(</span>read.multi.map.counts<span class="p">[</span><span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span><span class="kp">seq</span><span class="p">)]))</span>
  x
<span class="p">}</span>

<span class="c1">## Get needed annotations</span>

<span class="c1">## rRNA &amp; tRNA (from the repeats table)</span>
repeat.table <span class="o">&lt;-</span> get.ucsc.table<span class="p">(</span><span class="s">&quot;rmsk&quot;</span><span class="p">,</span> <span class="s">&quot;rmsk&quot;</span><span class="p">,</span> genome<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">)</span>
repeat.ranges <span class="o">&lt;-</span> table.to.granges<span class="p">(</span>repeat.table<span class="p">,</span> seqnames.column<span class="o">=</span><span class="s">&quot;genoName&quot;</span><span class="p">,</span> start.column<span class="o">=</span><span class="s">&quot;genoStart&quot;</span><span class="p">,</span> end.column<span class="o">=</span><span class="s">&quot;genoEnd&quot;</span><span class="p">,</span> seqlengths<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">)</span>
trna.ranges <span class="o">&lt;-</span> repeat.ranges<span class="p">[</span>elementMetadata<span class="p">(</span>repeat.ranges<span class="p">)</span><span class="o">$</span>repClass <span class="o">==</span> <span class="s">&quot;tRNA&quot;</span><span class="p">]</span>
rrna.ranges <span class="o">&lt;-</span> repeat.ranges<span class="p">[</span>elementMetadata<span class="p">(</span>repeat.ranges<span class="p">)</span><span class="o">$</span>repClass <span class="o">==</span> <span class="s">&quot;rRNA&quot;</span><span class="p">]</span>

<span class="c1">## miRNA &amp; snoRNA</span>
small.rna.table <span class="o">&lt;-</span> <span class="kp">subset</span><span class="p">(</span>get.ucsc.table<span class="p">(</span><span class="s">&quot;wgRna&quot;</span><span class="p">,</span> <span class="s">&quot;wgRna&quot;</span><span class="p">,</span> genome<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">),</span> select<span class="o">=-</span><span class="kt">c</span><span class="p">(</span>thickStart<span class="p">,</span>thickEnd<span class="p">,</span>bin<span class="p">))</span>
small.rna.ranges <span class="o">&lt;-</span> table.to.granges<span class="p">(</span>small.rna.table<span class="p">,</span> seqnames.column<span class="o">=</span><span class="s">&quot;chrom&quot;</span><span class="p">,</span> start.column<span class="o">=</span><span class="s">&quot;chromStart&quot;</span><span class="p">,</span> end.column<span class="o">=</span><span class="s">&quot;chromEnd&quot;</span><span class="p">,</span> seqlengths<span class="o">=</span><span class="s">&quot;hg19&quot;</span><span class="p">)</span>
miRNA.types <span class="o">&lt;-</span> <span class="s">&quot;miRNA&quot;</span>
snoRNA.types <span class="o">&lt;-</span> <span class="kt">c</span><span class="p">(</span><span class="s">&quot;CDBox&quot;</span><span class="p">,</span> <span class="s">&quot;HAcaBox&quot;</span><span class="p">)</span>
miRNA.ranges <span class="o">&lt;-</span> small.rna.ranges<span class="p">[</span> elementMetadata<span class="p">(</span>small.rna.ranges<span class="p">)</span><span class="o">$</span>type <span class="o">%in%</span> miRNA.types <span class="p">]</span>
snoRNA.ranges <span class="o">&lt;-</span> small.rna.ranges<span class="p">[</span> elementMetadata<span class="p">(</span>small.rna.ranges<span class="p">)</span><span class="o">$</span>type <span class="o">%in%</span> snoRNA.types <span class="p">]</span>

<span class="c1">## chrM</span>
chrM.range <span class="o">&lt;-</span> GRanges<span class="p">(</span>seqnames<span class="o">=</span><span class="s">&quot;chrM&quot;</span><span class="p">,</span> IRanges<span class="p">(</span>start<span class="o">=</span><span class="m">1</span><span class="p">,</span>end<span class="o">=</span>seqlengths<span class="p">(</span>read.ranges<span class="p">)[[</span><span class="s">&quot;chrM&quot;</span><span class="p">]]),</span> seqlengths<span class="o">=</span>seqlengths<span class="p">(</span>read.ranges<span class="p">))</span>

<span class="c1">## For each sequence that maps once to an anotated rRNA, tRNA, miRNA,</span>
<span class="c1">## or chrM, remove *all* mappings for that sequence.</span>
forbidden.ranges <span class="o">&lt;-</span>
  <span class="kp">Reduce</span><span class="p">(</span><span class="kp">append</span><span class="p">,</span>
         llply<span class="p">(</span><span class="kt">list</span><span class="p">(</span>rrna.ranges<span class="p">,</span>
                    trna.ranges<span class="p">,</span>
                    <span class="c1">## miRNA.ranges,</span>
                    chrM.range<span class="p">),</span>
               <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> <span class="p">{</span> elementMetadata<span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> <span class="kc">NULL</span><span class="p">;</span> x <span class="p">}))</span>

generate.null.ranges <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>y<span class="p">)</span> <span class="p">{</span>
  x <span class="o">&lt;-</span> GRanges<span class="p">(</span>seqnames<span class="o">=</span><span class="kp">names</span><span class="p">(</span>seqlengths<span class="p">(</span>y<span class="p">)),</span> ranges<span class="o">=</span>IRanges<span class="p">(</span><span class="m">1</span><span class="p">,</span><span class="m">0</span><span class="p">),</span> strand<span class="o">=</span><span class="s">&quot;*&quot;</span><span class="p">,</span> seqlengths<span class="o">=</span>seqlengths<span class="p">(</span>y<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>elementMetadata<span class="p">(</span>y<span class="p">)))</span> <span class="p">{</span>
    elementMetadata<span class="p">(</span>x<span class="p">)[[</span>i<span class="p">]]</span> <span class="o">&lt;-</span> as<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> <span class="kp">class</span><span class="p">(</span><span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>y<span class="p">)[[</span>i<span class="p">]])))</span>
  <span class="p">}</span>
  x
<span class="p">}</span>

select.nearest <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">,</span> y<span class="p">)</span> <span class="p">{</span>
  y <span class="o">&lt;-</span> <span class="kp">append</span><span class="p">(</span>y<span class="p">,</span> generate.null.ranges<span class="p">(</span>y<span class="p">))</span>
  y<span class="p">[</span>nearest<span class="p">(</span>x<span class="p">,</span>y<span class="p">)]</span>
<span class="p">}</span>

annotate.by.granges <span class="o">&lt;-</span> <span class="kr">function</span><span class="p">(</span>peaks<span class="p">,</span> gr<span class="p">,</span> annot.columns<span class="p">)</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>elementMetadata<span class="p">(</span>gr<span class="p">)))</span> <span class="p">{</span>
    elementMetadata<span class="p">(</span>gr<span class="p">)[[</span>i<span class="p">]]</span> <span class="o">&lt;-</span> <span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>gr<span class="p">)[[</span>i<span class="p">]])</span>
  <span class="p">}</span>
  nearest.ranges <span class="o">&lt;-</span> select.nearest<span class="p">(</span>peaks<span class="p">,</span> gr<span class="p">)</span>
  nearest.distance <span class="o">&lt;-</span> distance<span class="p">(</span>peaks<span class="p">,</span> nearest.ranges<span class="p">,</span> ignore.strand<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>
  in.ranges <span class="o">&lt;-</span> Rle<span class="p">(</span>nearest.distance <span class="o">==</span> <span class="m">0</span><span class="p">)</span>
  annot.data <span class="o">&lt;-</span> elementMetadata<span class="p">(</span>nearest.ranges<span class="p">)[</span>annot.columns<span class="p">]</span>
  <span class="kr">if</span> <span class="p">(</span><span class="o">!</span><span class="kp">is.null</span><span class="p">(</span><span class="kp">names</span><span class="p">(</span>annot.columns<span class="p">)))</span> <span class="p">{</span>
    <span class="kp">names</span><span class="p">(</span>annot.data<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">names</span><span class="p">(</span>annot.columns<span class="p">)</span>
  <span class="p">}</span>
  DataFrame<span class="p">(</span>overlap<span class="o">=</span>in.ranges<span class="p">,</span> distance<span class="o">=</span>nearest.distance<span class="p">,</span> annot.data<span class="p">)</span>
<span class="p">}</span>

forbidden.seqs <span class="o">&lt;-</span> <span class="kp">unique</span><span class="p">(</span><span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>read.ranges<span class="p">)</span><span class="o">$</span><span class="kp">seq</span><span class="p">[</span>read.ranges <span class="o">%in%</span> forbidden.ranges<span class="p">]))</span>
forbidden.indices <span class="o">&lt;-</span> <span class="o">!</span><span class="kp">is.na</span><span class="p">(</span>findOverlaps<span class="p">(</span>read.ranges<span class="p">,</span> forbidden.ranges<span class="p">,</span> select<span class="o">=</span><span class="s">&quot;first&quot;</span><span class="p">,</span> ignore.strand<span class="o">=</span><span class="kc">TRUE</span><span class="p">))</span>
forbidden.read.ranges <span class="o">&lt;-</span> read.ranges<span class="p">[</span>forbidden.indices<span class="p">]</span>
read.ranges <span class="o">&lt;-</span> read.ranges<span class="p">[</span><span class="o">!</span>forbidden.indices<span class="p">]</span>

<span class="c1">## Read the counts table</span>
read.counts <span class="o">&lt;-</span> <span class="p">{</span>
  x <span class="o">&lt;-</span> read.csv<span class="p">(</span><span class="s">&quot;./data/R21_R82_read_expr_matrix_no_blanks&quot;</span><span class="p">,</span>
                stringsAsFactors<span class="o">=</span><span class="kc">FALSE</span><span class="p">,</span> sep<span class="o">=</span><span class="s">&quot;\t&quot;</span><span class="p">,</span> header<span class="o">=</span><span class="kc">TRUE</span><span class="p">,</span> row.names<span class="o">=</span><span class="m">1</span><span class="p">)</span>
  <span class="kp">row.names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> str_trim<span class="p">(</span><span class="kp">row.names</span><span class="p">(</span>x<span class="p">))</span>
  <span class="kp">names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">&lt;-</span> str_replace<span class="p">(</span><span class="kp">names</span><span class="p">(</span>x<span class="p">),</span> <span class="s">&quot;_collapsed_sorted_with_zeros$&quot;</span><span class="p">,</span> <span class="s">&quot;&quot;</span><span class="p">)</span>
  x <span class="o">&lt;-</span> x<span class="p">[</span><span class="kp">row.names</span><span class="p">(</span>x<span class="p">)</span> <span class="o">%in%</span> elementMetadata<span class="p">(</span>read.ranges<span class="p">)</span><span class="o">$</span><span class="kp">seq</span><span class="p">,]</span>
  x <span class="o">&lt;-</span> x<span class="p">[</span><span class="kp">order</span><span class="p">(</span><span class="kp">names</span><span class="p">(</span>x<span class="p">))]</span>
  x
<span class="p">}</span>

<span class="c1">## Create per-sample bed files with score = count / multimap</span>
sample.read.ranges <span class="o">&lt;-</span>
  llply<span class="p">(</span><span class="kp">names</span><span class="p">(</span>read.counts<span class="p">),</span>
        <span class="kr">function</span> <span class="p">(</span><span class="kp">sample</span><span class="p">)</span> <span class="p">{</span>
          x <span class="o">&lt;-</span> read.ranges
          <span class="c1">## Score = sample count / multimap</span>
          elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span>score <span class="o">&lt;-</span>
            <span class="kp">as.vector</span><span class="p">(</span>read.counts<span class="p">[</span><span class="kp">as.vector</span><span class="p">(</span>elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span><span class="kp">seq</span><span class="p">),</span> <span class="kp">sample</span><span class="p">]</span> <span class="o">/</span> elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span>multimap<span class="p">)</span>
          <span class="c1">## Eliminate reads with zero score</span>
          x <span class="o">&lt;-</span> x<span class="p">[</span>elementMetadata<span class="p">(</span>x<span class="p">)</span><span class="o">$</span>score <span class="o">&gt;</span> <span class="m">0</span><span class="p">]</span>
          x
        <span class="p">},</span> <span class="m">.</span>parallel<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>
<span class="kp">names</span><span class="p">(</span>sample.read.ranges<span class="p">)</span> <span class="o">&lt;-</span> <span class="kp">names</span><span class="p">(</span>read.counts<span class="p">)</span>

<span class="c1">## Run blockbuster on each sample</span>
<span class="c1">## x &lt;- sample.read.ranges[[1]][1:500]</span>
<span class="c1">## y &lt;- run.blockbuster(x)</span>
<span class="c1">## z &lt;- calculate.block.entropy(y$tags)</span>
blockbuster.results <span class="o">&lt;-</span> llply<span class="p">(</span>sample.read.ranges<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> run.blockbuster<span class="p">(</span>unstranded<span class="p">(</span>x<span class="p">),</span> minBlockHeight<span class="o">=</span><span class="m">5</span><span class="p">),</span> <span class="m">.</span>parallel<span class="o">=</span><span class="kc">TRUE</span><span class="p">)</span>

<span class="c1">## calculate.block.entropy(blockbuster.results[[1]]$tags[1:5])</span>

<span class="c1">## Calculate entropy of blocks</span>
block.entropy <span class="o">&lt;-</span> llply<span class="p">(</span>blockbuster.results<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">)</span> calculate.block.entropy<span class="p">(</span>x<span class="o">$</span>tags<span class="p">),</span> <span class="m">.</span>parallel<span class="o">=</span><span class="kc">TRUE</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>blockbuster.results<span class="p">))</span> <span class="p">{</span>
  elementMetadata<span class="p">(</span>blockbuster.results<span class="p">[[</span>i<span class="p">]]</span><span class="o">$</span>blocks<span class="p">)</span><span class="o">$</span>entropy <span class="o">&lt;-</span> block.entropy<span class="p">[[</span>i<span class="p">]][</span><span class="kp">as.character</span><span class="p">(</span>elementMetadata<span class="p">(</span>blockbuster.results<span class="p">[[</span>i<span class="p">]]</span><span class="o">$</span>blocks<span class="p">)</span><span class="o">$</span>blockID<span class="p">)]</span>
<span class="p">}</span>

<span class="c1">## Plot entropy vs block length</span>
x <span class="o">&lt;-</span> blockbuster.results<span class="p">[[</span><span class="m">1</span><span class="p">]]</span><span class="o">$</span>blocks
y <span class="o">&lt;-</span> <span class="kp">cbind</span><span class="p">(</span><span class="kp">as.data.frame</span><span class="p">(</span>elementMetadata<span class="p">(</span>x<span class="p">)),</span> width<span class="o">=</span>width<span class="p">(</span>x<span class="p">))</span>
ggplot<span class="p">(</span>y<span class="p">,</span> aes<span class="p">(</span>x<span class="o">=</span>width<span class="p">,</span> y<span class="o">=</span>entropy<span class="p">,</span> color<span class="o">=</span><span class="m">..</span>density..<span class="p">))</span> <span class="o">+</span> stat_density2d<span class="p">(</span>geom<span class="o">=</span><span class="s">&quot;tile&quot;</span><span class="p">,</span> aes<span class="p">(</span>fill<span class="o">=</span><span class="m">..</span>density..<span class="p">),</span> contour<span class="o">=</span><span class="kc">FALSE</span><span class="p">)</span> <span class="o">+</span> scale_fill_gradient<span class="p">(</span>low<span class="o">=</span><span class="s">&quot;blue&quot;</span><span class="p">,</span> high<span class="o">=</span><span class="s">&quot;yellow&quot;</span><span class="p">)</span> <span class="o">+</span> scale_color_gradient<span class="p">(</span>low<span class="o">=</span><span class="s">&quot;blue&quot;</span><span class="p">,</span> high<span class="o">=</span><span class="s">&quot;yellow&quot;</span><span class="p">)</span>

<span class="c1">## Compute nearest miRNA/snoRNA to each block and add annotation</span>
block.annot <span class="o">&lt;-</span> llply<span class="p">(</span>blockbuster.results<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>br<span class="p">)</span> <span class="p">{</span>
  annotate.by.granges<span class="p">(</span>br<span class="o">$</span>blocks<span class="p">,</span> small.rna.ranges<span class="p">,</span> <span class="kt">c</span><span class="p">(</span>nearest.ncRNA<span class="o">=</span><span class="s">&quot;name&quot;</span><span class="p">,</span> ncRNA.type<span class="o">=</span><span class="s">&quot;type&quot;</span><span class="p">))</span>
<span class="p">},</span> <span class="m">.</span>parallel<span class="o">=</span><span class="kc">TRUE</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>blockbuster.results<span class="p">))</span> <span class="p">{</span>
  elementMetadata<span class="p">(</span>blockbuster.results<span class="p">[[</span>i<span class="p">]]</span><span class="o">$</span>blocks<span class="p">)[</span><span class="kp">names</span><span class="p">(</span>block.annot<span class="p">[[</span>i<span class="p">]])]</span> <span class="o">&lt;-</span> block.annot<span class="p">[[</span>i<span class="p">]]</span>
<span class="p">}</span>

<span class="c1">## write output</span>
<span class="kp">saveRDS</span><span class="p">(</span>blockbuster.results<span class="p">,</span> <span class="s">&quot;blockbuster_results.RDS&quot;</span><span class="p">)</span>
block.tables <span class="o">&lt;-</span> llply<span class="p">(</span>blockbuster.results<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>br<span class="p">)</span> as<span class="p">(</span>granges.to.dataframe<span class="p">(</span>br<span class="o">$</span>blocks<span class="p">,</span> ignore.strand<span class="o">=</span><span class="kc">TRUE</span><span class="p">,</span> include.width<span class="o">=</span><span class="kc">TRUE</span><span class="p">),</span> <span class="s">&quot;data.frame&quot;</span><span class="p">))</span>
write.xlsx.multisheet<span class="p">(</span>block.tables<span class="p">,</span> <span class="s">&quot;blockbuster_results.xlsx&quot;</span><span class="p">,</span> row.names<span class="o">=</span><span class="kc">FALSE</span><span class="p">)</span>
</pre></div>
</body>
</html>