-
Notifications
You must be signed in to change notification settings - Fork 3
/
dataDrivenEst.html
563 lines (465 loc) · 89.5 KB
/
dataDrivenEst.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<title>for reproducibility</title>
<script type="text/javascript">
window.onload = function() {
var imgs = document.getElementsByTagName('img'), i, img;
for (i = 0; i < imgs.length; i++) {
img = imgs[i];
// center an image if it is the only element of its parent
if (img.parentElement.childElementCount === 1)
img.parentElement.style.textAlign = 'center';
}
};
</script>
<!-- Styles for R syntax highlighter -->
<style type="text/css">
pre .operator,
pre .paren {
color: rgb(104, 118, 135)
}
pre .literal {
color: #990073
}
pre .number {
color: #099;
}
pre .comment {
color: #998;
font-style: italic
}
pre .keyword {
color: #900;
font-weight: bold
}
pre .identifier {
color: rgb(0, 0, 0);
}
pre .string {
color: #d14;
}
</style>
<!-- R syntax highlighter -->
<script type="text/javascript">
var hljs=new function(){function m(p){return p.replace(/&/gm,"&").replace(/</gm,"<")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}};
hljs.initHighlightingOnLoad();
</script>
<style type="text/css">
body, td {
font-family: sans-serif;
background-color: white;
font-size: 13px;
}
body {
max-width: 800px;
margin: auto;
padding: 1em;
line-height: 20px;
}
tt, code, pre {
font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace;
}
h1 {
font-size:2.2em;
}
h2 {
font-size:1.8em;
}
h3 {
font-size:1.4em;
}
h4 {
font-size:1.0em;
}
h5 {
font-size:0.9em;
}
h6 {
font-size:0.8em;
}
a:visited {
color: rgb(50%, 0%, 50%);
}
pre, img {
max-width: 100%;
}
pre {
overflow-x: auto;
}
pre code {
display: block; padding: 0.5em;
}
code {
font-size: 92%;
border: 1px solid #ccc;
}
code[class] {
background-color: #F8F8F8;
}
table, td, th {
border: none;
}
blockquote {
color:#666666;
margin:0;
padding-left: 1em;
border-left: 0.5em #EEE solid;
}
hr {
height: 0px;
border-bottom: none;
border-top-width: thin;
border-top-style: dotted;
border-top-color: #999999;
}
@media print {
* {
background: transparent !important;
color: black !important;
filter:none !important;
-ms-filter: none !important;
}
body {
font-size:12pt;
max-width:100%;
}
a, a:visited {
text-decoration: underline;
}
hr {
visibility: hidden;
page-break-before: always;
}
pre, blockquote {
padding-right: 1em;
page-break-inside: avoid;
}
tr, img {
page-break-inside: avoid;
}
img {
max-width: 100% !important;
}
@page :left {
margin: 15mm 20mm 15mm 10mm;
}
@page :right {
margin: 15mm 10mm 15mm 20mm;
}
p, h2, h3 {
orphans: 3; widows: 3;
}
h2, h3 {
page-break-after: avoid;
}
}
</style>
</head>
<body>
<p>Here we demonstrate a use case for Polyester where the means, variances, and fold changes of the transcripts in the experiment are estimated from real data.</p>
<p>You will need to run <code>polyester_manuscript.Rmd</code> (the original manuscript) first in order to get some of the dependencies for this code. You will also need to download a few more GEUVADIS BAM files (in addition to the ones needed for <code>polyester_manuscript.Rmd</code>. All needed files are listed below):</p>
<ul>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA06985_accepted_hits.bam">NA06985_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA12144_accepted_hits.bam">NA12144_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA12776_accepted_hits.bam">NA12776_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA12778_accepted_hits.bam">NA12778_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA07048_accepted_hits.bam">NA07048_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA12760_accepted_hits.bam">NA12760_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA12889_accepted_hits.bam">NA12889_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20542_accepted_hits.bam">NA20542_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20772_accepted_hits.bam">NA20772_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20815_accepted_hits.bam">NA20815_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20761_accepted_hits.bam">NA20761_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20798_accepted_hits.bam">NA20798_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20518_accepted_hits.bam">NA20518_accepted_hits.bam</a></li>
<li><a href="http://www.ebi.ac.uk/arrayexpress/files/E-GEUV-6/NA20532_accepted_hits.bam">NA20532_accepted_hits.bam</a></li>
</ul>
<p>You will need to run Cufflinks on these BAM files. You can do this with <code>cufflinks_pre_estimation.sh</code> in the <code>data_sim</code> folder.</p>
<p>First we get everything set up and read in the Chromosome 22 transcript abundances estimated with Cufflinks for the GEUVADIS data set.</p>
<pre><code class="r">library(polyester)
library(ballgown)
library(GenomicRanges)
library(limma)
library(EBSeq)
gtfpath = 'chr22.gtf'
seqpath = 'Homo_sapiens/UCSC/hg19/Sequence/Chromosomes'
ceusamps = c('NA06985', 'NA12144', 'NA12776', 'NA12778', 'NA07048', 'NA12760', 'NA12889')
tsisamps = c('NA20542', 'NA20772', 'NA20815', 'NA20761', 'NA20798', 'NA20518', 'NA20532')
allsamps = c(ceusamps, tsisamps)
m1 = read.table('data_sim/abundances/NA06985/isoforms.fpkm_tracking', header=TRUE)
ntx = nrow(m1)
n = length(allsamps)
fpkmMat = matrix(NA, nrow=ntx, ncol=length(ceusamps)+length(tsisamps))
rownames(fpkmMat) = m1$tracking_id
for(i in seq_along(allsamps)){
m1 = read.table(paste0('data_sim/abundances/', allsamps[i], '/isoforms.fpkm_tracking'), header=TRUE)
o = match(rownames(fpkmMat), m1$tracking_id)
stopifnot(all(m1$trackingid[o] == rownames(fpkmMat)))
fpkmMat[,i] = m1$FPKM[o]
}
colnames(fpkmMat) = allsamps
</code></pre>
<p>We will need the transcript lengths in order to get counts from FPKM measurements:</p>
<pre><code class="r">annot = gffReadGR('chr22.gtf', splitByTranscript=TRUE)
names(annot) = substr(names(annot), 2, nchar(names(annot))-1)
transcript_lengths = sapply(width(annot), sum)
o = match(rownames(fpkmMat), names(annot))
transcript_lengths = transcript_lengths[o]
</code></pre>
<p>I updated Polyester's <code>fpkm_to_counts</code> function to accept a matrix rather than a ballgown object. This change has been incorporated into the devel version of Polyester.</p>
<pre><code class="r">fpkm_to_counts = function(bg=NULL, mat=NULL, tlengths=NULL, mean_rps=100e6, threshold=0){
if(is.null(mat)){
tmeas = as.matrix(ballgown::texpr(bg, 'FPKM'))
tlengths = sapply(width(ballgown::structure(bg)$trans), sum)
}else{
tmeas = mat
stopifnot(!is.null(tlengths))
}
index1 = which(rowMeans(tmeas) >= threshold)
tlengths = tlengths[index1]
counts = tlengths*tmeas[index1,]/1000
counts = round(counts*mean_rps/1e6)
return(counts)
}
</code></pre>
<p>Next we take the GEUVADIS FPKM data, convert it into transcript counts, and calculate the “true” fold change for each isoform in the sample. Expression fold changes were calculated between CEU (Europeans living in Utah) and TSI (Tuscans in Italy) populations.</p>
<pre><code class="r">countmat = fpkm_to_counts(mat=fpkmMat, tlengths=transcript_lengths, mean_rps=5e6)
logcountmat = log2(countmat+1)
pop = rep(c('ceu', 'tsi'), each=7)
x = model.matrix(~pop)
fit = lmFit(logcountmat, x)
truebetas = 2^(fit$coefficients[,2])
</code></pre>
<p>We (arbitrarily) will consider a transcript “truly” differentially expressed if its fold change between the populations is above 1.5, in either direction.</p>
<pre><code class="r">isDE = truebetas > 1.5 | truebetas < 0.67
sim_info = data.frame(transcript_id = names(truebetas), fc=as.numeric(truebetas), isDE=isDE)
write.table(sim_info, quote=FALSE, row.names=FALSE, sep='\t', file='sim_info.txt')
</code></pre>
<p>Cufflinks didn't estimate the abundances for one of the annotated transcripts, so we add a zero row to the count matrix.</p>
<pre><code class="r">countmat = rbind(countmat, rep(0, ncol(countmat)))
rownames(countmat)[926] = 'NR_073460_2'
</code></pre>
<p>Next we put the count matrix in the same order as the annotated transripts will be read in the call to <code>simulate_experiment_countmat</code>. </p>
<pre><code class="r">tt = seq_gtf(gtfpath, seqpath)
names(tt) = substr(names(tt), 2, nchar(names(tt))-1)
o = match(names(tt), rownames(countmat))
countmat = countmat[o,]
</code></pre>
<p>Finally we simulate reads based on the count matrix we derived from the GEUVADIS data:</p>
<pre><code class="r">simulate_experiment_countmat(gtf=gtfpath, seqpath=seqpath, readmat=countmat,
outdir='reads', seed=4831)
</code></pre>
<p>Next, we'll need to process these reads to get simulated abundance estimates. You can do this with the <code>tophat.sh</code> and <code>cufflinks_post_estimation.sh</code> scripts in the <code>data_sim</code> folder. The rest of this code relies on the outputs of these scripts.</p>
<p>After processing the simulated data, we read in the estimated FPKMs from the simulation:</p>
<pre><code class="r">fpkmMatSim = matrix(NA, nrow=nrow(fpkmMat), ncol=n)
m1 = read.table('data_sim/abundances_post/sample01/isoforms.fpkm_tracking', header=TRUE)
rownames(fpkmMatSim) = m1$tracking_id
for(i in 1:14){
m1 = read.table(paste0('data_sim/abundances_post/sample', sprintf('%02d', i), '/isoforms.fpkm_tracking'), header=TRUE)
o = match(rownames(fpkmMatSim), m1$tracking_id)
stopifnot(all(m1$trackingid[o] == rownames(fpkmMatSim)))
fpkmMatSim[,i] = m1$FPKM[o]
}
o = match(rownames(fpkmMat), rownames(fpkmMatSim))
fpkmMatSim = fpkmMatSim[o,]
</code></pre>
<p>Now we can correlate the estimated FPKMs from the simulated data with the FPKMs that were used to generate the count matrix in the first place. </p>
<pre><code class="r">sapply(1:14, function(i) cor(fpkmMat[,i], fpkmMatSim[,i]))
</code></pre>
<pre><code>## [1] 0.02381512 0.85779791 0.87300172 0.47599074 0.87604302 0.85955803
## [7] 0.28769016 0.98353719 0.10043093 0.97357947 0.36858699 0.91965147
## [13] 0.35976548 0.15624747
</code></pre>
<p>The correlations are positive, but some are much stronger than others. Basic plots of the correlations make it clear that there are a few transcripts with extremely high estimated FPKMs in the simulated data, and those high correlations really bring the correlations down. Removing 5 outlying transcripts shows that the other 920 transcripts have very strong correlations between the real data used to generate the count matrix and the simulated FPKMs. Below we show the correlations: each box representes 14 correlations, from the 14 replicates, where each correlation is calculated between the real and simulated FPKMs for the Chromosome 22 transcripts.</p>
<pre><code class="r">outliers = c(725, 204, 842, 843, 580)
outlier_cors = sapply(1:14, function(i) cor(fpkmMat[,i], fpkmMatSim[,i]))
no_outlier_cors = sapply(1:14, function(i) cor(fpkmMat[-outliers, i], fpkmMatSim[-outliers, i]))
hasoutliers = rep(c('outliers', 'no outliers'), each=14)
boxplot(c(outlier_cors, no_outlier_cors) ~ hasoutliers, boxwex=0.5, col='gray', ylab='Correlation')
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWgAAAFoCAYAAAB65WHVAAAEJGlDQ1BJQ0MgUHJvZmlsZQAAOBGFVd9v21QUPolvUqQWPyBYR4eKxa9VU1u5GxqtxgZJk6XtShal6dgqJOQ6N4mpGwfb6baqT3uBNwb8AUDZAw9IPCENBmJ72fbAtElThyqqSUh76MQPISbtBVXhu3ZiJ1PEXPX6yznfOec7517bRD1fabWaGVWIlquunc8klZOnFpSeTYrSs9RLA9Sr6U4tkcvNEi7BFffO6+EdigjL7ZHu/k72I796i9zRiSJPwG4VHX0Z+AxRzNRrtksUvwf7+Gm3BtzzHPDTNgQCqwKXfZwSeNHHJz1OIT8JjtAq6xWtCLwGPLzYZi+3YV8DGMiT4VVuG7oiZpGzrZJhcs/hL49xtzH/Dy6bdfTsXYNY+5yluWO4D4neK/ZUvok/17X0HPBLsF+vuUlhfwX4j/rSfAJ4H1H0qZJ9dN7nR19frRTeBt4Fe9FwpwtN+2p1MXscGLHR9SXrmMgjONd1ZxKzpBeA71b4tNhj6JGoyFNp4GHgwUp9qplfmnFW5oTdy7NamcwCI49kv6fN5IAHgD+0rbyoBc3SOjczohbyS1drbq6pQdqumllRC/0ymTtej8gpbbuVwpQfyw66dqEZyxZKxtHpJn+tZnpnEdrYBbueF9qQn93S7HQGGHnYP7w6L+YGHNtd1FJitqPAR+hERCNOFi1i1alKO6RQnjKUxL1GNjwlMsiEhcPLYTEiT9ISbN15OY/jx4SMshe9LaJRpTvHr3C/ybFYP1PZAfwfYrPsMBtnE6SwN9ib7AhLwTrBDgUKcm06FSrTfSj187xPdVQWOk5Q8vxAfSiIUc7Z7xr6zY/+hpqwSyv0I0/QMTRb7RMgBxNodTfSPqdraz/sDjzKBrv4zu2+a2t0/HHzjd2Lbcc2sG7GtsL42K+xLfxtUgI7YHqKlqHK8HbCCXgjHT1cAdMlDetv4FnQ2lLasaOl6vmB0CMmwT/IPszSueHQqv6i/qluqF+oF9TfO2qEGTumJH0qfSv9KH0nfS/9TIp0Wboi/SRdlb6RLgU5u++9nyXYe69fYRPdil1o1WufNSdTTsp75BfllPy8/LI8G7AUuV8ek6fkvfDsCfbNDP0dvRh0CrNqTbV7LfEEGDQPJQadBtfGVMWEq3QWWdufk6ZSNsjG2PQjp3ZcnOWWing6noonSInvi0/Ex+IzAreevPhe+CawpgP1/pMTMDo64G0sTCXIM+KdOnFWRfQKdJvQzV1+Bt8OokmrdtY2yhVX2a+qrykJfMq4Ml3VR4cVzTQVz+UoNne4vcKLoyS+gyKO6EHe+75Fdt0Mbe5bRIf/wjvrVmhbqBN97RD1vxrahvBOfOYzoosH9bq94uejSOQGkVM6sN/7HelL4t10t9F4gPdVzydEOx83Gv+uNxo7XyL/FtFl8z9ZAHF4bBsrEwAAHv9JREFUeAHt3Q20VWWdx/GLyGtxEdNADSGtGMk3Ukua1HEmsQVUjoo1ztIGHbTSVCqzGocZXSxsTWkkjZhNOTUOttI0Syt1FcqY5EthSJov4AtqxKsXLvGi3Du/H5zTOp722Z577j77PPs+32etH3efffbLsz9b/2yeu88+bW00BBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQACBZgnsrg2PaNbG2S4CCCCAQM8EBmrxOcpKpUvpVjYry5TpCg0BBBCITsBXqyG0eerEKGWKskJxcW5XxitzlcHKfKW37VRtIJRj7u2xsD4CCOQjsFq7+UU+u3rtXvq99mXLXj2jPU9UViX04GjNu0w5MeG9nsw6RQt/RvlOT1ZiWQQQiF7gAgn8o/JI3hKhXE16KON45cYEgKmatyZhfk9n+Vi/q3yjpyuyPAIIRC3wDh39bq0QCKVAz9LBL1BmKsuVjcpw5SDFfZys0BBAAIGoBEIp0EukPkHxMMdYxePRvmr2uPMixb80pCGAAAJRCYRSoI2+VVmYoD9O84YqLuKv1yZpgeNqLPQezX9JubbG+8xGAAEEghIIqUDXgpmmN8YoM2otUDH/95rurHhdOXmgXgyrnME0AgggELJAEQr07B4APq9lnaTmQj8y6Q3mIYAAAiEKtOQ3kyFC0CcEEEAgNIEiXEGHZkZ/EAhN4G51aECDnfLvd1wHfOdUo80fAFvb6MqsV1sglALtD5Ck/QfmseUf1j4M3kEgaoETenH0J2vd/RR/mpcWmEAoBXqsXM5X/Ck/f8y7umXxQZXqbfIaAQQQCFoglAL9KSl5PNw5L2gxOocAAgjkJBDSLwkv0TH7AUlvzOnY2Q0CCOx6YuQvgQhTIJQraOv4/mU/kISGAAL5CTyZ367YU08FQrqC7mnfWR4BBBDo0wIU6D59eus+uEFa8mLlZsUPrfJDqmgIINBiAQp0i09AALvvpz7cp/hTln7u7VcVP6/kaIXW9wXeqUP0Q8poAQpQoAM8KTl3aar252+x+azih0k9pLhQX6TQ+r6AH0Z2ZN8/zGIeIQW6mOcty177udvVX+fztOZ5Pg0BBFooQIFuIX4gu/6d+uGr6CEV/fmgpvlwUAUIkwi0QiCk2+xacfzsc9dztu8UhL8P8pPKWxWPSZ6m0Pq+wDId4rN9/zCLeYQU6GKet6x7/XVt8GHlXcpKxb8oTPrIvWbT+pgA90EHfEIp0AGfnJy79ivtz6EhgEAgAoxBB3Ii6AYCCCBQLUCBrhbhNQJxCXAfdMDnmwId8MmhawjkIMB90DkgN7oLCnSjcqyHAAIINFmAXxI2GTjnzfvbzz/Qi336asqfKnylwW34Vr3bGlyX1RBAoEqAAl0FUvCXLqwdvTiGo7TuU73YBrfm9QK/RatyH3SL4OvZLQW6HqXiLONnaXyvF931Jwp/pLzYi22warEEuA864PNFgQ7s5LS3t//fgAED3tSKbm3atGnfoUOHvrt///6v5r3/7u7u/lu2bDlL4ds98sZnf8EKUKDDOzWHrVu3blirurV9+/aWPSRJfzm8X8cdY4F+s4772Fad8xbv99fa/zMt7kOwu6dAh3dqusPrEj1qpsDee+99w4QJE44dPXr0jmbuJ7Rt79ixo+3WW29d3dHR4ee/0BIEKNAJKK2c1dXVtePkk09e24o+rFy5ctDIkSNfGThwYFfe+1+6dOkbVq1atT7v/Yawv0GDBm2ePHnyoLFjx4bQndz6oCG1trvvvnu9CnRu+yzajijQgZ2xzs7OibfccstbWtStz2u/31TWtWj/97Rov+wWgSAFKNDhnZYn1CWnFW26duoxYO7iaIU++0SgSoACXQVS8Jf+steP9eIYDte6X1A6G9zGc1pvfoPrshoCCFQJUKCrQAr+8g/q/629OIYHta4/DdjoL6s29WLfrIoAAlUCFOgqkIK/fFn9f6Dgx0D3EUCgJMDDkvhPAQEEEAhUgAId6ImhWwgggEDIBdp9eyOnCAEEEIhVIJQC3a4T8Fnlx8rxyoeUPyp++M91CoVaCDQEEIhLIJQC7Q9IvFv5ifI15UuKi/RYZaAyTaEhgAACUQmEchfHh6XuAu3nCY9U9lIWK25XKFcp1/sFDQEEEIhFIJQr6McFfoLiJ6n5qV5HKOV2qCZ+U37BTwQQQCAWgVCuoK8U+LeVA5SrFT9u00X7t8r7lL9RaAgggEBUAqEUaA9njFf2VPygnkHKicoeip8PsUWhIYAAAlEJhFKgje7nIJeforZN0/7qJbdxylBliV+8TjtV70+usYzHuP0xZhoCCCBQCIGQCnQtMN/BMUaZUWuBivk/1/SvK15XTl6qFx46oSGAAAKFEChCgZ7dA8kNWtZJahs1c0jSG8xDAAEEQhQI5S6OShv/pTGicgbTCCCAQIwCoRRofxhljrJS2a6sV3xP9DLFvySkIYAAAtEJhDLEMU/yo5QpygrFxbld8Z0dc5XBynyFhgACCEQjEMoV9CSJn6ssVToV39Hhb5JcrFyonKTQEEAAgagEQinQHsrwQ5KS2lTNXJP0BvMQQACBviwQyhDHLCEvUGYqyxXfcTFc8XfsuY+17m3WWzQEEECgbwqEUqCXiHeCMlEZq3g82lfNHndepHjIg4YAAghEJRBKgTb6VmVhVPocLAIIIJAiEMoYdEoXeQsBBBCIU4ACHed556gRQKAAAhToApwkuogAAnEKUKDjPO8cNQIIFECAAl2Ak0QXEUAgTgEKdJznnaNGAIECCFCgC3CS6CICCMQpQIGO87xz1AggUAABCnQBThJdRACBOAUo0HGed44aAQQKIECBLsBJoosIIBCnAAU6zvPOUSOAQAEEKNAFOEl0EQEE4hSgQMd53jlqBBAogAAFugAniS4igECcAhToOM87R40AAgUQoEAX4CTRRQQQiFMgpG9UifMMcNTRC3R1dfVfvXp124ABA6Ky2Lx5c5uPPaqD7uHBUqB7CMbiCGQtsGXLlj3nzJmT9WYLsb3hw4ePKERHW9RJhjhaBM9uEUAAgdcT4Ar69YR4H4EmCwwZMmTdzJkz20aPHt3kPYW1+c7OzrYrrrhifUdHR1gdC6g3FOiATgZdiVNgt91262pvb28bMSKuf+3vvvvubf369euK86zXd9QMcdTnxFIIIIBA7gIU6NzJ2SECCCBQnwAFuj4nlkIAAQRyFwi5QA+WRlw3huZ++tkhAgiELBBKgd5fSN9VjlT2Vr6lrFJeVr6tDFRoCCCAQFQCaQV6D0ksUB5VnqjIXE1n3S7XBp9Xfqd8SvHdJQcrhyrDlH9VaAgggEBUAmm32X1OEsOVC5TOCpX1FdNZTR6rDf2Vsl35e+Uk5QXFzcX52p1T/IEAAghEJJB2Bb2fHK5RFioPVWS5prNuT2qDZ5Y2eo9+Ti5N+8dU5amK10wigAACUQikXUHfIoEzFBfn1U3WOE/bv105W3la+YpyltKltCu+wqYhgAACUQmkFeh9JeEr2WnKCmWH4vYz5aKdU9n94avy8coJyjjF49EbFF8536G8qtAQQACBqATSCrSvaB9WPA79ZsVFc5vSjDFobbatW7mrFL+mIYAAAlELpBXolyRzufJPiouy7+r4oXK6kmfzFfVQZUkdO/WwyGk1lvN2nqvxHrMRQACB4ATSCvS56u3bFA89PK4coHhs+BLFhTuv5iGWMcqMOnb4HS3zvzWW+w/N978EaAgggEAhBNIK9NE6gi8rLs5uHod2Yb7SL5rY3Cff++wxaLfZu37U9afHyctj5dUr+BeOHkahIYAAAoUQ2C2ll/frvWOq3vfrNVXzsnjpTwrOUVYq2xUPqWxWlinTFRoCCCAQnUDaFfRN0vCnCI9T7lOOUA4vvdaPTNs8bW2UMkXxlbqLs2+v8/DKXMXP5Ziv0BBAAIFoBNKuoNdJ4RDlBsXL+Xa3g5VHlKzbJG3QY95LFX9q0UMR/pqFxcqFykkKDQEEEIhKIO0K2hAu0lfnIOKhjOOVGxP2NVXzmjGskrArZiGAAALhCCQV6AfUvc8p71WSxn9/qvm+qs2yzdLGFigzleXKRmW4cpDiPk5WaAgggEBUAkkF+hwJPKv4cZ+/UCo/xTdQr5vxjOYl2u4EZaIyVvF4tK+aPe68SOHuCyHQEEAgLoGkAu3b6jzm/EnFH+teqJSbx4LPUFw0s25btcHKfWW9fbaHAAIIFEogqUCfpSMo3zFxQdXRbNLrz1fN4yUCCCCAQBMEfKVc3a7VDA9jfEbxkIOnHRdz3/p2jUJDAAEEEGiyQNIVtHfpceeraux7iOZvqfEesxFAAAEEMhKoVaC9+b0UX02/Xemv+GrbHxj5lZL3A5O0SxoCCCAQl0DSEEdZwLe8DVW+qbyg+FY43/42R6EhgAACCDRZIK1AH6h9X6n8t+Kvv7pZ8X3RHpumIYAAAgg0WSCtQL+ofe+v+KPXvv/5Tcp6xfNoCCCAAAJNFkgbg/6W9u1nYfg7An+k3K64UPshSjQEEEAAgSYLpBXox7TvccoOxYX648rLyvcVGgIIIIBAkwXSCrR37Y97l9vXyxP8RAABBBBovkBSgV6k3e6Rsuu79N5nU97nLQQQQACBDASSCrSLrz85WNn21As/n9nDHWsr32AaAQQQQKA5AkkF+sHSrnyHxxeV05R+igv3J5QZCg0BBBBAoMkCabfZnaN9/61ycqkPfvSob73zfBoCCCCAQJMF0gr0Mdr3V5SXSn14RT/nKi7aNAQQQACBJgukFWh/w7aLdGX7sF78oXIG0wgggAACzRFIGoMu7+mrmnhIOUHZR/G90GOV9ys0BBBAAIEmC6QVaD8YabzyEcUf7763FN/JQUMAAQQQaLJAWoGeo33/UflSk/vA5hFAAAEEEgTSxqCf0/KHKH4WNA0BBBBAIGeBtCtof2vKVMVDHf6FYXlo405Nf1qhIYAAAgg0USCtQPsbtk9M2Pe6hHnMQgABBBDIWCCtQPtTg4xBZwzO5hBAAIF6BRiDrleK5RBAAIGcBdKuoBmDzvlksDsEEECgUiCtQP9MC/62cuHSNGPQCSjMQgABBLIWSCvQvs3O8XcRvkV5SvmTklcbrB0NVHwXCQ0BBBCITiBtDNr3P1+v+PnPfpLdJuUHyiAlj3aKdnJVHjtiHwgggECIAmkF+lx1+G2KP+7tq+i3K/2US5Ssm6/ON1TlOr0+ozTPf1HQEEAAgagE0gr00ZL4svJ4SWSFfl6uHFd6neWP6drYGsUPaDq8FH9ZwK2l6Yv1k4YAAghEJZBWoO+XRPXjRv3ahTTrdp82eKTiK3YPa2xWPLTSqXgc3NM0BBBAICqBtF8S3iSJRxVfMbuAHqH46rYZV9Da7M5fBp6pn/6KrUXKA0r54+WapCGAAAJxCaRdQft2ukOUGxQvd4dysPKI0sz2fW18krKXsqqZO2LbCCCAQMgCaVfQ7reL9NWK7+hwkfbXXuXRXtBOPpjHjtgHAgggEKpAUoH2nRrnKyOVS0sdf49+flPx/IWleXn9GKcdDVWW1LHDc7TM6TWWO1Dzn6nxHrMRQACB4ASSCrSL8qnK2RW9fVDT8xSPS39AeVjJq03TjsYoM+rYoW/Nc5Ka7xDxXzo0BBBAoBACSQX6FPX8o0r59jofyKvKtYrHhV28m1mg3adhygbFbfauH/yJAAIIxCVQ/UtCF8cDlCdrMNyj+UfVeK83s/2R7jnKSmW7sl7xrXbLFN8jTUMAAQSiE6gu0L5SfkzxmHNS83wPd2TdPHzyTmWK0q64X/sqHtb4uOJnU9MQQACBqASqC7QP/nrlW8phflFqXs5jwR6f/nFpXpY/fFudP1q+VOlUupUOZbFyoXKSQkMAAQSiEkgag/6GBIYovlL2bXZrFD+HY5PiQnm/knXzUMbxyo0JG56qee4DDQEEEIhKIKlAG2CuMl95l+I7KHxl+4TSrE/2zdK2FygzleXKRmW4cpDiPk5WaAgggEBUArUKtBG2KR5icJrdlmgHE5SJylhllOKrZv8lsUjxkAcNAQQQiEogrUDnDbFVO1yY907ZHwIIIBCqQNIvCUPtK/1CAAEEohKgQEd1ujlYBBAokgAFukhni74igEBUAhToqE43B4sAAkUSoEAX6WzRVwQQiEqAAh3V6eZgEUCgSAIU6CKdLfqKAAJRCVCgozrdHCwCCBRJgAJdpLNFXxFAICoBCnRUp5uDRQCBIglQoIt0tugrAghEJUCBjup0c7AIIFAkAQp0kc4WfUUAgagEKNBRnW4OFgEEiiQQ0uNGi+RGXxHITKCrq6t/R0dH2xNP+Dsx4mmvvvpqW3d3d/94jrjnR0qB7rkZayCQqcD27duvu+aaazb379/fz0TPtW3ZsmWUCuXgYcOGPZvrjrUzF2f95XRf3vst0v4o0EU6W/S1TwqsXbv2dh2Y04p2sna636ZNm+a1YufsM12AAp3uw7sI9HWBn+kAqQOBnmVOTKAnhm4hkJPAn3LaD7tpQIC7OBpAYxUEEEAgDwEKdB7K7AMBBBBoQIAC3QAaqyDQhwQm6VjO7EPH06cOhQLdp04nB4NAjwXeqDWG93gtVshFgAKdCzM7QQABBHouQIHuuRlrIIAAArkIcJtdLszsBIFgBbgPOthTU4wb1P1Zff9Fsi1gR7qGQFEFuA864DMXyhDHaBl9V+lU7lbeppTbNE38T/kFPxFAAIFYBEIp0DMF/gflSGWxskh5h0JDAAEEohUIZQx6ss7ABGWLMkt5TLlTeZ9CQwCB5gn4PuhRiv8FSwtMIJQraBdkXz2X2/c04adr/VR5U3kmPxFAIHMB7oPOnDS7DYZSoK/VId2kXFJxaFdp+gfKVyvmMYkAAghEIxDKEMddEj9QOaBK/jK9vrf0XtVbvEQAAQT6tkAoBdrKm5VHE7jv0TyHhgAC2QtwH3T2ppltMaQCXeugxumNocqSWgtUzD9H06dXvK6c9BX6M5UzmEYAgTbugw74P4IiFGjfBz1GmVGH43VaxklqHssemfQG8xBAAIEQBUIs0O7TMGVDCWx2iHD0CQEEEGi2QCh3cQzUgc5RVirblfWKx6SXKdMVGgIINEdgkjbL86CbY9vrrYZyBe17nn2z/BRlheLi3K6MV+Yqg5X5Cg0BBLIV4D7obD0z3VooV9D+W/xcZani53F0Kx2KP/Z9oXKSQkMAAQSiEgilQHso4/ga8lM1f02N95iNAAII9FmBUIY4/PyNBYofmrRc2agMVw5S3Ec/q4OGAALZC3AfdPammW0xlAK9REfkhyVNVMYqHo/2VbPHnRcpHvKgIYBA9gLcB529aWZbDKVA+4C2KgszOzI2hAACCBRcIJQx6IIz0n0EEEAgewEKdPambBGBIglwH3TAZ4sCHfDJoWsI5CDAfdA5IDe6Cwp0o3KshwACCDRZgALdZGA2jwACCDQqENJdHI0eA+shgEDjAtwH3bhd09ekQDedmB0gELQA90EHfHoY4gj45NA1BBCIW4ACHff55+gRQCBgAQp0wCeHriGQgwD3QeeA3OguKNCNyrEeAn1DgPugAz6PFOiATw5dQwCBuAUo0HGff44eAQQCFuA2u4BPDl1DIAcB7oPOAbnRXVCgG5VjPQT6hgD3QQd8HhniCPjk0DUEEIhbgAId9/nn6BFAIGABCnTAJ4euIZCDAPdB54Dc6C4o0I3KsR4CfUOA+6ADPo8U6IBPDl1DAIG4BSjQcZ9/jh4BBAIW4Da7gE8OXUMgBwHug84BudFdUKAblWM9BPqGAPdBB3weGeII+OTQNQQQiFuAAh33+efoEUAgYAEKdMAnh64hkIMA90HngNzoLkIs0B4XH9HoAbEeAgj0SID7oHvEle/CoRTogTrsOcpKZbuyXtmsLFOmKzQEEEAgOoFQ7uKYJ/lRyhRlheLi3K6MV+Yqg5X5Cg0BBP5S4K81q9GLrYO07l7KMX+52brnPKAlfWFFy1gglALtcbCJyqqK4+vQ9GLlQuUyhQItBBoCCQJ/p3mN/r/s9bYpJyRst95Zj2hBCnS9Wj1YrtGT2oNd1LWohzKOV25MWHqq5q1JmM8sBBDYJXA5EH1TIJQCPUu8C5SZynJlozJc8T+/3MfJCg0BBBCISiCUAr1E6hMUD3OMVTwe7atmD2ssUroVGgIIIBCVQCgF2uhblYW91O+n9Wv9ssTv0RBAAIHCCIRUoGuhjdMbQxVfZb9eO1sLfLTGQm/V/CdqvMdsBBBAIDiBIhToaVIbo8yoQ++/tIyT1D6imXwAJkmGeQggEKRAiAXafRqmbCiJzQ5Sjk4hgAACTRYIZVzWnyT8d+UMZT/F/fJjEJ9RrlSuV3rbDtcG7lDqGSrp7b6Kuv7R6rjtu4t6APS7xwK+IPIndx/r8ZrxrHCADtX3ib8YzyG/9ki/oZe3KYcqfjaAi8RwxXd1+FNKn1BozRe4QbvwX5C0eAQ+pEP17a00BGoK+ErZt9YlNV/V3Zn0BvMyF6BAZ04a/AYp0AGfolq3pOXd5fInCZP2yycJk1SYhwACfV7A408hND5JGMJZoA8IIBCUQCgF2r+445OEQf2nQWcQQKDVAqEUaDtsVRa2GoT9I4AAAqEIhDIGHYoH/UAAAQQQQCBIgb3Vq/5B9oxONUvAj1HwLa00BBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAAEEEEAAAQQQQAABBBBAAIEYBS7TQfvbbtxeUvwFvkcofPmuEAraOKcFPXF0G4FKAX/k2199Nbg0s1ygXbD3Kc3jR7EEOKfFOl9tPCwp/BN2iLr4HeVryjrlEcVfDebm83ee4se1+vvS/k1JOqdpyy3SOm9Vyu0+TYxRvlea8Vv93Ks07R9vV64vve6nn5cqLyje/78onue2UPm88kflA4q/b/J5xcdwkzJCoTUuwDlt3I41EchM4ChtaYfi4usr1/nK7Yrb+Yq/jebdynuVJ5V/Vqpb2nLlK+PyOqs04SK8h9Kt+KvIXHTLy1UOcZyp+b9X/Cxv99N9eY/itlK5S/E34uyvbFL8xb3e7k+ULyi0xgU4p43bsSYCmQm48G2o2NpxmnYhdLtfuWjn1K4/vqgf91a8Lk+mLVcuvOVlywW61j+HKwv0z7XSJUp7KZ/Wz9mKmwv05J1TbW2D9HOz4vdHKgMVWu8EOKe98yvE2v5nEi18gdUVXXSh2730eox+Lq54z9P7VrwuT9a7XHn5en/6G8AvVvxLQ8fTvpouNxdpt23KacrHFA+F3KH4F460xgU4p43bFWZNCnQxTlV3jW6u1fx3Vrzn8eoVFa/Lk2nLefjEV7huvmrec+dUfX88pMU8VLFPKR4aOV0pN2/bzf+d/UY5rJSN+vmfCq1xAc5p43aFWZMCXZhTldjROzX3HxQ/cN2FdZryS6W6pS3nX+IdXVrhVP0cUJp2cd2meNu12m16Y7riX/h5nPoGZaZS3fbSDA/LvEX5nfJThdY7Ac5p7/xYG4FMBI7SVvyLuHI7UhPl1/6F2+3KesXDIC6Q5eEPTf65pS13ipbqVJ5RblYeVXwl7OZC+oriq/TyWHXlGLS/jeMWpUN5WvmxMkRx8/DG+J1Tu/5w4X5WcYF+TvFx0RoX4Jw2bseaCOQq4F/SlQtj2o5rLeerZl8FJ7U3JM2smudl6lnOq+1dtS4veyfAOe2dH2sjgAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCCCAAAIIIIAAAggggAACCAQj8P+mQu3jvHxvHgAAAABJRU5ErkJggg==" alt="plot of chunk nooutliers"/> </p>
<p>So we have shown that true and simulated FPKMs correlate with each other, which we also show in the main manuscript. We will now show that the estimated differential expression fold changes correlate between the real and simulated data, and that you can use Polyester to evaluate statistial methods for differential expression. Below we estimate differential expression status and fold change between the two populations (CEU/TSI) using EBSeq:</p>
<pre><code class="r">Conditions = rep(c('CEU', 'TSI'), each=7)
IsoformNames = rownames(fpkmMatSim)
iso_gene_relationship = read.table('data_sim/abundances_post/sample01/isoforms.fpkm_tracking', header=TRUE, colClasses=c('character', 'NULL' ,'NULL', 'character', rep('NULL', 9)))
iso_gene_relationship = iso_gene_relationship[match(IsoformNames, iso_gene_relationship$tracking_id),]
sum(IsoformNames != iso_gene_relationship$tracking_id) # expect 0
</code></pre>
<pre><code>## [1] 0
</code></pre>
<pre><code class="r">IsosGeneNames = iso_gene_relationship$gene_id
IsoSizes = MedianNorm(fpkmMatSim)
NgList = GetNg(IsoformNames, IsosGeneNames)
IsoNgTrun = NgList$IsoformNgTrun
IsoEBOut = EBTest(Data=fpkmMatSim, NgVector=IsoNgTrun,
Conditions=as.factor(Conditions), sizeFactors=IsoSizes, maxround=20)
</code></pre>
<pre><code>## Removing transcripts with 75 th quantile < = 10
## 572 transcripts will be tested
</code></pre>
<pre><code>## iteration 1 done
##
## time 4.84
##
## iteration 2 done
##
## time 3.1
##
## iteration 3 done
##
## time 1.12
##
## iteration 4 done
##
## time 1.27
##
## iteration 5 done
##
## time 0.88
##
## iteration 6 done
##
## time 1.31
##
## iteration 7 done
##
## time 1.11
##
## iteration 8 done
##
## time 0.84
##
## iteration 9 done
##
## time 0.94
##
## iteration 10 done
##
## time 0.78
##
## iteration 11 done
##
## time 1.03
##
## iteration 12 done
##
## time 0.92
##
## iteration 13 done
##
## time 0.97
##
## iteration 14 done
##
## time 1
##
## iteration 15 done
##
## time 0.94
##
## iteration 16 done
##
## time 0.84
##
## iteration 17 done
##
## time 1.06
##
## iteration 18 done
##
## time 0.84
##
## iteration 19 done
##
## time 0.78
##
## iteration 20 done
##
## time 1.08
</code></pre>
<pre><code class="r">fold_changes = PostFC(IsoEBOut, SmallNum=1)
fold_changes$Direction #CEU over TSI, so need to flip (CEU was reference in my estimate of "true" fold change)
</code></pre>
<pre><code>## [1] "CEU Over TSI"
</code></pre>
<pre><code class="r">o = match(names(fold_changes$PostFC), sim_info$transcript_id)
true_fc = 1/sim_info$fc[o]
sum(sim_info$transcript_id[o] != names(fold_changes$PostFC))
</code></pre>
<pre><code>## [1] 0
</code></pre>
<pre><code class="r">plot(log2(true_fc), log2(fold_changes$PostFC), xlab='True fold change', ylab='EBSeq estimated fold change', main='True vs. Estimated Fold Changes (log2 scale)')
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWgAAAFoCAYAAAB65WHVAAAEJGlDQ1BJQ0MgUHJvZmlsZQAAOBGFVd9v21QUPolvUqQWPyBYR4eKxa9VU1u5GxqtxgZJk6XtShal6dgqJOQ6N4mpGwfb6baqT3uBNwb8AUDZAw9IPCENBmJ72fbAtElThyqqSUh76MQPISbtBVXhu3ZiJ1PEXPX6yznfOec7517bRD1fabWaGVWIlquunc8klZOnFpSeTYrSs9RLA9Sr6U4tkcvNEi7BFffO6+EdigjL7ZHu/k72I796i9zRiSJPwG4VHX0Z+AxRzNRrtksUvwf7+Gm3BtzzHPDTNgQCqwKXfZwSeNHHJz1OIT8JjtAq6xWtCLwGPLzYZi+3YV8DGMiT4VVuG7oiZpGzrZJhcs/hL49xtzH/Dy6bdfTsXYNY+5yluWO4D4neK/ZUvok/17X0HPBLsF+vuUlhfwX4j/rSfAJ4H1H0qZJ9dN7nR19frRTeBt4Fe9FwpwtN+2p1MXscGLHR9SXrmMgjONd1ZxKzpBeA71b4tNhj6JGoyFNp4GHgwUp9qplfmnFW5oTdy7NamcwCI49kv6fN5IAHgD+0rbyoBc3SOjczohbyS1drbq6pQdqumllRC/0ymTtej8gpbbuVwpQfyw66dqEZyxZKxtHpJn+tZnpnEdrYBbueF9qQn93S7HQGGHnYP7w6L+YGHNtd1FJitqPAR+hERCNOFi1i1alKO6RQnjKUxL1GNjwlMsiEhcPLYTEiT9ISbN15OY/jx4SMshe9LaJRpTvHr3C/ybFYP1PZAfwfYrPsMBtnE6SwN9ib7AhLwTrBDgUKcm06FSrTfSj187xPdVQWOk5Q8vxAfSiIUc7Z7xr6zY/+hpqwSyv0I0/QMTRb7RMgBxNodTfSPqdraz/sDjzKBrv4zu2+a2t0/HHzjd2Lbcc2sG7GtsL42K+xLfxtUgI7YHqKlqHK8HbCCXgjHT1cAdMlDetv4FnQ2lLasaOl6vmB0CMmwT/IPszSueHQqv6i/qluqF+oF9TfO2qEGTumJH0qfSv9KH0nfS/9TIp0Wboi/SRdlb6RLgU5u++9nyXYe69fYRPdil1o1WufNSdTTsp75BfllPy8/LI8G7AUuV8ek6fkvfDsCfbNDP0dvRh0CrNqTbV7LfEEGDQPJQadBtfGVMWEq3QWWdufk6ZSNsjG2PQjp3ZcnOWWing6noonSInvi0/Ex+IzAreevPhe+CawpgP1/pMTMDo64G0sTCXIM+KdOnFWRfQKdJvQzV1+Bt8OokmrdtY2yhVX2a+qrykJfMq4Ml3VR4cVzTQVz+UoNne4vcKLoyS+gyKO6EHe+75Fdt0Mbe5bRIf/wjvrVmhbqBN97RD1vxrahvBOfOYzoosH9bq94uejSOQGkVM6sN/7HelL4t10t9F4gPdVzydEOx83Gv+uNxo7XyL/FtFl8z9ZAHF4bBsrEwAAQABJREFUeAHtnQW4HEXatn+I4QnuJEAI7u7BYXH3xWWBxWWRheCLs7gGl8XdCSG4a7AACUkgOMESLPDfd9LF9s43M0fnzMzJ+1zXnaqu7q6ufrr77Xdqzjn5f/8vFA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA7814EJ/lutydpSjGqGMiN7l3XSXrQIJzJbmZO5m3V/lFmfXzUhC79nDdNTLg0j4IWsrVJF/rjNOcby7DQ1DICRRTpYh7ZORdptegRGlViXb26MH7Ozw4LwAQzM71xQ78Ky96nbemz9zW/v/ev6T+BFaI/alZMaBg/CeuA9cBe0tiaiwxXhR9DLX6Aaasz9kx/XDix8AfflG9tD/Q5OwoBUij7t4SRz59C3zLnqQefctqWqvnS3g5tzG/jQuP+tubbWrnalwzNh7xZ2bGB2rL5QismgXep+6FFshyJtjfHj79lxTi+yf2ryhfpmtl0aky/Fy6EDqA3AdfnrYXt70ZKciOe3fnZCo7Nlg3RrSS99Nn6A5LOJWU+ohhpz/+TH1YeFr2CKfGNj6h0bs1EVt/FGfyw7/raU3gx35tqezda1t8JM+dEiJ/VbkbbCpsVpuAaeya0YRP0kMJhUSkfS8QEZlTpGvt8TWTAryevr/EKF693o/2GYBnzxmUwYsHeDnUH/L4P2Lu+rL+GBCp6onu4Eb8MZ4MtgQzgH/gK1rmsZ4DFwIPSBRqvWA7SBKmkxKgZoM6x/p0bKBcGLZCA3OM0Gl8JCMDfcAh+A2h2mhAvge1C9YGvwQXsJDG5joFDL07ACuM0juZUGpc7gmH6FrcBxWvcF4riaGjieZJ/8ObL4P7L/jcEg8So8DmYUk8OOoGaGw+BCMLB/B2Y3yv1XhXuhO6wJQ+BimBp8GCaFW+F5SHK/lWAOeB/ugUFg+9KgVgPP92oXUEP+LsA2+u+0xXXQWF3Jho6hlBzjWmD/budDUhjQafpTk1HbA2aHJ/5sLV05nFXeM3fCZtlmHsN7zXXzZG35wgDu/WGW7bkOhKRZqHhNnA75HLx3HgI1F2wCz8APsCl0gOvhdUiaj8q24LobwT7nh6vgU1Dlrof7NeX+9XlbHc4F7/dSasy18H5eBYbCZbAreK4XgT6/CcfDTWDyYoBeFErJe7LYM5K2n5CK951+j4CH4WX4A1S56zFui//9d2IWdwGvwWDw+n4C6n3w2h0AZ8G30O7kTaZ5Bxac2c5Z+7tZ6TYGZwOz9fUgyWBi26xZw6qUP2Zto7LycUpv1EJ5M7rv27kVy2VtPkzqOHAbL/hHWf15yo7QGPkxzv29QdcqYEGWlefmNo7bsfwGP4GByPNyXZ7uLOuBbQZcpYcuvwbO4xm4XX4CfJDznqTj+kC4jcf7PquPpJwTUn+uF89ZNeTvymyTjm3peXwJ9rE0FJPHdL37ds8xBfUkXyLpHNxWDM4+jKrQDx+ut8HtfoYx8HG2fDplMT1No9vrSzltwEq3Gwqe43fZsuPrAWpK0He3+yYrrf8N1Ebg8gugR8l/6wZvtQKke9jjWE/PxJLUVUPXo6n3r+NzXPqZ5LFtmzBraOhauNkp4D7i/m+A99lgKKaNaXTbAcVW0lbuGXEXn2+DvH14HF+Y1vcE1dD1KLx/JmefgWAf6Rp8Sd1xJPWh4vrVUkN7KxsK0L7BfRjSzXILdQ1Jy1THZnu2pQD9DnWXN4OJwIzE5b9CMfmAuD7d8BdnyztTekN+DZ9BN1CHwLEwswuNUF+2sf9iXJPtf2m2PgWG3iz/G1YBb7x1wf1fBI9rmx7YVhigP6LN7GS2bL3b/BN8odyWte1PqQ4DMzq39VzT+p2oTwYXgvv3gWlBNeRvCnKnsa3+7wj2IUtDMaUAnbZLZXpx28+3YPveMD38K1t2PI690I8DsvV6NhPMA8OzttMpi+kLGj3GYsVW5tq8J93uR1gIJoB7wLYUgL12d0C677bN1t9HqTaC1Mci1D2HFGA8R/UUuI3n2gXsy2VJ92u569Gc+/eirP+FKZNGU/GY9teYa9Gd7X4FA6UvmSkgPYdDqBfK7YeCx9DbYir3jLj9zuD+JhLGgp7gS+9VmBgauh6F98/x7GN/V4PnvFW2fDdlUjrmvqmhMeWEjdmoTrbxhr0LvPlLyYcjaToqc8NPMBlsCZ+D8gIVkwFZ7QBeCPcxI/oP/A5DwH7t5wnoDH3BbKwp8sE8tIAbsg4+zEofaDMkXy63w2MwBtI5mBl7XNtKSc++hKHgg6V8OHxYXnQBzTCuGJvlGDhWgvNhmax9UsofwBtcjQSDV0P+ei1ScDuHutfhGjC4Nka+bK7I8Wa207KUPuRvg+P0hXkMOEavtw9joRbPGjz3T+Ad8GVUTgYVNcm4osF/X2aL1+EPeCHbesqs9NptASPgBNgPlN7m9RILBhHvtSezFfahl+kcLqD+M+hluheoNng97HNItl1j79+F2F55/xRTY67FIuzYEbzfPCefp4uhmOagcQDMCufCXVBM5Z4Rt/dFoG6GYfA+TA2OxeegsdeDTcfKZ0J53xsTDPJeg3wcGc6yWmBc0bh/Naa9yJu7mDrkGjvl6pNl9QkpDYZJPtjp4UttqbyByhmwFTwNXeFCMDtSG8KxsDF4E8g/wQvvQ99YPc6Gp5XY+BTafch3gV4Ze1MeBqdCU5R/gA2Q3lgGNDVqXPHnv/+gdiIY6B4Ez2dG8MEupob87cJO4gvEG1tZ96Wirw3J8fhgFWq6rCG/zofFfueGabL1+cKArj4dV4z996NcvVj1XRo9/3nBwJJk0DoAboO7UyNl8tWm37L29PyZ4d4HU8FTYBC2rdDb5BOr/vwRM+/fTuA94fYGOOXy15D8aOh6uE9T71/79jjfunMRpWOXuxbp+RyZ2/+rXD1VDc79weB8AewHpXQKK8o9I/qs9CfJhCapsdcjbZ+83YiGdbLGD7NyUkrjQzqnxtzb2a7jPob8uVDnFR/CvDRFTTGuGBt80g1jkwYOgQlgBZgPNHgPOBSKyT6vA9+2J2Ub+HFKeSF6wh3g+qXgdjAI7QytpUXpyLe+ZXc4CNTfxhV/Zszpxs+aixbexIUqDAr64w14HHwPPiBbgC+yvNJ+Keg05K8vhE/AcS6WdTQ95bxZvbnFc9mOq1D6KUfNDgZnXwBvQaEcq/KaJbl/Od2SrTQYT5nVDZYnwI6wN+RVzOu0/p9UfHFsBSvB9VBMKbAXrjO4DACPv0a20vOdJ6tbNHQ9mnP/+pLy/vB+L6bGXIvH2PF3WBEcg9p4XPHnv95/94L33rmgt+X89Nko94ykF4bPfZIv1JvBe7Cx1yPt+0hW8dobR7yfD4G/QGEcGkRbo+UFbS8qvGApgJhp7Qn3QOeCk/Wim314YXaHO+FxWBJK6eJsRXfKF+GVbNnju84AfQT0gG6ghoz9d9x0hy+SzbPlUsVurPCiF2Lw2g6uhWvAcc4I6qNxxZ9TDT6gp0Jan61uVuGYzQB82W0N+8COoNIL0OCtNoXDxtbGPVTl/L082+5+ylPg0Wy5JcUQdn4YfKhfhzPhKVBnw8ixtf/95woWDRJmZQZH74cVoZwuYuUL4DUx6LvPUFgf9KsPNFYfZxtuSOm94ZhV8nbcUvl/78pWX0epny9DPit09b1Q6no05v61j7zSy26mfGOuPoR6Q9fiG7Z5AiaBD6A/HAV57cVCetmsRX1ghv4XU0PPyGXs9Cu4nX55zX0peM98Bk29HvexjzoWDgDvJ+ONsSdpjqzSpACddq6H8ioG6U10YMFgd87aLyxo9+b2wrvPaPCmfyhbnpVSeUG8OD5QbmeWsT80pGfZwO13K9hwcZYfgB/A9fZrUEgvwquz9i0pi6kvje5XimVYNxG43bDcdv2pzw5qAhgA9vEbLA3rZcu3Uio9dP3JLmT6mtI2HxSlDy4bONVf4SOwbSgcn9VTn/OzbBB3vT6qhvydmG28rr/DGPBaiH047mIywLq+Z7GVWZseXQlp21HUT4TOoAr9sG1H8CVj3+9A8uh06qWkV+fBt+B+8gasDkkbULHdAJ50JBXbjs4aelE+CfrwExwHn4Pjnhw2Arf/DyQV9mF7H/AaDQFfNvbpfguAauh6NHT/juvlv/9uS9X+9/hv09hnzbYJs7aGroWbTQsmRwbr52F9sI+3Qb0ELhfi9Sqmhp4R9/G6GIjt035uhblANXQ9it0/e7Jfeoa8H7xW3t9J3ide37lTQ5TjHJiOoksDZnRifaksoIFdi662vx5gWUn5ovEBLqZZaMzfIMW2aWrbBOxQzicD4GyQHs7Uf0P+TsGG3dLGrVg6Dj3q0Mg+O7JdufMr1Y2+dIcpS23QyHbv1fQSaeQuf25msPTTTY8/W8a9aAxAheNq6Ho09v71/voUbswds1S11LXoyg5/BwNmSg4WpO64H4GWqtwzYt9e71KeN+d6zEx/3keFepeG9CmncF0shwPhQDt34CDOz6A2GE4Fg4HLL0MldTCdfwelEoXGHPs5NkoB+TTqb2bLhzdm5zrYZtHsfJavg7HGEMOBcKACDkxKn33hM/Cj9JdwN8wFlZTTJh5r7xYcZHH2fRwM9E7LfQSnQEOffNmkLuR1eaIuRhqDDAfCgYo7UDjFVOkDOiVS7GN9c47b1mNvzhibuo9TKK3lT1OPHduHA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDrRvB/zv4scnbcbJxv8NNj5d8TjXcKDlDnxOF/1a3k3TexifAvSm2ON/TX9V022KPcKBcGA8dmBfzn1beLWtPRifsknP9Wq4uK1NjuOFA+FAXTvQi9FX5X8br8pB6/pSxeDDgXAgHGgjByJAt5HRcZhwIBwIB5rqQATopjoW24cD4UA40EYOjE9z0G1kaRwmHKhZB6ZlZHtBVxgEfh/zO4Rq1IHIoGv0wsSwwoFWdmAS+rsNhsM10BuegkjSMKFWFQG6Vq9MjCscaF0Hjqa7c+FyeAW2hAFZSRGqRQciQNfiVYkxhQOt78CUdPlyQbdPsjxdQVss1pADEaBr6GLEUMKBCjrg1MaaBf3vyfKwgrZYrCEHYv6phi5GDCUcqKADTm8MhRnhEdgcTNBugVCNOhAZdI1emBhWONDKDoykv2nga1gVnoX1IFTDDkQGXcMXJ4YWDrSyA7/Q31mt3Gd0V0EHIoOuoLnRdTgQDoQDLXEgAnRL3It9w4FwIByooAMRoCtobnQdDoQD4UBLHIgA3RL3Yt9wIBwIByroQAToCpobXYcD4UA40BIHIkC3xL3YNxwIB8KBCjoQAbqC5kbX4UA4EA60xIEI0C1xL/YNB8KBcKCCDkSArqC50XU4EA6EAy1xIAJ0S9yLfcOBcCAcqKADEaAraG50HQ6EA+FASxyIAN0S92LfcCAcCAcq6EAE6AqaG12HA+FAONASByJAt8S92DccCAfCgQo6EAG6guZG1+FAOBAOtMSBCNAtcS/2DQfCgXCggg5EgK6gudF1OBAOhAMtcSACdEvci33DgXAgHKigAxGgK2hudB0OhAPhQEsciADdEvdi33AgHAgHKuhABOgKmhtdhwPhQDjQEgciQLfEvdg3HAgHwoEKOtDYAD0DY+hYwXFE1+FAOBAOhAMFDpQL0K47Cl6Hh2E1uAOmhVA4EA6EA+FAhR0oF6B359irwibZGPpRfgy2t4Um4iBTtMWB4hjhQDgQDtSiA+UC9IoM+HT4JBv4r5Rng0G7LbQpBzmzLQ4UxwgHwoFwoBYdKDevPIwBG6T75wa+IfURueXWqg6io2kKOuvMsuMzUDu1shM0pOXYYPESG61E+1cl1kVzONDeHNieE1oHfI4GwDkQqjMHygXosziXF2ANmBGegR6wOrS2DL594Vq4Kut8I8pl4TD4MWtrqPiWDZyGKaYuNM5UbEW0hQPtzIF/cj6bwV/gDzgcDoHTINSOHJiMc9kFjgWnNjpApeR889VwK5hNbwuXQWvJF871rdVZ9BMO1KgD0zGugZBPviZg+W6YFUJNd+AMdlms6bu1fI/8RSzW2w80Xl5sRQXavqPPv8IW4Eey52AMhMKBcKDxDkzOpi/Bb7ldzKKdmowv3XOm1EN1wjKDPJJ1gwv4kOU34RaYGSqhm+h0TTCL/rQSB4g+w4F27MBwzm0SWDl3jstQ97scn+dQHTlQLoPuz3lsB04NmM0uDQeBP1nhnPQ9sChUQt5k61ei4+gzHGjnDvzM+Tnn/B4cDKPBacreMApCdeRAuQC9A+dxHNyQnc9rlD/BErAv+CXe9PAZhMKBcKB2HPCnoqaG1bIhrUsZn0YzM+qpKBegv+JEehacTC+Wf8naOlH6dg6FA+FA7TnwNUO6ufaGFSNqigPlAvSVdPQArAIvgvNYc8GK4DTHEPCLvVA4EA6EA+FABRwo9yXhuxzPoOxb2O38GcrZ4X3wJzuc4giFA+FAOBAOVMiBchl0N47pF4QLgr+NlL60u5/6/hAKB8KBcCAcqKAD5QL0oRy3K/iF4A+5MTi3FQoHwoFwIByosAPlArQ/53wBPFbhMUT34UA4EA6EA0UcKDcHfRvbbw/TFdkvmsKBcCAcCAcq7EC5AO0fFvKPrfgrov5c5TsZZ1OGwoFwIBwIByrsQLkpDn9T0B+vK1TMQRc6EsvhQDgQDlTAgXIBehjHk0JNXNgQy+FAOBAOhAOt70C5AD0Nh7sI/OWUDuB0iP8N1bOwDYTCgXAgHAgHKuhAuTnoAziufxXrUhgOR4O/OXgShMKBcCAcCAcq7EC5AD0nxz4DrgR/5M4/MboT+BftQuFAOBAOhAMVdqBcgP6YY88G/pKKv0noX8fyC0LbQuFAOBAOhAMVdqDcHLR/b+MZ8G9v3AX+VIeBOv5CFiaEwoFwIByotAPlAvRbHHxuGAMG6j1hJPg/noTCgXAgHAgHKuxAuQDtofN/5Pu8Co8lug8HwoFwIBzIOVAuQE/OdufAYjBpbp/7qO+bW45qOBAOhAPhQAUcKBegD+V4s8Ih8GXu2N/k6lENB8KBcCAcqJAD5QJ0d47p34N+qELHjm7DgXAgHAgHyjhQ7sfs/Gt2W0K5bcp0HavCgXAgHAgHWuJAsQzan9iYMuu0J+WmMBz+yNoeoIz/USUzI4pwIBwIByrlQLEAvRcHK9aexhB/zS45EWU4EA6EAxV0oNj0xSsc7wXwT42aQfvX61xeEhaAIRAKB8KBcCAcqLADxQJ0OuQmVPyDSelnoQdQ3xp2gFA4EA6EA+FAhR0oF6DX4dhHwnvZGN6kNGBvli1HEQ6EA+FAOFBBB8oF6I847loFx16Z5e8K2mIxHAgHwoFwoAIOlPsysC/HewTWBf9I/0IwPZhZh8KBcCAcCAcq7EC5AO2fG10GVgf/V5XLwB/B+x1C4UA4EA6EAxV2oFyA9tDfwq0VHkN0Hw6EA01zYBY2Xxx+Aj/l+hcna0X+vfgToRv0hBPgOgg1w4Fyc9DN6C52CQfCgQo7sAL9Xw3zgL9E9iFMAbUgx+F3V/5J4vXBT+CO0WnSUDhQ1gH/rsj1ZbeIleFAbTswE8PzN3rNTJP8SavT00KVS/8zaf+4Wl6O9cp8Qx3W/a///Kueba7IoNvc8jhgONBsBxZlz6Pg/VwP/ifOS+SWq1mdiIMX/qbxaNr808WhZjjQscg+/Whz/qiUHmbFYaVWRns4EA40yYFJ2drpilHwdgN7us10Bdv4rNpHMU1LY1cYDs5XV1oDOMANYIwYmh3sYMoXsnoUTXSgWID2I5Pt/uSGb+sL4GmYD/w7Hf4qeCgcCAda7sAcdHEafAJzwy/gb/BaFpMB0N/k3RP6gsH3PPg3FMptNoZPYc2MNygrKTP7w+EdcEx+YTgY/gWhVnbAwLx9QZ/LsXxvQVulFztwgC6tcJCYg24FE6OLVnPAzNf55N65Hq+nboJUTj4LN0J/8FncAgpl2+XQKVuxLOVAmDFbrmTh+Hwh7AROybQH1eQc9Pc426PAXf9Y0pcFba2x6P/ccjX8AH486glJm1O5Ji1EGQ7UuAPLM75dwPt2wjJjXYp1x0L/3DY7Ul85t1ys+jONW0FvWBf8iYlC2X4y/JqteIbSBGW1bLlSxbR0fAWsBfpwB5jlh5rpQLkbyI9QO8NT4MUdAEdndYpW1QH0NgL8ssObyWP1glA4UE8OOCV4IBgYN4CXYSIoJrcpfP7MeDsX27iJbWbmhT8bPQFtUik57g/gETgIdoVL4VQoPE+aQq3hwHR0shto8u4wE1RCzln5Z02TzBAGw8xgvViWQPP/kWN9rARDaTfwh8KBSjiwMp0OA4Ns0slUSn2h7v1+K5htqg5gInSkCy3UNuzvF3P2qeYHg/YsLlRIHsOkrlA301CpuFF4rEotV22Ko9iXhAaxUj/FsQ/rHgK/mW1NvUVnZs9PZJ06x+ZFvR8uztoaU/jGlmLy5p++2IpoCwdawYGF6cPM8ddcXxdSPyW3nK+OZsGk50l4D8y0H4CzoaW6ng56wevwFJhoLQ7DoVJy6qVYpjwN7XlPKnX8dtlvsQBt8M1nAYUnXok56Is4iG9ag2i6oc+k7s9P2uZcVigcqGUHfmBwBsW8urNgIC6lr1jh9zoG0FHwLbSW+tDRZdAV/PTod0qV1Pt07vj/Br6YDNb3gVOXX0CoQg5MTb9mB5NUqP/U7aRUFkwLubI39V1yy82tGujNLELhQCUc8Pm4F3aEqWAFeBr8Anx8URdO1E8Bj4DB+SgolgTSXFeq2hRHOZecv7oCnLvyTe+XDreCF6EeFQG6Hq9afY3ZeeXz4DZwPrZYwkFzqM4cqFqALvd22wMTe8J88DbMAafDYXAchMKBcOB/HXA6w+9pWkNm5DvBlDAcrgKTpdB45IDzRKW0DCtOA4Oz+hAMzCu7EAoHwoGKOWDi9C5MAX5xvj48DuUSKlaH2psD5QK082crFpywyzHhX2BKLIYDrezAfvR3C5wMBuZN4TnYCmpVftl5DfhTXgNgXgi10IFyb+Sb6fsNMGN+EhaHRbJlilA4EA5UyAG/WCz8mWK/fFu6Qsdrabcz0IGxYnkwsVsILoJdwJ/uCDXTgcIM2uUrsr403eB8LdjuN9S+JV+FUDgQDlTOgS/permC7jdjuVY/ve7O2LYGg7Py56/9Ym17F0LNd6Awg/YnN7wRzodd4VH4D+TVlYXW/HnNfN9RDwfCgXF/QXI4RvjjeoNhIzBg7wPNlZ9+D4Gp4RcwqH4KrSF/yWZEQUdfszxpQVsstoIDJ9KHb2p/+8dvpX8swHmmelT8mF09XrXxd8wGvfvhPbgJTJbuAJOopmomdvgDlgJ/CW01eAymg6bIMe0Lp4M/4+yyWg8ehJTw+Yn7K/DF0h5UtR+zK2fe0axcodwGdbYuAnSdXbB2PNxlOLer4U44BYr95q5B1Cy0MyRdQcXMt6kyoPqTIHntwsKB+YYG6r4Y3oJjYFE4EkzgzPLV8TAQzNJvhH9Ce1FNBuj2Ym46jwjQyYkoq+nAYhy8HzjlYAa7P1wDKfukOlb/4N8Nsnoq/MmIy9NCE8oL2Xbhgu3/wrJJWGPllOclBRvvybIBO8ljrAMG8FqWL5V/Q394HjaEcqpagPajSCgcCAfazgGnEP8Or8LncHZWGjDzcmpx1nwDdb+4t72peo4dzGyTzNjN4G1vrKZhQ3+yK68nWXBOO+k1Kk7LvJIaarD0E4k/WSKrgJ8sdoZ1oeYUAbrmLkkMqJ07MAHnN6TgHAexPAW4rgf0ghtgUzAj7QJLwglwCjRVBmOf9cdhL7genPZw3rixGsaGaxZsvAPLtteTDMRXwrnwB3wGzqvvBDWnwo9VNTfAGFA40M4cMHPeHc7KzmsJSoNlHzCb7godYFVYEVy3P3wLA+E8+B2OAH/bsDFy+22gN0wJ/4KXoCn6DxvvBreB41wO7LMH1JMmZbAfFgz4a5Z9QdacigXofoyyW5mRPsy6w8qsj1XhQDhQ2oGTWWXG3B0mhk3gTjgI1JwwCnzGnN/dAsz0PoBL4AjoAZeCgftlaIx6sJFB32D9A5iR+3LwWAbfn6CcfmPlKrAD9IaRYKb/KzRHE7DTLrAeTA5vgeczBiqpZ+m8L/ii+SQ7kHPpjX3ZZbtUr1iWQ/vm3hl80xwMvi13BW+GraAeZcZyfT0OPMbc7hxwDtgA/AusDH5p5dzt8dAHki6gMh9sCdbzWpqFK/MNZeru70f5J7LSAG3A3xY85o+Qn0tmseLqwxEuBD8xmNX2AV9ebaG1OYgvAl92+noJdIBSOoMVfrlbU3Lg2xeMyEB9b0FbvSxGgK6XKzV+jNNMdI/sVP3Eejs4/ZB/vm5l+UgwW/4b5DUdC2aBDclM2Qx8F9gMDDZmwwNgBVAmX2ePrY2bB/e5N8s0Bsyftbdm4fm+Ch0LOr2O5XkL2iq1OBsdrw+9YUIop6oF6HID+54R9ygY9QIsf1nQFovhQDjQdAfMWmfNdhtJ+Tb4Cc9gqi4Hv9D6FL4Dg6XZdNKmVIakhTLlmqx7Ctx+VTBrfw8WghQMr6K+IBiwXoAj4GK4AgzcZpytqYnpbDD4osjrGxYmyzdUsD6Uvu+G/vA71J3mZsSa6MU1+xwAw2ERqEdFBl2PV639jrkLp3Yr7AR+xN8HDNqfgNmlgTT/rJ3P8mgwEz4ZbobO0JAOY4PXwcCrDNhm3p+DmbnqBT7rPuPOxe4If4AxYBZ4AAyqu8OhsBG0RE4nXA35T+grsTwK2ipAc6hGq2oZdOFHjPyIvVBLw4YwF1wL98AnEAoHwoHGOdCVzQw6X4BBN+lnKtuBmfGpYJA2o5wuK62vDgZrtTcYSLuBbcdAvj8W/4/moeUDmB18blcDs/GN4UuYEJYAg78B+ThYDq6Ej8Bg7AvBdbfAs+D3UEeDccGXS3M0hp3s+2PwJeCY1oPFwPnxUOZAuQDtJukt6w/Ie0G9aULhQDjQOAcMYH4R9xUsCU4VvAdJZsQ+U5PAHWCgMkC7nQFyeTCrfR7MYLvDJeA+x8KiMDk8B4eAgU9NBNfBFOCLwWX7PQ8mAzPnWWFZMEjfDPb5PvwdfAn0h8PBJG0BOBEuAnUvXAEbwF3QHH3KTo5lRegEfcEYE2qkA164o+B1eAPWAm+iaaEeFVMc9XjV6nfMGzF0P3UaHNXC8DBM70JOA6lvDyZDBsCFwB+7M2s1sL8Km8FNsAeop+F8mACcLjgZDgPl8lAwsBtgDbxnwmfwD1gEPJ79m7WrpcCx2p9z1a47EN6GR+E/MAfkZXDeP9/QjutncG6LVeP8DMKltDsrVoVNsg36UfqRxPZQOBAOjMto/YhuBlgopxGOgZ+yFa9R3gJOM+Tl+jXA6QqzSBMiA6YZpoHyVzC5MJBfDPPDCNgbXG/27fEPATPp+8G+VgID96ywJNwLs8OeMB/MAD9CZ7BPtzOD9Rk3OO8FN8K2YFthgFqPtm8hVCUHruO4fwFvADNoNRf4Rq1HRQZdj1etNsdslno8GHDNPL8AA15ePj+z5Ruo7wo7FrT9jeVRYMZrUP4ADLw3gVmv6gF3WUGLwzngc+n88XA4Gu6DqcDgbdbtGJP2ofIYuG9eJmgvw4mwGjwOHvMk6AlJnsezsCH0Ao83CPLHYLHdqmoZdMcylg5j3YrQP7eNF8gbIBQOjM8OOJdrYHTqQW0Pl8PmYLBVTmdcCQY+A+7ccCmYxaqpYUEwQDr3uzx8DTPD9+BUyN2g7NP91oSXYE54FNRkYDbdGQziXcC+jwEDqdoK5gEDt1oBTLaWgoFwJCj7PAk+AseUNJTKWmB/68O7sBCMgVCVHPAG8cK8CF/CMzAC/DhUj4oMuh6vWm2OuT/DKswej6NtpdxwJ6Bu5uXUhhnvHWBAVIuDUw5mrg/DD3AF+AXi7vAWGKy7g4HV5TfBYxjsEwZy903LBmiD7O/gs+tLRJzymBXUKdAfTgXnmB2fwT1pGSonpIUoxzrgdVysGl50LHPQUaxbGDaB2eBx8GOOgTsUDozPDvzEyRuA1bRg4O0F/SHJoHkQmO12gw9gJMwEJj0mOgZeM9EbYBDcCT5zz8GOsD2Ytb4A68B34IvB9hnhV5gI+sGKGR7DAD0LKMf5EHwE+8C+4IvBwC/D4AjoA8oM2ZdDqEYd8G3qRf83eFNYT/hRyTd/PSoy6Hq8auPGPDPF1XA7GNQMMtXUoRz8flgQDI7/gd/gG5gKysmEx8Cd1JXKPTAgNVAuAX5aPR4ehElhblgSDM4GZjNn9xsDvgw8voE5ZdOWD4B+Gdj7gJn0KXAUOG4Ds3257c5wDnwLHi/0XweqlkH/dwj/re1JNX+R83Uv9F7/3bSuahGg6+py/TnYyah5D26WtZgsXANbZsvVKMxiDcoGRBOWa2EG8MVxESg/aRqMN4ApIGldKn3SQlauR+k5bgt/hUGwMMwDD8PlcDuYNBmQDcajcnXHYbulAdfyJ/gK3oaP4Qrw+X0SvoB5QZldvwc/gsc6DwzSU0NonAM1FaAdUkc4EJbJ6i57U9azIkDX59XbiGGb8eXlx3sDZDW1Age/EOaETrmB3E99PngG9od/gMG3B6iJwWC7vgvIn8bw08HBsDvsArOB2/WFYWCAtY9V4GcwACdsT/VU2iZu+zoY0A3Wrjf4m627fj8waI+G5yFpbypXwISpYTwvqxagDbzF5AU9s9gK2rxxvKChcKAtHHDKzYwuL7PDfFaaX9cadedtd4W/gPe7c8ZHgwEuyWegC3yQGigd60xgsLsFnCq4Ct6Ak2EHcD+nE8y8DeDKAO12SfZjEH0KPOYvYGb8KCQ5RuWzavLkskFXpdK2Bce2jDsPXxq9wBecAdvA8xwMhzkg6XwqvkRmgE9SY5Rt78CEZQ45Deu8yV6DN8G5rA/hcgiFA23lgEFqW1ggd8DDqL+SW27t6il0uDvsDJuDge4oyOslFpwWOAQMqL4wDGxfw/dwFbwNZsCfgkHT4K26wftg4J0VtoL74G7wfA2Yi8M+YBI1W1ZSjFUKxvbpepfzLw83ss11Y8DjPAI7gc+80yNTg5n59NAbPNckt3GsvghDVXTAC1hKB7DCj1+XwnpwGXiTngShcKCtHDCI7AsGRAOg2Z+Jwj+hteRzsAP0AuvbgIHS+dh/wJFwG8wBHjvJaUATlofBqQL3WQvMoB8DtRi4/4pwLxicDX6bwwi4BwzAZrQGS7Ph37LSutsqg3E5GZCV21l3esNzsY9PwPaD4BzwJeKxLwHHbLBeG0zC3P5keBrcLlRFB8oF6DkZV/oItAd1s+kPwIvsmzgUDrSVAwY8P27PDd/DQOgJq4B6AIaNrTX9H4PZo/Ai+LHevpzWMBnxU+SNsAsYrJyyMGnZC+aCb8D5WsekVoIuMAZOBTPxj2FD+BLehcFgwnMzqF/B424AKTCbwebVUHBO6x3PlPA7OA5L1Rk83jUuoD1hN/BF5Ll6LguDPvwID8G5EKqyA4U3Qn443li+2X8AL/DU4E1qWygcaGsHDD7PgsF5ebgLvDd/gaGwODRHW7CTwfMg6AWXgsH6WHgSjgEDrRnvh/A6+LLoBI7jO5gZ1GewIPwLbL8SbgMD7wg4GDyG5+D5GMivhoXAgOpzZrD1pZGg+qdsy8tltxeDscFZ+VynNsspYFb4AHyROf4j4AkwcHue20FvWBf+DfYXqrID5TLoyxnbM/A++DDcA95A6c1PNRQOtLkDBhuD3srwTnZ0pz+Oh23B+dWmyGB1VbaDfRus7H8A+JF/IlgGloCtwPXLwWngHLMB1hdHDzBDdr/R0BsM9CeC692vLwyBzWBSMHjuAG6vUsAtVrre7fMy4M4Jbu9LIJU+12PgV7gVNoRp4EpYB9YAM+WjYQ94GorJ530/WCVbeQnlHVk9ijZwYMIyx/DmnBt8CA6H6+AsOBNC4UC1HJiJA98O3pdJb1L5BmZMDU0oR7BtCkD2sw30AAOwUyt+crwMBoEZqtnyPyGNwZfC17AJqHNgY9gQloU9oRP8He6GrpCCM9Wx88STURpcDcCWSfmAnK+73u0cS8p0DcrKQP0rnA4/gQH5cfCTgtnzU7AwDId/gIlYKZmYzQt6YiD3XLeElmpyOugD14PHnwFCzXCgA/ssBt5oiV7N6KcWdjmLQXhDhOrbAR/m+8DsLmliKoNg2tTQhNJ7vB8YjJaBB2AU7Aj/gpvAbZTByWDsNveD49gPvK8MYKXkM2OANXC+Az+BgVXMdA2qrk9t1vPk262nbd33hwxfNC5/BtfCRfApmCl/D7YlOYfu85BXVxZWg+XBF8rS4DnmZWB9JN/QjLp9Oy6nlHzZesynoSfUqs5gYMbBNlfHMkf0Qt0KXlxvqKSHqGhuKByohgM+3HeDmddJ4KdA79Oz4QtoqgxqBomdYAPoD+eCCYkB9Z/gNsptPN4C8A3MCWuAwfEKUIvDXOD6B0G9B/fCCtAdfO4MsvZlwO0Ayqw4L7exLd+e6j/TPhp8YcwIZuUe8zLYEwbD9HAD9IBHIelSKvnga3A8D54DX3K3gH34CSKv71lIY823N6W+AxvfCAY99Qn4gt03gyKUHCgXoLdno3/AlWnjNiod0+TgzRYKB4o5cCGNBkUfaoPYIWDQbqxWZkOD8i9wM7wLfSHpX1RmA6c8+oNB2mzZQL05PABOSxiAVwIzwbXBwO2ySYzPz2HwNxgG9r8u/A7KQOvYDXiWKrWl+thG/knr8+324xh8XgzCs8CqsDV4vAVB+QJxf8ecNBuVTtmCfQwC/eiXtQ2h3AgM/N3hI1COv1zMGLtRA/90Zf2Agm0GsrxFQVssNuDA4aw/oIFtWmt1Zzo6CbyxvPG8oX4EH5CdoDV0Fp1c3xodRR917cAujP416A1rgveagS1pWyoGf4OlMoA9Cc4rXwNmfwY0g5lB+xk4GQbDr2Cw+RScMvF+fjpb3o7S6Yh0f1smHIOkdWk535bWpX2+YnuP5z3tdq9CPxgMfWE0+IWeAdpseQgsAwvDo7A8KNt8Ngp1Lw0rgn3vDUeBL8GpoCVah51vKOhgf5aNN7Uqs/2am+JwUK+AF9ibLsmb+9q00ErlufQzA/iG/hAMzlPAfHA2TARmTQ3Jj0puW0xdaOxQbEW01aQDBo5DwCx3CXDe9GFoicx0j4Y+sBt4rxwLBohnYRSsDN7zBiZlUD0N5oGfoTd0A4PWNmDAMZEwcDlWg+NssCt8Cmaf58B5oH7P6DR2adw/Hiu9EErV0+ZuZx9Twm2wIfwGvcApCMezOlwBeqb6wXUZT1IeDU+B8px8CeXlczI9PAGzw3IwAnwW9aMl8mWxNrwOJ8IiMDdsDqECBzoWLOcXvQEngnfhp9wK39qtLTOZZcEbOulbKs/AfnAsNCZAb8p2m0Ax+ZFveLEV0VZzDvRgRF77eeEdMFhcA1/DS9Bc2c9HcCAYyAxsR8CsYEJggJaukNecLJhxfgwzgUHqb+A9uhqYkRroPoEvYAE4Hgw6rv8GOsMkYIDtAOkF4HJSvm5bftl62scA/RB4v6vL4Ho4Dm6Cq8ExJfWgYoAdDbdACs5Ux2bePt/7gi+RTnAqPABqSIb11pLP9BowF7wJx8AYCDXBgavYduMmbN+STe9m561LdHAC7deWWNeU5rPY2Js4VPsOHMwQC1+0K9J2YguHPj/7fwf5ADwjywbcXqAM3p/DU+AzsBkYGG3fDUwiXG+bGexzYNJi8LOfz7K60x9uY6BbPau7bHDNY1shrk9taVuXfaF8CWbKj8GTYCLluqlgYngVLoUbYUKYB/qDQVxGwuGQV2cW7oAnwPHaZ/7lwOJ4rZqc4riLS7Iu3AneJJXU0XRu8DwAPoDvwIdoXugIf4HQ+OPARJyqwS4vg1OnfEMz6mbHr4FBalWwP6dR+kEXMADvBGany0I32Aj2ggthNvBe/AY+hNnBIGhGbHZ/MbjfKrAOfAVrgdm3MpAqg591y2Iq1v47G9ruc2G/PcFnY3tYCDy+UxVuNxwmh4HgPgbeqcGMfk+4FhaBV0H9Ap5nqI4c8EIaKL2pB8E7GWdTVkI+lKvATuAbfldYGbzBWkNn0YkvgVDtO7AcQzSQTpEN1QzPgGawa4ns52b4OxwHR8DWYMCdH74F73uDsMfzHvS+T1MJVMf+XQ7XnQhm9cPgeTgfDHQGyo+yupm1nw4Nmu6Tx7aE7YXbpOVU/sw2PounQ0d4AYaCPvmMGLgNvB5zNVCrw1uwHXjuSb6U1k4LUTbowBlssViDW1VgAy90Kd3LipeKrPQGrIR+otPHKtFx9FnTDkzK6A4DHwAD1THwNPhQDITrYW4wu30QWiIDqC9/E46DwYBn2RsuhffhEjAoOi4Dnts4tltBXQ67wzzgdgZKg9+bYKA2M50exoDTEStBXp6jyiceqc32YnU/PXgMpS+Oy77Xh3dhJNifQfxTWA/6wSPwKrwHnrvymd8BtnUhFA7UigORQdfKlfjvOJwaMHj0gW6wABhYVgA1EywO3V1oRU1FX5vBluAxlHOvBr9DXch0EqUBU3xBLAEXg8lL0rJUTFpeAbNq56DNYg2WfcHA6P6/Z2XqL7Wl9nxpcE/L1t3WcjQ4zhMgBezu1F1/IMwPBuoBcAe8CAZy178EnvOD4KeHUOMd8L7wJV0Teo5RrAyHgw9PIf+mrR4VAbr2rtrGDOnygmEZZP5T0NYWiwY0XwYj4AIw4H0PTnnMBa5/GMz2fbHk9SQLBlA/BRqYzXg/hm+yegq2lqlu0Cysp/Wp/Jxt3M4sfV1wbAb8Z2FHUH+FPcbWxv1jgPZl4fSH4xoCb8AwsC/PK9Q0B6oWoDsWGefutA0BPyr1A2+2JN/andJClOFACx1wjtmEIC8DS9d8QxvVT+c4ZpxrwTywA9wEj8L5cAsYHK8AM9RnIGkaKj4vR8BCYKBfEQyWBkWVytRmKSqtG7c0Lth/xcJE8DR8ASZGk2bLBuw1wfH4QpgMkvz08QnMDN1hO3Ce2nNx3GbRvhS/hVAdOmAQ9sbwhlgnq7ssW0H+4x2LdaPIoGvvUi3MkO6HfIDZjeWLqjTUOTju6XAeeK/7LLwPjjPJbcymp04NlAZEM2afkTXAKQ6DbsqELRtadpshYEJk0HXu2wDsNMybYJ/3wD/AvjaCB2A6MLP3xdINNoBX4F24HpIMzgeB5Z0ZBuq89yyGijhQtQy6yFjGfqTzBijGd7TvVWynOmiLAF2bF2knhmUwMjCfBLeAga5SMukwCB8Pc5c5yDysc6plBBgYZ4GkM6ksATPBIWCANpDLe2CwHZOVKTCnMh+084HcoGyQt3Tbt2E0DAf7MhibcTuWfeAueAqWBQP09xnu9yYYnN+CLrAF9IMfQa83B4O55UPgJ5lQaQdqKkA7TKc+nKtaJqu7XDjvRlNdKQJ07V6uRRjajrAZdIZKyXvaKRWD2qpgIFwNCjUtDa5bGx4FtzHzNEgbJJ0yWBEegx3AQG6wNYFxauI3cH9Lg671FJjzddvMktM6S4Ox2/wEBmf7lZfhKjAYK8dlgHbbXqBv88NQMIDbbt/OY78Oz8PXcAPkZVZdr0lX/jwqWa+5AJ1OdmYqnWBy2A82hnpVBOh6vXKtM+4edPMBmFEmdadiRlqoY2nYNmvcntIAbHA38+4H3ktPwCagJoE3wAD4KRhkDbYGRQOldct8PbVZpqz5KOpXg/sb7A3uK4J9vgjufyQ4/WjGfAw4BZLXUiycBrOCmfNt8AzcBb5k5oK8nCo5NN8Q9f/jQNUC9IT/Zyj/bViO6nswPTjAHeB42AVC4UC9OeB9fBMYDJM+ouLH/kKZkMwAB4CZ6TWwOphR3w7TwOywd7b8K+VCcDd8AwZXg6lB8gdQLksKyBNky8Mo/YRqtuvz5aeIz8H9voUn4HBYHNYHpygMqgeC682087If+7Nfz3cT8BPDBnALmF0nud258EpqiLJ+HLiAoW4NZtAjYQHwo+j9UI8y67m+HgceY24VB5yeMJucLNebQfat3HKqPkJlCGwI58NgcOrg72CAvxF8DhYD+zwalD9d8S44tZECsqVBNS0bQL8GA7WB/E1wnX2a6Xpsg+5nMBregAXhSTglYzVKZfvb4LkpPx3Y1xouFJHPsln0rbAzGMD7QGO0AhvtAHriy2V80hmcrNe6pnQdo1kdVgUzabU0+HGvHhUBuh6vWuuOeUu6M4CtCWaU/WBJyGtdFl6Ci+B52AteBds+gI/hMbAfMWA9De/DVzAIDIJuOwYMxL9ldYPunWDWbmB+CIbCK+D+34MB3O3t2+fOutsvAcX0FxqdAjkH7oFdoZwMro75r7BCuQ1z646lbjDfEW4Fs3o/WYwvOoMTrbkAvQ2D8ibyDX0EzAfvwPZQj4oAXY9XrfXHbFA6Bg6DuYt077SGUwhqedgfnBow2Jo9G4BHwKMwBD4C1w0Gs+keYCbtdpeBn9rMiM2iB8IQ+AEeg9EwEgwAv2bYZlAWXwYeyyDvM+iYbwaTp16QNBWV+WHm1NCK5Zr05UukY65Pn6X9csvtvVqTAVrTvThbwITQE1aBelUE6OJXrgPNPmxmdvdBfo6SxZrSLIzGzNbgaBZXiaxmJ/o9FtTy4LFMVMxo5RcwYBtYDcwGXoNrWn8o9WfhDlgLXP8SmB2nbcyI0z5+IjUoPw0Ggm/AoPwKOJYVwADpy8H1BuGlYEBWUlRUB9P7RgVHMBZcVdDWnhdrNkAn02egkn+DpvZ6KiNAF79aV9B8L3SDaeB82BdqTZMzIAPcJjABzAP9YRFoTU1EZ3eBL6pBcAkYiD8Hn4P7wCDtWMxsfbH1BV8Yb8FXYNB9AB4GA7FzzmbNBmL3+w7eB4Oc+zwHBmX3c/uRsCP8C54EXxC2LQ1JvpyuTwsF5UYsXwCO3amcluhv7GyQzms1FrxPxhfVZIA2az4KXge/qDAbuAOmhXpUBOj/e9XMhJ4paO7E8oMwdUF7tRcNmAcVDGJFls8taGuNxQ50ciQMgdfA+94gaXB12WBqkJYp4Fo4GwbCU2AWPBI+AAP4q/AZDIMP4UsYDvpsn2bcnkcPOAZ8IfwGBudHIB13cepJk1G5Ny3kyv2pO1ZfYPOCY9sdmivPz3FuCVPCSuA9MyOML6pagDYIl5IXdVXYJNugH6VzYi252FlXUdSIAz7kfpTOywzODG/SfGMN1CdmDAa3vL5gwXNobRkgT4QXwYDsvT8LqLnBoOWL7A34EfaDXaEnGBiHQFcwY14EfobBYPCeCoaD2fCn0Ae6gcsGAs/JAG9Q9nufReElWBB8ISVtSMVt8zLD3ws85jvwNuwMm8J00Bx5L7j/mtAX7M8Y4EsqVGEHOpbp35vhdPgk28YH92y4CLx5Q/XvwCBOwYd6OXg6O50tKJeEWnsAzTKPhfvhJ1B7QBr32IZW/udU+nsOVgADsoH7B5gavofFwKmPc8AX2mg4H16B28BnZwFYH66Fb8HA/R74wjH4/wtuB481Pbh+LbgcpoXn4UuwP7cxmPsses28Vnm57+Ng9p1klm/m7n6OtTnynHdpzo6xT+Uc8MY5GSaBN7LDHEzpjVaPOotBX1+PA6/wmOehf7M1r+1h8AzMDLWo/RjUcDgAroB/Q6W1PAf4Ggy+ZsEGq9thOzAzNlM1kzXT7A1JU1K5E14Dx+y+Bk/7MmieAgbad8F+DoK8fM58AcwIE2YrjqS8BXYAn8tCTU6Dx5w3t8IXxKfgCyTUPAfOYDevRU3Jt/lQeBF8g/vgjoD5oR4VAbr0VZuGVRvAutC19GY1scaAsyGsVOHR6MmOYHLiy+AbuAluhpGwLzwGC4NB9DYo1PE06Ovv4CfPP8BnylIeBjP0PmC/ZtpJbu/LIa9/srBHvqFI3U8/9r19xtuUS0Go+Q7UZID2dCYDP9r40XJV6AD1qrMYeGTQ9Xr12m7cTgVcDF+BUxEfgQFvIJioWDfoWZ8d1MwwANI8dWfqq8NgcCrC6YwVYR04FAzMr8B5sBGonnDF2Nq4f8zY+oOfcHzuzNjHQLHMmeb/UXeWjoIjIY3xfzaIhSY5ULUA3bGBYfpx7vIGtonV4UB7ccDA+iF43xs4nRowGL8In8DZYMDdGgzgJi0G0DnAOeoPwEC6AywIQ8AvBheFJ0DdD5+DQddny0B+NUwHs8IE8Ae8DAeALwunRJxK8UUwCsppclYanM3aJ4J+sAh8C6FwoGYdiAy6Zi9NzQzMwGu2dF9uRE5xGDDvhQPhafgRzKzNml1n4L4dXgADo9MWW4JaGB6EqV3IZLZs0N8Z7oFe0BceAvtqrgzu70H+S/yjWb4EJoRQ8xzwnqjKHHRctOZdsNirfTowKaf1FpgNbwx3wyzwMywDPqQLZcsXUBqcr4Q1we2d/zU7tp9HQb0GBuz8A/4+y4eAgdNAfwIMAfsxA14BmiMz8FfBqY2k46h0hh6pIcr6caChKY76OZMYaTjQOAd2YbNtwGzT+39d+B6UGfCZYHB1amF7MAivBs79Lg7uZxA8BXrDdPABJP1EZQAYZO/IGueifCirp+IxKs/CteD0gy8G9SEYpJujMez0W5Edu9H2a5H2aKpxB8oF6IMYuzdoOS3HylHlNoh14UANObA3YzGTNSgbSLeAy+BCMGg7ZzsM+sAt4LRYFzCILgU+L84Frw7qcdgPnLNW24HzxGbb/cHt/wXzwVNQqIE0zAn2r6YDs995XWiGPmafT6EfmMlPBrOALyDPK1RnDngDlZI31M5wOTwBC8Df4VYwQ1B+9AuFA/XiwFYMNAVnx3wTbASnwe5gJrsvDAXv+WvAzNfMdFuYErqC2fXJcB8cCrZdCZ3AZ+ZqMMHZE5zeWBbMxAt1DA0joCcYQN3X43wEzdUv7JgSJ+tD4HUItTMHLuB8Nis4JzMDb8p6lNnQ9fU48BhzqznwKD0ZRPMy4zQo5vU8C3fnGgzABtg0j7wHdYPwRbAGnA3vwA0Z3Skbq0nY0OPvBvM1dqcS2xmYvwPnyacBM/p74SrwJRFqngNnsFu69s3roZl7lcug/dg2e0G/DvLHgrZYDAdq3YFtGKDZ6UxgQF4SnOtdBKaHhyAvg9paYPsjsDI4DXEcTABXwF+zuoHVDPgr6AA7wGhorEax4XWN3biB7czaB8Lr2XaO2SmOJWCKrC2KOnKgXIC+jPN4EDaBF8CLPBusDaFwoF4cWJOBeg9vCj+D9/LT8B8wsN4FS4FBWTmNsR8sBrOCX7Cle96M1IzXT2OdoRcY4LcG+zZAHg8HQ2trajp0mvFHeDHrvAulY/CF8D5MCINgR7gS1C+wIezvQqi+HCgXoN/jVJaG9cAb4yjoD79DKByoFwecPjgCnF9WC4I/gfEbbA+TwxAwuH4CBrK9YUiGiYlBzm3VSNgF7gSDpuvtS50MTid4jDegteTL4lgwO7Zvn9s94V/geHrAT/A0+AXjZeDLx3XOr+8BwyFUZw6UC9C+jXeHLcCPdWbTt4Ef6WzpwugAABL2SURBVL6AUDhQDw50YZAzgsHq82zAb1F2gq8zulFuDtPB3+AVSDKAv5QWsvIHyjEwBFJwpjpWHmeirN4ahWP3+POD41bXglM0fjK4HZSfCCYGs/77wXnw9WEfuARC7cwB39D9wDex2YA39PlwJNSj/Fh6fT0OPMbcbAcWYc+XwezxXTgQDFx/QC9ojAza94LPQdLaVJ4FM1WTmKSVqPwKk6WGVigNwmb3ea3FggE6L+e/H8oaelNuAHNmy1G0zIEz2N1PMW2ujmWOuCLrTodPsm288c4Gv7k+MWuLIhyoVQdmZmBmwouCn/qWAac6tgTv7fegMTLrPgEGwb4wOSwLZqdqMBjszZxXBacDzbBbS7/QUWHAL/bcGqBT5t6/tQ4e/VTXAacxSmkYK7yR89qQhRH5hqiHAzXqwJqMa1d4FZxT3g72AjPpJ6EpeoaNZ4ev4G3YGpzmE+ehH4AXYTOw/9aUn2J9yayXdeqXlB7fF8LfszaLY+C53HJU24EDxd7E6bScEngB1oAZwZu0B6wOoXCg1h1wGiM/P2xgHQNODzRHQ9hJCvUzDY8UNrbi8ij62hkeBDN4k6prYT/w+ewNZs5Pw9EQakcOlAvQn3Ge84EfCWeDxzO8yUPhQK070J8B3gdm0K+B9/rhcD/Um75mwEsWGfRctDmVMxrM5kPtzIFyAboT5+pc2uUwD3SHzuDN0JZybs1xmqmEwoHGOjCEDXeAh8EfiTOQ+WXfTdBeZLI0tL2cTJxH4xzwRnY+7dRs810pfwS/cBkEzrm1tmalw6vBF4IPVE9I2opKazxUTtlcnzqNcrxxwC/Y5oUZx5szjhNtbQfOoMOq/BSH81mFMmN+Hk6CSeBs2BT8osKPjIdBa+sAOhwBS4Bz3QOgF4TCgZY64Ev/bfD+qqb8sTd/DO52eB3mhlA40CQH/IWU72D6bK+1KIdndYvVwW+VW1vv0KE/ZJ9k1jwYzOYjg06uRFnMgSlpnAOmKLayRtpScpOeq4UZl8+R93eo9h2oWgbdscCbNN87Oms3IDvlkNSNigG8tfUWHZo9P5F1fCPlTHA/XJy1Naawj4VKbDgf7b+WWBfN9enAZgx7FzCJWBN2hkeh1rQTAzoG/OJdvQZXg1n1hRAKB4o6MGFBqz+W5M2zK8wCZq+3gjK79mFwLrq1dREd3gyH5To+k7rHdu64sfJLE4NwMX6n/Y/GdhTb1bwD/lSD98w24C+i9IZLYX6oNTlV6Pc4eTn10iXfEPVwoDEOLMBGQ8Fg1xcMzGazb4PBeVKohOx3wSId96bNF0NLZaCPLwlb6mLt7O93JMsXDMdf5jikoK0WFn2JPAwpIfKLdpOFqnzxxHFDTXOgZqY4HPabMBs4p5emMyz3h8fgF6iEzDDeyHV8CfWDoX8GRSgc+NMBEweTiLwKl/Prqlk3MVgc/Jnsu8AkyJfJyxAKB1rsgD8QvzZM3uKeGt/BT2w6XeM3b3DLyKAbtKiuNliV0frdxWTZqGegNCst9iks26Tqhd+DmPXPXPWRxACa4kBNZdAGxSvAHwM6FTqA88Ffw6+wHHwCoXCgmg74UxCngZ+67gU/9a2TLVPUpHyhhMKBRjtQ+FMc7ngO+G3zKeA83xzQEz6GG8D54OOh0rqKA5hFh8KBUg6YSNwDU8Hn8A2EwoF240CxAO1c2Qbgl4KXws5gcFb3wepja5X/Z4/KHyKO0A4c+IJzkFA40O4cKBag/cH/lIl8S90vCJOc4pg4LUQZDrSiA+vS1wowGvyxSzPiUDgwXjuQfuyn0ATbDd6WfltuXZyPDoUDre3AcXT4d3gABoFTbH6hFgoHxmsHimXQGpKmNJI5Zs5JN6dKlOFAKziwCH1sBguBvyilhsE/4K8uhMKB8dWBYgHaB6ZUZq1Po8ZXs+K8K+LATPR6OaTg7EGehKOthMKB8dmBYgHaOeYPodgP/XenfV54AkLhQGs44HTG9mBS8HvW4cKU8WvQmRlRjL8OFMuUn8GOaTJLlqe8LGfPstSdKwyFA63lwEt09C6MAH/G3p8gOhf2hFA4MF47UCyDzhvib2nNlm+IejhQAQf60OfrsCn4K/87w/sQCgfGawcaCtDjtTlx8m3qwG0cTULhQDiQOVBsiiPMCQfCgXAgHKgBB0pl0M5B+4dnukFn8O9zqK7jivg3HAgHwoFwoNIOlArQbxYc2G/ak+LnoJMTUYYD4UA4UEEHigVof5TO3x4spfzPq5baJtrDgXAgHAgHWuhAsQDtt+ihcCAcCAfCgSo7EF8SVvkCxOHDgXAgHCjlQAToUs5EezgQDoQDVXYgAnSVL0AcPhwIB8KBUg5EgC7lTLSHA+FAOFBlByJAV/kCxOHDgXAgHCjlQAToUs5EezgQDoQDVXYgAnSVL0AcPhwIB8KBUg5EgC7lTLSHA+FAOFBlByJAV/kCxOHDgXAgHCjlQAToUs5EezgQDoQDVXYgAnSVL0AcPhwIB8KBUg5EgC7lTLSHA+FAOFBlByJAV/kCxOHDgXAgHCjlQAToUs5EezgQDoQDVXYgAnSVL0AcPhwIB8KBUg5EgC7lTLSHA+FAOFBlByJAV/kCxOHDgXAgHCjlQAToUs5EezgQDoQDVXYgAnSVL0AcPhwIB8KBUg5EgC7lTLSHA+FAOFBlB4r9p7FVHtKfh5+IWmf47s+Wyldm4xDrgMd9FN6CUDgQDoQDVXGgljPoTXHkzDZ0ZUGOdQX8AJ/CQFgVQuFAOBAOVMWBWsmgB3H20xQ4YBbr+AzUd8BOUCl1oePbYV14NzvILJQXwOvwZdYWRft2YAJObx9YH7wnfoX14CcIhQNt7kCtZNAG3y/gLFgk4whKg6bLh0AlNSWdvwgpOHusj+F96O5CaLxw4HjOck7wRb0y/AfOgQ4QCgfa3IFayaCf5MyXgPPAaY09wKzV6YaPoLHakQ03K7HxfLQPLbHO40wF3WBkto0PpfPRjik0fjiwBqe5PPyWne6llAvDYvBC1hZFONBmDtRKgPaEv4O/whYwAJ6DMdAUXcfGN5fYwcA9eYl1Bui+GcdR/gyXgxnUYAiNHw58y2mm4JzO2OkNv7AOhQPhQOaA8793wwmt6MiW9LVnA/2ZPZ0LF4MvitD45cDZnK7THEmrUfkDSr3Y03ZRtm8HzuD0/BTV5qqlDDp/8sNZGAGn5hvboP4Ux5DQ+OnAUZz2y2CC4E/yLAhzwvcQCgfa3IFa+ZKw2Ik73REfLYs5E22VcsCprnngQrgftoYPIRQOVMWBWs2gq2JGHDQcwIHf4flwIhyoBQdqOYO+CoPi509r4S6JMYQD4UBVHKjlDNoftQuFA+FAODDeOlDLGfR4e1HixMOBcCAc0AF/tXV8kb+ReC+80owTnpp9/Ea/PUy5+FL215hHQ72rPZ2Lz6JfireH6+K5TAyjoN6VrksvTuTjej+Z9jp+f5vMn5FtD5qXk/BvjLQH+SNwl7WHE+EcZoWr28m5TM953NhOzsU/A3Fbtc4lpjiq5XwcNxwIB8KBBhyIAN2AQbE6HAgHwoFqORABulrOx3HDgXAgHGjAgQjQDRgUq8OBcCAcqJYDEaCr5XwcNxwIB8KBBhyIAN2AQbE6HAgHwoFwoLYd8L/f8meh24M6cRLTtIcT4Rz8Tdj2ci4dOJdp28l18VymayfnYhLrjw2GwoFwIBwIB8KBcCAcCAfCgXAgHAgHwoFwIBwIB8KBcCAcCAfCgXAgHAgHwoFwIBwIB8KBcCAcCAfCgXAgHAgHwoFwIBwIB8KBcCAcCAfCgYo4MD79JwcVMTA6DQfCgXCgtR04kA7fhjfhPPC3pepdl3ACF9fxSfhbkafBixknU/pbn/Wmwxnw6zAYrNer2sv1KPS/3p+TwvNpd8urckYvw2TgrxdfB9tCPWtdBv8V1HOA3pXx+79dGBjkTrCtnrQ5g30SusIM8CqsA/Wo9nA9Cn2v6nMSfyyp8HIUX96NZv+bqB9hEjA4G6TrVf5dkaPg3Ho9gWzcr1EeAr9mvEW5PNST1maw18K38CncABtDPao9XI+871V/TiJA5y9H6fpsrJoFPgGzTrO2evwozbDH6kL+7QM/jF2q339eYOgfZMOflHIbuCdbrpfCe2tEbrAG6Xr94zzt4XrkLsX/q/pzEgE6fzlK1/3LXFvD6jAP9IS1oB5lEBsND9bj4EuM2ZfljWCAuLXENrXabJbmJ7OkUVR82dSz6vl6JN9r4jmJAJ0ux/+WX7P4S8balCPhGhgIZmxXweZQD8qfi2P2fyfvB+uB/8N3d1gW6kH5c/G6KIOBn2j80taHqt70JQOeIjdo65/kluutWu/XQ799adbEc9Kx3q5+G413FY6TfkrjferD4PvcsQ3e9ZLl5M9lDOP2fPbIzmUmyolge3gma6vlIn8unof3r5mz18p5W69LvWk4A/YlmdSDivdbPao9XA997wb1/JzU473TojFvxt5+ATIx+JMcT8E+UO86mBOo55/i2J/xe11mhKkyvD71JD8JeA6+LHvAIFgC6lHt4XoU873en5Ni59Su2szQzgG/zPHjpz/B0R6mh+r9xhvCdfijgHtZridNwGD7wjfg/dUH6lVDGHi9X49i3tf7c1LsnNplmxn05O3yzOKkqu2Ac89dqj2IOH44EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPt34HOnOJ07f804wzDgXAgHGi+A9uz63cZP1P+klvenXoltC6d+h+vvgAdyxzA//R37iLrV6Tt5SLtpZqWZMU7pVZGezjQGg60h/+ZujV8iD5a14Fr6M7/BFXOh0tyy9YroQ3p9HQwcP5WiQNEn+FAWzsQAbqtHY/j6cBh8E8YDv+GI+BvkHQUlT2yhWkpb4OR8BqsBIX6Bw1bwL5wCnhf7w2vwMdwDBS71zeh/Q0YAhtDKXnMp8Ds+wKYCFQnOBe+hhchZeaTUL8IPLbrbob0v8E/Sn1HeB9GwF6Q5HheBX05FB4BNQHoie32eSTYFgoHwoFwoEUOnMne5xX0cAbLX8JGsBScDQbsJIOeQVvdDVfADLATfACFmpiG68CgZnDcB94E+14O3oNdQaUpjjmopzEsSN3g/zIUyuN+AR57anA8+4OZ+h9gsHSbq+F2UIfDwzAdGLTfhXT8IdQfgnlgPXAKqCvMCR7HF8W88AIMBvVXeAcWBY/ruS0NoXbuQLGsop2fcpxejThgoLsDni8znqlY9xc4FUbBrWCAXQjyGs3CL+A2sg1cBvb9NFwJ20Nea7AwEBzDG+BLoJjMnj8F138FZvr9Qf0AJ0Ja39NGdD0YVD+Hn2AQGMSTTqNiwL0HvoWZYB0w4zfIvw0XQtIOVDy+LyeDfV9YH0Lt3IEI0O38Atfw6flxvSHNwgZ/QD8wMIlBcDkop+6sfCa3gXWDYF7281Ku4blcPV81AzbQJznuV7MFXxZJI6mkqY8x1P8Nn8G9MBd0gCTbk36k0gnM6PMZvBl00sxUDoHkgXWz6VA7d6BjOz+/OL3adcAglvQ7lS5pgXJaGAFmmWaYC4LTEcp1tpWT284PKei6/4cFOwxled1c2+y5er76NQtr5Rpmpb4EGKh9eRTTRTS63uOaRfulaX7OuNh+vgR2hKSFU4XSYD0ALs3aJqPMB/ysOYr25kBk0O3titbn+ZhROqdqEJsReoP6BR6FvcF71WmCt8D523J6kJVbg3O7TpNsDk9BXo+zsAz0AjPfLaCYHqNxMZg3W2n2mg+eWfP/FM5VPwkGZz8FrAFmyeX0ECuXhVXBfXaBpDup7ARTgh5dCwdAqJ07EAG6nV/gOjk9A45B6WNwOuN+SDqZyrYwGJ6G0+F1KKeTWPkzuI9Z+EdgW172cRTY5/vwExSTL4QjwCzWueS54Fwop9NYeQI8C7eC88ruV05m6vtDX3gGPJYvKHUffApDwHaz51MgFA6EA+FAmzng9EUpuc7ssSmago0nbmAHM9vJG9jG1R3BjLyxcqzTNHZjtpsdVsttb0bfP7dsddKMguZYDAfCgXAgHKikA5PQ+TA4Bg4F58zXg9B47EBTM5Lx2Ko49XCg4g44x74mmPk79z0QQuFAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA+FAOBAOhAPhQDgQDoQD4UA4EA6EA/XvwP8HifOxITpyqfoAAAAASUVORK5CYII=" alt="plot of chunk ebseq"/> </p>
<pre><code class="r">cor(log2(true_fc), log2(fold_changes$PostFC))
</code></pre>
<pre><code>## [1] 0.5536229
</code></pre>
<p>So the estimated fold changes are correlated with the fold changes we observed in the data from which we generated the count matrix.</p>
<p>Finally we show that you can make an ROC curve based on this simulation, where true differential expression status is known. Recall that earlier we defined “true” differential expression as transcripts with a fold change of more than 1.5 between populations, in either direction.</p>
<pre><code class="r">reallyde = sim_info[sim_info$isDE,]$transcript_id
notde = sim_info[!sim_info$isDE,]$transcript_id
ppde = IsoEBOut$PPDE
sens = spec = NULL
qaxis = rev(seq(0,1,by=0.01))
for(i in seq_along(qaxis)){
sens[i] = sum(reallyde %in% names(ppde[ppde>qaxis[i]])) / length(reallyde)
spec[i] = sum(notde %in% c(names(ppde[ppde<=qaxis[i]]), setdiff(notde, names(ppde)))) / length(notde)
}
sens[i+1] = 1
spec[i+1] = 0
plot(1-spec, sens, xlab='False Positive Rate', ylab='True Positive Rate', main='ROC Curve', xlim=c(0,1), ylim=c(0,1), type='l', lwd=2, col='purple')
</code></pre>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAWgAAAFoCAYAAAB65WHVAAAEJGlDQ1BJQ0MgUHJvZmlsZQAAOBGFVd9v21QUPolvUqQWPyBYR4eKxa9VU1u5GxqtxgZJk6XtShal6dgqJOQ6N4mpGwfb6baqT3uBNwb8AUDZAw9IPCENBmJ72fbAtElThyqqSUh76MQPISbtBVXhu3ZiJ1PEXPX6yznfOec7517bRD1fabWaGVWIlquunc8klZOnFpSeTYrSs9RLA9Sr6U4tkcvNEi7BFffO6+EdigjL7ZHu/k72I796i9zRiSJPwG4VHX0Z+AxRzNRrtksUvwf7+Gm3BtzzHPDTNgQCqwKXfZwSeNHHJz1OIT8JjtAq6xWtCLwGPLzYZi+3YV8DGMiT4VVuG7oiZpGzrZJhcs/hL49xtzH/Dy6bdfTsXYNY+5yluWO4D4neK/ZUvok/17X0HPBLsF+vuUlhfwX4j/rSfAJ4H1H0qZJ9dN7nR19frRTeBt4Fe9FwpwtN+2p1MXscGLHR9SXrmMgjONd1ZxKzpBeA71b4tNhj6JGoyFNp4GHgwUp9qplfmnFW5oTdy7NamcwCI49kv6fN5IAHgD+0rbyoBc3SOjczohbyS1drbq6pQdqumllRC/0ymTtej8gpbbuVwpQfyw66dqEZyxZKxtHpJn+tZnpnEdrYBbueF9qQn93S7HQGGHnYP7w6L+YGHNtd1FJitqPAR+hERCNOFi1i1alKO6RQnjKUxL1GNjwlMsiEhcPLYTEiT9ISbN15OY/jx4SMshe9LaJRpTvHr3C/ybFYP1PZAfwfYrPsMBtnE6SwN9ib7AhLwTrBDgUKcm06FSrTfSj187xPdVQWOk5Q8vxAfSiIUc7Z7xr6zY/+hpqwSyv0I0/QMTRb7RMgBxNodTfSPqdraz/sDjzKBrv4zu2+a2t0/HHzjd2Lbcc2sG7GtsL42K+xLfxtUgI7YHqKlqHK8HbCCXgjHT1cAdMlDetv4FnQ2lLasaOl6vmB0CMmwT/IPszSueHQqv6i/qluqF+oF9TfO2qEGTumJH0qfSv9KH0nfS/9TIp0Wboi/SRdlb6RLgU5u++9nyXYe69fYRPdil1o1WufNSdTTsp75BfllPy8/LI8G7AUuV8ek6fkvfDsCfbNDP0dvRh0CrNqTbV7LfEEGDQPJQadBtfGVMWEq3QWWdufk6ZSNsjG2PQjp3ZcnOWWing6noonSInvi0/Ex+IzAreevPhe+CawpgP1/pMTMDo64G0sTCXIM+KdOnFWRfQKdJvQzV1+Bt8OokmrdtY2yhVX2a+qrykJfMq4Ml3VR4cVzTQVz+UoNne4vcKLoyS+gyKO6EHe+75Fdt0Mbe5bRIf/wjvrVmhbqBN97RD1vxrahvBOfOYzoosH9bq94uejSOQGkVM6sN/7HelL4t10t9F4gPdVzydEOx83Gv+uNxo7XyL/FtFl8z9ZAHF4bBsrEwAAN0VJREFUeAHtnQm8HFWZvk/f7CsJZGdJ2EQWkSUw4igQBlDZRQKDDirKvhokCYuDF/8YQNYRBxhGwAUCEiAo+xoICA77juyEBAKELUAgJLn3/p/3pjp22uq+fW96qe56vx9vqupU9alznrp8/fVXp06FYDMBEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzCBfxDI/GPVayZQFQJDOctWeWdazPYH6En0ad4+bfZGX47UwvJp9DDSer71omBL9CWkuh5Cz6BSbAMO2gwNQs+hu1HcOSi2mYAJmEDjEfgGXWoroLcp/9e8Ln+V7Vdijv8bZV/MO3YTtuW8c+tvZfti1A0VsoHs+D3K/ZzWH0PrIpsJmIAJpIJA1kHPprdHR/opy/uQnOILKGursfIhUvmd6EfoQHQ/Upnq6I9kinrnIZVfjf4DnYnmI5UdgArZVHbomOfRIegnKOvob2XdZgImYAKpIJB10I/k9XYY23KS0srRvgujbTnJ3HScouEHo32/YCk7Hemz12kjxw5l/XUkZx1nW1Coz32EVsk5YA3W30I69/CofA+Wk9Ha0bYWByGVDUA9o/X9WOqXwLlofFR2BMtc24kNfW6TqLAPSx1zPpqIRiGbCZiACVSVQJyDlvP9LpKjfCqnNVknLGeWbz+gQMffHe24P9rePdoudTEh+txlJXzg6ujYXXKOfTEqW51l/2j9TZYfR+unsHwpWt+IpawJ6UujBY1Gcu7Kk6s/yptr+S7aGNlSTEB/KDYTqAUBRY6KWj9Bi9DlSCmKk5FMTnvD9rWlDi5aXbZ4OVr7QrTM5oqV9uiMdfVzxc4xkp0Xo22QfgVchGQ/WLoI27GUQ78NzUKT0Aboj0i/HvZFiuZ/iWwpJmAHneKLX+OuK8LUjb7PUHeklMea6GokUxSpCFPWd+liuX8VrcoUcco0EkQWd+zSPfH/6stB1tnPLf3U8qmXbJmWJ6GZaA66FOk830NKz2QdddZxb02ZTF9Q+yClOz5H45AtxQTsoFN88WvcdUXAO6L10PNoc5R1WKy2m8plWy5dLPfv2Ggre0x2uf5yRy1NE8hB7ppXnt18IVrJ/5z+39DnDkUrRcdkF3KyWeuRXclZLmBdvw6yJsc7HY1EaseeSPnt65Es+2WTzXErB/0Keh31QzYTMAETqAqBuBy0nLMi4Db0o5xWyFGpTGkLRddZUzrgPaR9+0WFh0fbz7IcHJXJyf4lKr8lKstfjKZAUXgr2iZn5w9ZV/0L0VAk+z3KPaciXf0CUJlSFnK0Wn8f5ZuiYe17LVpOYZm101nRvv2jgt4sd0Zjom0vTMAETKAqBOIctE58MpKT0rC4UUim6PRRpHJFpHK2NyI5VJXdjrKmqPZBpPK5aBpSekHbcrJfQYXsOHbouEVI9T8Qbass15Fmj3uG8kPQnagF6bhcB/0u23GmKF/H6stgrZwDtonKFTHrpuUV0faFLG0mYAImUDUChRy0nPHjSA7s6pzWKEo9C8npaZ/0Nvo5yk8vKI/8GyQnnz32Kda3Rx3ZdzngNZT93Aes/wxlUNb0QMu9SMcocj4b3RZtl+Kgj4mOzf1ioajd5PAVeatutf9PSH23mYAJmEBdEBhBK7PphmINllMdjQYXO6jAvlUol7NVeqSQDWNHr0I7V7B8VT6vm6Y2EzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABEzABE0gBAb2ePk22F531K+3TdMXdVxNYcQLvUMVdK15N52tIooOWAx2APuh8d4p+4jvs/Sn6fdGjvNMETMAElidwFJvfQ48vX5yerZ50dQqajVpRG1qAnkb7o3LYPlRySDkqch0mYAKpInAWvd2sFj1Oys/98+j8CLQzegXJOQ9EG6BzUW90AbKZgAmYQGoINCWkpzvSjoPRk+gTpAh6PnoAHY32QDYTMAETSBWBpDhopTLGFSC/C+XzCuxzsQmYgAk0LIGkpDhOgvBUNAG9jD5CK6H1kdq4E7KZgAmYQKoIJMVBPwb1TdFWaAxSPlpRs/LOM5FSHjYTMAETSBWBpDhoQV+IZqSKvjtrAiZQKwKZ08O8zT8Li59tDqM+rVUjOjpvkhx0obaux46+SFF2R3YQB3y3wEFrUf4iurDAfhebgAk0OIEzwtsbh5DBRzR9LxPCan1Dz8PocmJHiNWDgx4PwNHowBL+di7iGCnOzqFweNwOl5mACTQugTPCW2u2he774JC/j3Rfa3FbaLuDvOnkl8KsaUnueT046FOSDNBtMwETSB4BnPKwTOi+F474u5mQ+Vcccxvr9xM9H74oLPzTCWG195LX6n9uUT046H9utUtMwARMII9Ac3hxYP8weHcc8b44ZT1b0Y3lUzyafPyi8PnlPwur6knlurKkOGjNkdGjCLm/s++6Ivu9ywRMIIUEjgwv9lojDP4mAbLuPe2GehMt62nk01vCossmh1HP1TOWpDjoMUA8Av0e6THvfPODKvlEvG0CKSXQHJqb+oTDt+nGgIC2kBmPQ9YzE2+z/tu20DJ1YhiuJ5AbwpLioI+Epp5qlA5vCLLuhAmYQFkJnMmwOG7skVMO+5JLHknleqDtuiU45UfCPXdMC3u3lPWECagsKQ5aKCaj/0H9kebjsJmACaScAGOV1yNq+3ec8n/glNfJhLbPQXJza2id+mlYcH1zWFPPTzSsJclByylrzlWbCZhAigmcEmav2jv03ntptJwZC4oWbvzdw0zEp7WEj68+LqytidRSYUly0KkA7k6agAn8M4FTw6zB3ULf7xAtK6+8DUew2vYwOqY1tF0xKQx7658/1fgldtCNf43dQxNIJIHm8GbffqE7s1Xqyb7Mt0hj9CRqfo40xsmkMC6fFIZr4rRUmx10qi+/O28C1SXQHGZ07xs23J7xybrZ920cc3+c8hxSGOexnDopDH20ui1K9tnsoJN9fdw6E2gEApnTwttbdV8aKe+DUx5Cp97HKV9BCoNhcUPJL3vGyrgLbQcdR8VlJmACK0yAx603Yjp3PUCiiYnWwCF/yvL6Fpzyy2HWzReFsYtX+CQNXoEddINfYHfPBKpJIDsxEXf4GBYXNkTRxESZExmMMf3YMCLuQbRqNrGuzmUHXVeXy401geQR0MREzBY3XiMwaN1XiZJJJ7f9lVRGXU1MlDyyS18nlcR2uU0mYAIJJsADJAO6hcwepC00W9wOOOVurPPS58xxn4fPp9bjxERJxO0IOolXxW0ygQQSKDYxEZmMyyeGkc8msNl13SQ76Lq+fG68CVSWQO7ERETHe3G2QaghJyaqLMmu1W4H3TVu/pQJNDSBuImJSCxf19LAExMl8YLaQSfxqrhNJlADAtHERMwUp5t9mXXTNjFRDZB3eEo76A4R+QATaFwCmpioV+it9/XJKWs6z1Yc8wwetT6tNXx8TZomJkriVbaDTuJVcZtMoIIE8icmwjk3MQLjIYbGTWB55cSUTkxUQeRdrtoOusvo/EETqB8CcRMT4ZD/romJFvNk3/Fh2Ev105v0tNQOOj3X2j1NGQFPTFT/F9wOuv6voXtgArkEPDFRLo06X7eDrvML6OabgAicFt76UvfQxGT3TfuSUx5NLtkTEzXAn4YddANcRHchnQRODXPHdA/dGBbXpLmVN2IExhLyyrfx5tQTm0LLdZ6YqP7/Luyg6/8augcpIpA/MRFdb5+YqDVkDvs0LJrWHEa9myIcDd9VO+iGv8TuYL0TyJ2YiLHK2zNrXPfsxEQtYdEVx4VRr9d7H93+eAJ20PFcXGoCNSXQHJ7p2S8M/xYBstIXu9KYPuhV1n/liYlqemmqenI76Kri9slMoDABTUzULxyxLWOTySsvNzHRxUx2z6uhhj9Q+NPe04gE7KAb8aq6T3VFgBTGWCJjRcp65HoUzvkjEsvTmZjoikfCPXdMC3tz38+WRgJ20Gm86u5zzQlMCfO+0GPpG0g04f26pDI+bwuZm5gKY+onYcENzWHNhTVvpBtQcwJJdtDcCwl90Sc1p+QGmEAZCHhiojJATFkVSXHQA+F+ENoGnY0GIPJuoRe6Eh2D7KiBYKsvAsUmJloUwp9ODEPn1leP3NpqEkiKgz6OTq+D+IkX/gupXbuh55Ec9nh0KbKZQOIJTAiz+4wKvRh5kdEUnt8ir9yTFIYnJkr8lUteA5PioHcHzZZIr2Qfjoag7B3rU1mXk7aDBoItmQQ0VpnRF98kn7wnDnkXHHN/bvTNYbzyeSynTgpDH01my92qJBNIioN+Dkg7oBloa6Qxn1nbmBX/cWdpeJkYAlPCnFV6hJ674Yy/jVPekUeuScm18b6+cDmO+Qqc8kway6bNBLpGICkO+iyafwlaC/0aKQctp/0E+hraFtlMoOYETg9zVsuE3rsTLRMpZ7ahQd3Qa62h7XxGYEz/NJz/V8Yzt9a8oW5AQxBIioNWOmMDtDJ6DxGJhG8gvUF4f/QZsplATQhoSFz3kNmDKFnpC6XiMgyJexYvfCoB8nSnL2pyWVJx0qQ4aMHWT0E5Z9nn6C/ta/7HBGpA4FfhnU2Wpi7ac8ob0YQ20hYPMynRCUtC27UnhKEv1KBZPmXKCCTJQRdCvx47NB76sUIH5JRrqB53zmNtbUpfjd3jwtQT0GPWvcOhX2kKTYqSlb5YEygtOOV7SV8c1RYWTZ8cVpuTelAGUFUC9eCgNcRuNDqwBDIXcYwUZ+dQqBEiNhNoJ3BQeLjH2mF15r5od8qkMDIjCJR5oi/c0RZaf7kgLPmzp+/0H0stCSTRQatNukn4QQTmlFoC8rkbi4Benton9NyxaelNvl3p3SCi5I9Z3sxyOnnlGyeHodq2mUDNCSTFQTOQPzSj/dCqiF+Z4VOklMRZ6FJkM4EuEWgOrw7qE/rv3MQYZSJkxiqHvtzke49IeTrb018PH952XlhX9z1sJpAoAklx0OdBhZ+XYWf0CtIDKwORRnaci3qjC5DNBEoi8Mvw9vBeoWl3vutxymE71AOHPIftS+SYHwx33eNZ4kpC6YNqSCApDppB/mEr9FYOi/msP4CORicjO2gg2AoT4HVQa2ZCtz04Yk8i5K+ybCJCfoHUxdkaozwxjHiQMvyzzQTqg0BSHPTT4BqHrojBxmOzYV5MuYtMIJwZ5vIrqxtRsiLlzKZCgkN+nAdJTm4LLdfilPW3ZTOBuiSQFAd9EvSmognoZfQRWgmtj9TGnZDNBEQgc2Z4dywRcdYpf4GQmHt7bQ/glI+NnPKrRmUCjUAgKQ76MWAq+lGaYwxSPlpRs9IaM5F/lgIhrTY+XNVtyzDu65qIiD+Eb8NhNbIXi1mfgVM+h4j5uklhWG56LK2o3O8GI5AUBy2sC9GMBuPr7nSRwJHhxV5rhEHbR06ZuS/CKjhiRvZkbqXKEz4JH1/PW0c+7GL1/pgJ1AWBJDnougDmRlaOgKbsZHwy8ye3R8o7sdR4eJxw2w26ybcgtNzCgyMafmkzgVQQsINOxWVObid5cGRIv9B9VzllWsmUs029iJTfImK+Ygk3+V4Js++6KIxdnNweuGUmUDkCdtCVY+uaCxDQlJ1NoacerdaNvq05rBtO+VVyyue3htZrF4bz7/eUnQXguThVBOygU3W5a9fZpW+xbp/YXpHyFoggOTzD8ItTW3DKx4dhulFsMwETyCFgB50Dw6vlJaApOxUlR+mLDamdIDk8hFP2lJ3lRe3aGpRAqQ5aw97eRUsalIO7VQYCmrKzbzh0K80OR3W8BiqzJssW/PJMpuz8H0/ZWQbIriJVBIo5aB6TDSegvZF+jh6LDkWa9tNP9gHBFoKm7FwnrD4uO2UnfyrDccjtU3aSTz5lSVj05xPCau+ZlQmYQOcJFHPQmvx+O6RoaDq6C/GCTP6fDOGXyJZSApqys1/o9g1GXCh9oUfx26fsJH/BlJ2t1zIPxk2esjOlfxzudlkJFHPQX+dMZ6I3ozNqqJNmlrsQ2UFHUNKy0JSd/cIAOWO+sNu+gWPuy/JdnPJ0RmBcOzt8eLun7EzLX4P7WS0CxRz0bBohJ313TmN2Z31uzrZXG5gAN/lG4Ii55u3v5RtHV5dN2dmCU344zJjpKTsb+A/AXas5gWIO+hxa9xDi4YEwEmnqzzFoe2RrUAKaspNhyZrvgkg5ww2/f0zZ2cqDI5PDSP1NaDSGzQRMoMIEijloPVL7ZaQc9BroHvQ3xE0gWyMROC3M3bBb6K6XpeKY/zFlJ364eUlYcu1xYeQzjdRf98UE6oVAnIPuSeM1guMUdAu6HGVNk6Hvh3bOFnhZlwQyp4e5WzQxjzIpDDnl9ik7ySXzK6ntp3LKx4eRr9Vlz9xoE2ggAnEO+kf0T9N8yo5aulj2r16medyyLa/UDQFN2Tk2jNt66Xv5Ao9Za8rO0D5lJ0757EWh9boTw/C366ZDbqgJpIBAnIPWKI3fIjnn+9HDSKa8Iw8d2OqFgKbsXD0M2oEoWWmq3XDKy6bs5GIevyB8fIOn7KyXq+l2ppFAnIMWBz0xyHvcYq0PpZ/F7nFhzQlcRaQ8K2w3ie9THrPWm2gy/WnUh6QvbtAYZabsvNVTdtb8MrkBJlASgUIOWh8eghRNr4u6IeWl9XZt3Sj8LrIlkMCssO3ZOGZ+/WTewSlfrldAvRRmz/CUnQm8WG6SCXRAoJiDnsBneRgh/C/SAwpKe/wMTUG2BBLgYRK+QDN7kb54aEH4zVc8ZWcCL5KbZAKdIKCouJCtzY6z0O/QquhqtD/6KbIlkEDf0O9g8s2jSG9MsnNO4AVyk0ygkwSKOeg3qEvjnz9BGnq3CnofqcyWMAKKnhmhMZm0xt0Tw9C7E9Y8N8cETKALBIqlOC6mPj09+BL6C+K9cO2OehpLW8II9Av9DyG9MRIHvW/CmubmmIAJdJFAMQf9LHWuhzS0To4aB6AXeIYXkS1BBBQ9c2NQIzdmTApD70lQ09wUEzCBFSAQl+LQjcEzkZ4i3BW9i5agO9DB6BhkSxCB/qE/83RnmC+lrTlBzXJTTMAEVpBAXAStSZL0zjg94q1J+mehQeh3SKmOo5EtIQQmhNl9GLWh6PmuiWHYzIQ0y80wARMoA4E4B70V9SpSVlpDk+ScgYah76HpyJYgAqNC70MZuTGCt5fsnaBmuSkmYAJlIBDnoDW16GNR3S+z3ABtip6Oyqq1YExve/57cbVOWG/nUfQc5Z7vnBSG3Vtv7Xd7TcAEihOIy0Hz/3xojT6m3PPrqNLOWUP3/oDGoqFII0jeQropeQnSMD9bHoFVQ+/DyD0PX+Lccx4Zb5pAYxCIc9Dq2ZBIK0fdzG5rOSAqK+fiF1SmLwKlVI5Eiuw3Qhsjne8/kS2HgN4LyOZE8s93HBeG3Zezy6smYAINQiAuxaGu6SGVXJuXs6Fx0OXOd25NnV9EixDzE4c90Bwkk3PWnCC2HAJ9Q4/26Jm5Nr6TU+xVEzCBBiIQF0GvSf9WKqIfVKD/L1Dn96N672a5U7SuxS7IY69zgCh6Jg81kaLbJ4Xhf83Z5VUTMIEGIhAXQX9cg/4dzjlvQD9GenJR47B/hFrRQKQI2xYR6Bd6HM7IjWFLQkuzoZiACTQugTgHXYveZkeL7MDJ10PKR3+AFDnfiHSz0gYBXuraj4Wi51uPC8P1QgWbCZhAgxJIioMWXu53hdsiaTtrcti6IZYd+pctj1vuQ6GefoyzzSicG7ejnsraQvfDyUsNbQlLmuup3W6rCZhA5wmU6qBHUHX2ke/On2XFPjGej49GB5ZQzc0cc1+B435OuZ6IrFtbGj23Hcs32S2Tw4i/1W1H3HATMIGSCBRz0LqBeALSiA2NjT4WMedDu6PMHdVBUUXtlE7U/hHHSnG2gML+cTvqp6z7EVyIoW2OnuvnkrmlJrACBOSEC9lB7NgO7RkdcBdLDb9TeSVNXxqDK3mCeqxb0XMmKHpuu3liGPF/9dgHt9kETKBzBIo56K9TlUZTvBlVqUeuz0Vy2uW2nlQ4Bc1GGgutFwMo4n0a7Y9Sb+SeeYAnM6TV0XPq/xYMID0EijloOUs56VzbnY1K3Gg7j3o3RDsjDatTu0Yh5Z0PQUqtpNaawzv9m0LbT4meb5ocRj6YWhDuuAmkjECxHPQ5sHgIaeibJlB6AI1B26Ny245UqFn0NP9G1uazonMejU5GF6BUWr/2x98zQ8g/N6cSgDttAiklUMxBvw2TDZCGrq2B7onUwrLcplTGOHRFTMW7UDYvpjwVRYqe+UGh6PlG3jWoL0ybCZhASggUc9CHw0Djj/+IciPbSqA5iUqnoglID61oJMZKaH2kNu6EUmlcgKOInFfhW7E5lQDcaRNIMYFiDvpeuGhmOUW3SjVcim5AuolXbnuMCjdFSnOMQRp3rahZaY2ZiKG/6bPTw7wBvKn7p3T/hslh6MPpI+Aem0C6CRRz0E+CRjfp5KR3Qz9Av0anov9G5baFVDij3JXWc33cGDyKkRsr04fmeu6H224CJtA1AsVGcWRr7MaKhsFpqTkxPC8GECptip4zoekYcs/XHxuGPlLp87l+EzCB5BEoFkFvTnOVE94Z3Y0uQjchO2ggVNq6hTZGr2RWJrfTXOlzuX4TMIFkEijmoLekyY+iY9A7yWx+Y7aqObw4sI3omScH/zIpDNU1sJmACaSQQDEHndpxx7X+O+gXBh3NyI3Bi/2uwVpfCp/fBGpKIM5Ba56HSeiraP+Y1t1MmR4esVWAgKJnUhukltquOz4M0+gWmwmYQEoJxDnog2DxGtLYZ02QlJtz1s3CHshWIQL9w6CfUPWg1tB2coVO4WpNwATqhEDcKI7naPvn6DCkIV7P5Gg065ORrQIETgsvr6TomRuD100Kwx6vwClcpQmYQB0RiIugf0T7s/lnxuEuZ3pf4XHLlXijbAS6hYGKnldqCS2OnstG1RWZQP0SiIugL6Q7SmPwBFv7k31al+TMyY+G85GtzAQUPfMi2J8QPU/nXYNPlLl6V2cCJlCHBOIiaHVDeeez67A/ddtkomeNOSfF4ei5bi+iG24CZSYQ56A9iqPMkDuqrjm8Oojcs6LnayeG4U92dLz3m4AJpINAnIPOjuLQrHK3xWD4IKbMRStAoG/oP4FxzwOXhCXOPa8AR3/UBBqNQJyDzuY/59NZgrr2pwh7s9TNw9eR54UAQrns1DBrMLnnowF9zXFhxFPlqtf1mIAJ1D+BOAed7ZUeVLkdrYc0X/NYpHHQGnp3MbKVgUD30Lc9euZN3Y6ey8DTVZhAIxGIG8WR7d9/sHIA0ptV9kbfRyrbC9nKQCAbPVPV1bypW/Nu20zABExgGYFiDpoRBe2T5n+dpSZLkgPphT5CtjIQ6Bb6HUM1A5aExY6ey8DTVZhAoxEoluK4kc6eizQP9O/QBuj36JfItoIEmsPslfl2PIr5nqcdF0bqaU2bCZiACSxHoJiDnsqR7yKGgIWr0VroUDQD2VaQQL/Qh8n4Q/+Mc88rSNIfN4HGJVDMQavXGma3CvoSehG9hGwrSEDRMwNk9DLYaceGkc+uYHX+uAmYQIMSKJaDVmrjUqQoWrPaaR6Oa5Dy0LYVINA/9OEx+ky/FueeV4CiP2oCjU+gmIM+mO6vg5R7VhS9LiLo82x2MOiyTQlzViHvfCRvS/nT5DBKMwfaTMAETCCWQDEH/RU+cQbKOpFXWP8F2gbZukigR+idjZ7/Xxer8MdMwARSQqCYg74fBhpil2vanpdb4PXSCSh6JnI+kvzzlY6eS+fmI00grQSK3SScBpSnkCLm+9DmaJNom4WtswR6hp7HtoVM38UhOHruLDwfbwIpJFAsgn4PHhq9cRnScRoXvRF6HFXDNP+H5p9uCGsObw4hhX8EEfQVJ4Qhf2+ITrkTJmACFSVQKIIey1m3QLeiX1e0BYUr/w67xqEDCh9SP3v6he7H4qD7tIQ2R8/1c9ncUhOoKYG4CFpzbtyL9kN6wu10VGnTGGtNY5qri9hWG1Sm4X51a1PC3KGKnhm9ccXkMPT5uu2IG24CJlBVAnEOmptY7ZMiaTa7f0GMOqj42Of9OYduPp6DlOeWTkDTo/WJLOvWehI9M6Vob+ee6/YSuuEmUBMCcQ56VVqit6rInkRynKtro4Kmm5BKq6yD9KqtBUgPyHyCZkXrLOrPeFtKb24MHsjIjatOCENfqL8euMUmYAK1IhCXg9aczwR7y0wpBt2wq7R9xAmUXtHUpjORviRaUF1b/zBgPB0YjIPWy3htJmACJlAygbgIWh/W6IlBkbqxHJCz3Zf1StpVVL4jGoLequSJqlE3eWdFzy8cG4bpS8dmAiZgAiUTiIug9eH8n+J6aCVrGh+tKLeSNofKd41OoDe66EvhsWi72GIHdn69wAF6MlJRetVsSnj3izwbT3vaGMFhMwETMIHOEYhz0Nk5NwrVtKjQjgqVK0UwGhGJdmgaDbKwwFFy9ErfVM16ED2Tf160OCz5Q9VO6hOZgAk0DIE4B62cc5LslE405jWOleJsTwqHx+2oRFlzeKYn0fP3SXFcd0IYqRutNhMwARPoFIFCOehOVVLmg/WlwU21+rZ+YRhfCBny6K3/W989cetNwARqRSApDlqphyloNlIK5X2koXZPI42RrkfTzcFXJobhd9Zj491mEzCB2hOIS3HEtWoEhRqXvCRuZxnKzqMOnWNnpGlN5Zw1kkRzUZ+LNMzvAlQX9qvw9to8mDKuLbSeSIPb6qLRbqQJmEDiCBSLoLXvZ+hJdDv6N3QdGorKbTtSoV4QoHN9guTU5qMH0NFoD1Q3hnPm5mBoQZfWTaPdUBMwgcQRKOagD6K12yHdXJPdhd5AKi+3KZUxrkClu1BeNzfZDgoP98BB/5A2Xz8pDKv7cdwFromLTcAEqkCgWIpD44nPRG9G7VjMUumGC9Evo7JyLU6ioqloAnoZfYRWQusjtXEnVBe2Tlh9N24ODs/45mBdXC830gSSTKCYg9YNOznpu3M6sDvrc3O2y7X6GBVtirZCY5Dy0YqalXeeiZTyqAtrCt2U3nh9QThfU7XaTMAETKDLBIo56HOo9SG0AxqJlA8eg7ZHlbCFVDqjEhVXq84zw5ujcc47MCn/yc2hubVa5/V5TMAEGpNAMQf9Nl3WKIp90BronkgtLG0xBFpDzwNI6re1hM8vidntIhMwARPoFIFiDlrzQO+XU5vSG7Lb0KT2Nf+zjMD4cFU3nPP+PDl4y+Sw2pxlO7xiAiZgAl0kUMxBT6fOB6N6MyxHoaPRTVGZFzkEtgzjdCNzVfIah+UUe9UETMAEukygmIN+hVqlXNP2seju3EKvtxM4lOj5zYfDjBvNwwRMwATKQYBf5Z2yNTl6pU59IgUH/yq8sy5jn79JVy+cFvZ2jj4F19xdNIFqECgWQStS/n5OI/qwvjraN6fMqxDIhCZeCKu30LT8j4GYgAmYQLkIFHPQN3ASDa3LmubhUIqjbp7qyza8ksvTw7wBJOh/SHrjTxPDiHcqeS7XbQImkC4CxRz0waDQULvT0oWkc71lzPMPiaEHtobw68590kebgAmYQHECxXLQs/jol5DeSWiLJ5BpChmlNx6YHIY+HH+IS03ABEygawSKRdCfUaUmKtK8GHrsO3vzS48wH4NSb9wc/AbR8xd4NsV5+dT/NRiACZSfQDEHfQuneyLmlO/FlKWyiJEbR2lo3Yth1jWpBOBOm4AJVJRAnIPWfM9ywkpxSLYYAtmhdTjon18UxmqmP5sJmIAJlJVAXA76Oc5QiUn5y9rwWlfWFJqOJLXB67k8tK7W18LnN4FGJRAXQTdqX8vWLw2tI3L+IflnD60rG1VXZAImkE+gkIPeigM/yD842tY46GcL7EtFMUPr/pWHUwa0hLbfpaLD7qQJmEBNCBRy0OfTmuyojfyG6QGWQ/ML07Wd2Uj9bQ2fPp6ufru3JmAC1SRQyEHr7SZ6SMUWQ4DRG19i7POc48PoQr8yYj7lIhMwARPoHIG4m4SdqyGFR5Pi+BI3CJ9KYdfdZRMwgSoSiHPQd3L+z6vYhro6lSbmJ3pen0Y/XVcNd2NNwATqjkBcikOvuLIVILBp2I6pRUNvdjuCLsDIxSZgAuUhEBdBl6fmBq2le2glvaEbhK120A16jd0tE0gKATvoTl+J9hEcLbPDfD3QYzMBEzCBihGwg+4k2mgExwvnhXWdp+8kOx9uAibQOQJ20J3jxdtTwqaM4niykx/z4SZgAibQaQL14KA1H3WvTvesAh84JbyhV36NaQ1tf61A9a7SBEzABJYjkBQHLcf3B/QJuh2tg7I2npU/ZjdquewVem2j8+OgZ9ayHT63CZhAOggkxUFPAPdcNBY9gOQAmQg/ada2NWOgPzguDHeKI2mXxu0xgQYkEDcOuhbd3ImT6vFyvcXlJKTJmG5FX0OJMfLPRNBt99Ig/LTNBEzABCpLICkRtByyouesXcnKeehmtEq2sJbLM8PcDfR6KzzzPbVsh89tAiaQHgJJcdAXgnwampyD/mzWr0Hn5JTVZFWPd4fQfSpzQM9DU2vSCJ/UBEwgdQSSkuK4DfJro7XyrsDJbCti1b6a2djwNfLhmS8zvO7giWHYWzVriE9sAiaQKgJJcdCCvgDFPT6tm4fzdUAJthnHtD+KHXOsJjjq0rsDu4Xu7a8AawmZl2PqdZEJmIAJVIRAUlIcxTqnYXaHFTug8vvahugcTaFFb5OxmYAJmEBVCCQpgi7U4VMK7Ygpf5QyKc42oXB43I6OytpC01BGcIRFockOuiNY3m8CJlA2AkmMoPWlMbhsPSxLRZnRVLP48TDjnbJU50pMwARMoAQCSXHQPWnrFDQbLULvI+Wkn0b7o5oaNweVv35hWti7paYN8clNwARSRSApKQ6NeR6BdkavIDnngYixx+FcpAnyL0C1MrXjiVqd3Oc1ARNIJ4GkRNA7gv9gpEeoP0F6Uk8jN/TY99FoD1QTOzK82Ishdmsz/lkP09hMwARMoGoEkuKglcoYV6DXu1Bes5tza4RBW3H+bpnQ6vk3ClwgF5uACVSGQFJSHJp/YyqagDTW+CO0ElLuV23cCdXIMicSzi9aEFpm1qgBPq0JmEBKCSTFQT8G/02RotUxSPloRc3KO8sxKuVRE+MGYW9Ofl9zGPVuTRrgk5qACaSWQFIctC7AQjQjaVcC51yzL4eksXB7TMAEqksgKTno6va6E2fjARX+s5mACZhA9QnYQRdhvnQWu8xqeOhPixzmXSZgAiZQEQJJSnFUpINdrbQ5PNOzfxj2HJ8f0xZaj+lqPf6cCZiACXSVgCPoAuR6hf79yG6sxfjnKyeG4dcVOMzFJmACJlAxAnbQHaJt08MyvlHYIScfYAImUG4CdtDlJur6TMAETKBMBOygywTS1ZiACZhAuQnYQZebqOszARMwgTIRsIMuE0hXYwImYALlJmAHXW6irs8ETMAEykTADrpMIF2NCZiACZSbgB10AaLdQ99vFNjlYhMwAROoCgE76IKYMz9m+PNnvMtbc1XbTMAETKDqBOygY5A3h1d7M//GOHbdPzEMuSvmEBeZgAmYQMUJ2EHHIu6nOUq6tYZwS+xuF5qACZhAFQjYQcdA5iWEemFAAI4f8Y7h4yITMIHqELCDjuHcPWS+rWK88wsxu11kAiZgAlUhYAedh7k5zCC9kZlCMf655f/ydnvTBEzABKpGwA46D/V7YTXe4N3+otpTJ4YR7+Tt9qYJmIAJVI2AHXQe6tGhX3v+mRuEH+ft8qYJmIAJVJWAHXQe7pbQ42sqagotzj/nsfGmCZhAdQnYQRfgvThkniywy8UmYAImUBUCdtBVweyTmIAJmEDnCdhBd56ZP2ECJmACVSGQRAetp/gGV6X3PokJmIAJJJhAUhx0Txhp7PFstAi9jxYgTVS0P7KZgAmYQOoIKFpNgp1HIzS8bWf0CpJzHog2QOei3ugCVHFjDPTmOkm30ObHvCtO2ycwARMoRiApEfSONPJgpJETnyA5x/noAXQ02gNVxQCyJief9VC457WqnNAnMQETMIECBJLioJXK0PSecbYLhfPidlSurG3+tLB3S+Xqd80mYAIm0DGBpKQ4TqKpU9EE9DL6CK2E1kdq406oKtYWMgTQzm5UBbZPYgImUJRAUhz0Y7RyU7QVGoOUj1bUrLzzTGSPCQSbCZhAuggkxUGL+kI0Iwb/epT1RXLiHdkYDli9wEGjKO9ZYJ+LTcAETCBxBJLkoAvBGc+O0ejAQgfklK/L+tdztnNXNbb6zdyC+PXWh5hu9K34fS41ARMwAROoBIF9qPSQSlTsOk3ABBqawFn0brNa9DApozhy+66oXtGuzQRMwARSTSApDlq54SnITxKm+s/RnTcBE8glkJQcdGKeJMyF43UTMAETqCWBpETQOwIhEU8S1vJi+NwmYAImkEsgKQ46YU8S5iLyugmYgAnUhkBSUhyJeZKwNpfBZzUBEzCBfybA5G2JMc1Yl/8k4YuUzUTleJJwE+q5EXX0wIvGXEua9rSeTb+OeqHP6rkTUdv7sdQMh/Vu7keyrmAfmvMg6uhvay2O2QG9gWw1JrAb59ecIPVuG9EB3Xytd9NTpPpibQS7uxE6QR8apR+X0ZdVk3xNkpKDTjIjt80ETMAEakLADrom2H1SEzABE+iYgB10x4x8hAmYgAnUhIAddE2w+6QmYAIm0DEBO+iOGfkIEzABE6gJATvommD3SU3ABEzABLpCQMO69LqtercedGCVeu8E7ddYfb1hpxFsZCN0gj40Sj+G0pduDXJN3A0TMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMIEqEtDDW7Y6IXA87XwSvYq0XshKPa7Q5ytdXkr79Id5Bno40qkse6IkWSn9yG2vXjT8fm5BQta3pR33If1dTUeDUZz9C4W6Hs+gG9D6KEm2LY0ppR/f4rg7kN5U9Fukp/OSaPvSKPEuZLpOVyG9xekp9FVkqxGB8ZxXf3x6pFuPET+O9IeWb6Uel/+5am2X2r4DaNC1SI5a+jNSWVKs1H5k26v/mfTl+kG2ICHLIbTjTbQxEuez0SUo3/QqspfRV6Id/87y6mg9CYtS+6FXRs1BX4wafTrLs6L1pCz0t/Ib9A56pEij5Jx/hjS1wLboLaT+2WpA4GLOeUjOeSezflHOdna11OOyx1d7WWr7tqBha+c0ThH0pTnbtV4ttR/Zdl7Gyv4oaRH0N2nTndlGslwTfZiznV39Niu3RhsKEpJmpfZD71hciPQOT5kc3BXta8n5Zy+a8iukAKyYg/6I/SujrCna3iG7UatlWmezWwPgc3Og69tyeM52drXU47LHV3tZavseomEvR43T/1TfRTdE20lYlNoPtVXRtpxCriNUeRIsvx9v0yg5YEXMuSaHpi+XmWge0rXZECXFSu3HAhp8JLoHXYO+h05GSTL9MpmEir08WVG2rlHuF758wjBUU0urg9Ysb/rjytqnrMhx5Vupx+V/rlrbnW2f8s5XooeQ/odKipXajxE0+CR0bFIanteO/H5knULfvOOGsr03uhDpM7cg/YpLipXaj940WPcClG9/HimtsymqN8vvr9qva9e/1h3pXusG1Oj873LegTnn1rpyh/lW6nH5n6vWdmfaJ+d8LdKXsiLoJFmp/fhvGn0f+hpSdKM+7YJuR5+jWpv6sXFOIwawrmg/P1eutMczaCqSKXernLr6swjV2krtx7/R0E3QulGDb2apL34FAW1RWT0s8v/+1OZCPqGq/Umrg54D5WzeTMDHoNlaybNSj8v7WNU2S22frrP+p9Hct8p/JsEJ0IxlVmo/1O4vR9JPUkVwJ6D7URIctPoxBmVtDCuF/q4+zh7EcjHSDamk/KIttR+r0+a/oaw9yoocmzQ/W1gHS31hKmJeDanvsjHoda3Yqk/gm5zyCTQKjUEaWjMWydZHG7WvhVDsuOiQmi6KtS+3Hz+hlervSKQbIVJ/lBQr1o/BNHL7mIYqT/p+THkti/Sl8TZSZKn136NTkSy3H4qsFbVtqR3Yieie9rVk/FNqP5SqkRPT35XsCHR9+1ry/tmWJj2S1yyVDYnKLmb5a9QdfQc9h5SysdWAQIZzXoL003MuakZZO4+VS6ONYsdlj6/lslj7cvvxGo3UT85c3VjLhuedu1g/tubYBXnHazOJDlrtGo8+RorE7kLZL8L8fuzGPv3tvYT+jtZESbJS+3EIjdY9DTk03XjeFCXRtqVR+Q76Lcp2iho7huVT6A2ka7ItstWYgH6KKVroyEo9rqN6KrU/6e0rtd+N0g9FYYqYOzJ9MQ3t6KAa7i+1H2qirl0jWJKvRyPwdR9MwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwARMwASSQ2BfmrIEfZSnwzto4rns/88Ojil1950c+BlSG+YjvS17KirlnYAc9k+2OSXPR6WrsTwyWs8tj4q6vPgxn1yIZkd6k6XeLP9z1IQ6spM5oGdHB3m/CZhAugnIQf9fFxCU20EfmNOGvqzfiH6bU9aZVTm+kdEH1L/ro/Xc8qioyws56Bl5n/4i2x+ib+aV5292o0BvZ++dv8PbjUeglG/rxuu1e1QtAj/gRM+hT9CjaAuUb/tR8Dp6D01D2chXb1e+FslpPYG2RqXYpxwk56foV7YqOhMpSn0SbYdkemP1JUj1z0KTkWxddCkahc5C26LLULa8P+tqzyCUtfNZ2RPpLd0/Q3PQG+hEpLJS7O8c9BJaKTp4A5bqx3yk9k1AsiuXLtrbMIT1rnKKqvHCBEygUQkownwVKYLNSo5KJocmx7wpGowuQrci2blIKQ5FgR+jTZAc3k3oeCRT5CpHOQLtj15GcaYUh84tkzPcDD2BDkaye5DqkaM+COl8cr77oJloFbQ++gitgzZHz6MmdAC6BckpZ8tZDTcjfbHIFLHPR3KW30dytOqzvoyeRv+C8k0R9ONo+0h7sPwN+gD1QzJ9oU1E2hbTJWhlJE5tSFzU3+tRKZw4zGYCJpAmAnLQck635ejCCMBAlhtG64oKf4KejLbPZSkH3QstQMeg4UhpBJkcUQuS41Q90r1oY5RvctCKmj9EqkvO6w4kR7ZatC3nnLXHWDkc7YleRbsitUOS5Tpi9U8OUJZbLkd8XXtpCHuxlBOXqS2TUbbN6tcpKN/koD9BTyG1QW2Wg1a/szaWFaUzuqPN0MdoPaQyHa8vNx1fKicOtdUbAV18mwmsCAFFjDvGVCCHsk8kOSz9fFdUmmufs7E3moJ+hWagI5CcpZzQXSjXvsrGk7kF0foZLK9BclZzkL40ZPqCmI3e0EZkf2OpCPoCJKd7MdIXwx+QItZSTM75v5AiaznobNpBXwSq4ycoa49mV/KWD7E9LirTZ/Sl0RxtazEU3Yu+iB5Hcsz5/LJfQKVyogqbCZhAWggowix0k/AH7HsWbRTB2I3lM9H6uSz/E8nhjIzK5EzlZBX9ymG+h4agrMlhqTzf7qQgm+LI3/cFCpagQTk75PQOQH2QIlAFKTuhV9HBSE77eSRT/65vX1u+XEVq6/fQPJTNG/+R9dy2yIFn97G6zH7Mmr6Mck39+HNUoHbpV8F3kCJlpTI+Q8pLy1Hry0vlneHE4bZ6I5D/jVxv7Xd7k0tATuZF9DSSg/kh6oFyTQ5Y+xUJynnfjGSLkByWokr9jSrfKmevaLIz9jIHy/HKIasNmyA54L+hf0dXITk7nTfrlFldZgtYi3OwOuBK9At0H8pG7HKw+6PBSOe7DE1ApdhhHPQttAuSY5fpy2oh0heFHLL4taDPkdpVLk5UZTMBE2g0AnIchSLo4ex7BD2B5FxPQB+jvigbQbPa7sBeYykHPQttgWSboheQyl5Bk1GcyZHnRq35x2xFgZz06+gdJAcqU+T8Z/Qamo1uQgOQHPjzSDYGvYseR7nlbLZH4OrPPtqITH27Fs1HL6HrkSL1fIuLoHXMGUhfakrxXIBeQ2L4R/QAGo9kN6PFSL86SuXEoTYTMAETWJ7AKmyW8ktt6PIfW7al8syyra6vDCtQjyJT5cgLmdoe52QLHa/yfpGKHVPKPtUjpx9n2pdr5eKUW6fXTcAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETMAETKD6BP4/0lGmcafEU88AAAAASUVORK5CYII=" alt="plot of chunk roc"/> </p>
<p>This analysis illustrates the use of Polyester to simulate data with known FPKM, fold changes between groups, and differential expression status between groups, using simulated transcript abundances based on real data.</p>
<h2>for reproducibility</h2>
<pre><code class="r">sessionInfo()
</code></pre>
<pre><code>## R version 3.1.1 (2014-07-10)
## Platform: x86_64-apple-darwin13.1.0 (64-bit)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] limma_3.22.1 knitr_1.8.2 EBSeq_1.6.0
## [4] gplots_2.15.0 blockmodeling_0.1.8 Biobase_2.26.0
## [7] GenomicRanges_1.18.3 GenomeInfoDb_1.2.3 IRanges_2.0.0
## [10] S4Vectors_0.4.0 BiocGenerics_0.12.1 ballgown_1.0.1
## [13] polyester_1.0.2 colorout_1.0-3 devtools_1.6.1
##
## loaded via a namespace (and not attached):
## [1] annotate_1.44.0 AnnotationDbi_1.28.1
## [3] base64enc_0.1-2 BatchJobs_1.5
## [5] BBmisc_1.8 BiocParallel_1.0.0
## [7] Biostrings_2.34.0 bitops_1.0-6
## [9] brew_1.0-6 caTools_1.17.1
## [11] checkmate_1.5.0 codetools_0.2-9
## [13] DBI_0.3.1 digest_0.6.4
## [15] evaluate_0.5.5 fail_1.2
## [17] foreach_1.4.2 formatR_1.0
## [19] gdata_2.13.3 genefilter_1.48.1
## [21] GenomicAlignments_1.2.1 grid_3.1.1
## [23] gtools_3.4.1 htmltools_0.2.6
## [25] iterators_1.0.7 KernSmooth_2.23-13
## [27] lattice_0.20-29 markdown_0.7.4
## [29] Matrix_1.1-4 mgcv_1.8-4
## [31] mime_0.2 nlme_3.1-118
## [33] RColorBrewer_1.1-2 RCurl_1.95-4.5
## [35] rmarkdown_0.3.12 Rsamtools_1.18.2
## [37] RSQLite_1.0.0 rtracklayer_1.26.2
## [39] sendmailR_1.2-1 splines_3.1.1
## [41] stringr_0.6.2 survival_2.37-7
## [43] sva_3.12.0 tools_3.1.1
## [45] XML_3.98-1.1 xtable_1.7-4
## [47] XVector_0.6.0 yaml_2.1.13
## [49] zlibbioc_1.12.0
</code></pre>
</body>
</html>