-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathR_vcf.html
762 lines (735 loc) · 57.1 KB
/
R_vcf.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
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.2.269">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<meta name="author" content="Solomon Chak">
<meta name="dcterms.date" content="2023-03-08">
<title>R VCF</title>
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { display: inline-block; line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
color: #aaaaaa;
}
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa; padding-left: 4px; }
div.sourceCode
{ }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { color: #008000; } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { color: #008000; font-weight: bold; } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
</style>
<script src="R_vcf_files/libs/clipboard/clipboard.min.js"></script>
<script src="R_vcf_files/libs/quarto-html/quarto.js"></script>
<script src="R_vcf_files/libs/quarto-html/popper.min.js"></script>
<script src="R_vcf_files/libs/quarto-html/tippy.umd.min.js"></script>
<script src="R_vcf_files/libs/quarto-html/anchor.min.js"></script>
<link href="R_vcf_files/libs/quarto-html/tippy.css" rel="stylesheet">
<link href="R_vcf_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
<script src="R_vcf_files/libs/bootstrap/bootstrap.min.js"></script>
<link href="R_vcf_files/libs/bootstrap/bootstrap-icons.css" rel="stylesheet">
<link href="R_vcf_files/libs/bootstrap/bootstrap.min.css" rel="stylesheet" id="quarto-bootstrap" data-mode="light">
</head>
<body>
<div id="quarto-content" class="page-columns page-rows-contents page-layout-article">
<div id="quarto-margin-sidebar" class="sidebar margin-sidebar">
<nav id="TOC" role="doc-toc" class="toc-active">
<h2 id="toc-title">Table of contents</h2>
<ul>
<li><a href="#load-libraries" id="toc-load-libraries" class="nav-link active" data-scroll-target="#load-libraries"><span class="toc-section-number">1</span> Load libraries</a></li>
<li><a href="#load-the-vcf-file" id="toc-load-the-vcf-file" class="nav-link" data-scroll-target="#load-the-vcf-file"><span class="toc-section-number">2</span> Load the vcf file</a></li>
<li><a href="#examine-the-vcf-file" id="toc-examine-the-vcf-file" class="nav-link" data-scroll-target="#examine-the-vcf-file"><span class="toc-section-number">3</span> Examine the vcf file</a></li>
<li><a href="#analysis" id="toc-analysis" class="nav-link" data-scroll-target="#analysis"><span class="toc-section-number">4</span> Analysis</a></li>
<li><a href="#assignment" id="toc-assignment" class="nav-link" data-scroll-target="#assignment"><span class="toc-section-number">5</span> Assignment</a></li>
</ul>
</nav>
</div>
<main class="content" id="quarto-document-content">
<header id="title-block-header" class="quarto-title-block default">
<div class="quarto-title">
<h1 class="title">R VCF</h1>
</div>
<div class="quarto-title-meta">
<div>
<div class="quarto-title-meta-heading">Author</div>
<div class="quarto-title-meta-contents">
<p>Solomon Chak </p>
</div>
</div>
<div>
<div class="quarto-title-meta-heading">Published</div>
<div class="quarto-title-meta-contents">
<p class="date">March 8, 2023</p>
</div>
</div>
</div>
</header>
<section id="load-libraries" class="level1" data-number="1">
<h1 data-number="1"><span class="header-section-number">1</span> Load libraries</h1>
<div class="cell">
<div class="sourceCode cell-code" id="cb1"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(<span class="st">"dplyr"</span>)</span>
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(<span class="st">"tidyr"</span>)</span>
<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(<span class="st">"ggplot2"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
</section>
<section id="load-the-vcf-file" class="level1" data-number="2">
<h1 data-number="2"><span class="header-section-number">2</span> Load the vcf file</h1>
<p>See the <a href="https://datacarpentry.org/genomics-r-intro/02-data-prelude/index.html">Data Carpentery</a> for a description of the data.</p>
<p>We will load the data and call the dataframe <code>variants</code>.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a>variants <span class="ot">=</span> <span class="fu">read.csv</span>(<span class="fu">url</span>(<span class="st">"https://raw.githubusercontent.com/datacarpentry/genomics-r-intro/main/data/combined_tidy_vcf.csv"</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
</section>
<section id="examine-the-vcf-file" class="level1" data-number="3">
<h1 data-number="3"><span class="header-section-number">3</span> Examine the vcf file</h1>
<p>Let’s get an idea of what <code>variants</code> look like. Here’s a <a href="https://datacarpentry.org/genomics-r-intro/02-data-prelude/index.html#how-variant-calls-are-stored-in-vcf-files">brief description of the CVF format</a>.</p>
<p>We can preview the first few lines:</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb3"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a><span class="fu">head</span>(variants)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code> sample_id CHROM POS ID REF ALT QUAL FILTER INDEL IDV IMF
1 SRR2584863 CP000819.1 9972 NA T G 91 NA FALSE NA NA
2 SRR2584863 CP000819.1 263235 NA G T 85 NA FALSE NA NA
3 SRR2584863 CP000819.1 281923 NA G T 217 NA FALSE NA NA
4 SRR2584863 CP000819.1 433359 NA CTTTTTTT CTTTTTTTT 64 NA TRUE 12 1.0
5 SRR2584863 CP000819.1 473901 NA CCGC CCGCGC 228 NA TRUE 9 0.9
6 SRR2584863 CP000819.1 648692 NA C T 210 NA FALSE NA NA
DP VDB RPB MQB BQB MQSB SGB MQ0F ICB HOB AC AN DP4 MQ
1 4 0.0257451 NA NA NA NA -0.556411 0.000000 NA NA 1 1 0,0,0,4 60
2 6 0.0961330 1 1 1 NA -0.590765 0.166667 NA NA 1 1 0,1,0,5 33
3 10 0.7740830 NA NA NA 0.974597 -0.662043 0.000000 NA NA 1 1 0,0,4,5 60
4 12 0.4777040 NA NA NA 1.000000 -0.676189 0.000000 NA NA 1 1 0,1,3,8 60
5 10 0.6595050 NA NA NA 0.916482 -0.662043 0.000000 NA NA 1 1 1,0,2,7 60
6 10 0.2680140 NA NA NA 0.916482 -0.670168 0.000000 NA NA 1 1 0,0,7,3 60
Indiv gt_PL
1 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 121,0
2 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 112,0
3 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 247,0
4 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 91,0
5 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 255,0
6 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 240,0
gt_GT gt_GT_alleles
1 1 G
2 1 T
3 1 T
4 1 CTTTTTTTT
5 1 CCGCGC
6 1 T</code></pre>
</div>
</div>
<p>What other codes can be used to check the structure of each column? For example, to find out if they are numeric or categorical.</p>
<div class="cell">
<details>
<summary>Code</summary>
<div class="sourceCode cell-code" id="cb5"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a><span class="fu">summary</span>(variants)</span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a><span class="do">## sample_id CHROM POS ID </span></span>
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a><span class="do">## Length:801 Length:801 Min. : 1521 Mode:logical </span></span>
<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a><span class="do">## Class :character Class :character 1st Qu.:1115970 NA's:801 </span></span>
<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a><span class="do">## Mode :character Mode :character Median :2290361 </span></span>
<span id="cb5-6"><a href="#cb5-6" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean :2243682 </span></span>
<span id="cb5-7"><a href="#cb5-7" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:3317082 </span></span>
<span id="cb5-8"><a href="#cb5-8" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :4629225 </span></span>
<span id="cb5-9"><a href="#cb5-9" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb5-10"><a href="#cb5-10" aria-hidden="true" tabindex="-1"></a><span class="do">## REF ALT QUAL FILTER </span></span>
<span id="cb5-11"><a href="#cb5-11" aria-hidden="true" tabindex="-1"></a><span class="do">## Length:801 Length:801 Min. : 4.385 Mode:logical </span></span>
<span id="cb5-12"><a href="#cb5-12" aria-hidden="true" tabindex="-1"></a><span class="do">## Class :character Class :character 1st Qu.:139.000 NA's:801 </span></span>
<span id="cb5-13"><a href="#cb5-13" aria-hidden="true" tabindex="-1"></a><span class="do">## Mode :character Mode :character Median :195.000 </span></span>
<span id="cb5-14"><a href="#cb5-14" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean :172.276 </span></span>
<span id="cb5-15"><a href="#cb5-15" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:225.000 </span></span>
<span id="cb5-16"><a href="#cb5-16" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :228.000 </span></span>
<span id="cb5-17"><a href="#cb5-17" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb5-18"><a href="#cb5-18" aria-hidden="true" tabindex="-1"></a><span class="do">## INDEL IDV IMF DP </span></span>
<span id="cb5-19"><a href="#cb5-19" aria-hidden="true" tabindex="-1"></a><span class="do">## Mode :logical Min. : 2.000 Min. :0.5714 Min. : 2.00 </span></span>
<span id="cb5-20"><a href="#cb5-20" aria-hidden="true" tabindex="-1"></a><span class="do">## FALSE:700 1st Qu.: 7.000 1st Qu.:0.8824 1st Qu.: 7.00 </span></span>
<span id="cb5-21"><a href="#cb5-21" aria-hidden="true" tabindex="-1"></a><span class="do">## TRUE :101 Median : 9.000 Median :1.0000 Median :10.00 </span></span>
<span id="cb5-22"><a href="#cb5-22" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean : 9.396 Mean :0.9219 Mean :10.57 </span></span>
<span id="cb5-23"><a href="#cb5-23" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:11.000 3rd Qu.:1.0000 3rd Qu.:13.00 </span></span>
<span id="cb5-24"><a href="#cb5-24" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :20.000 Max. :1.0000 Max. :79.00 </span></span>
<span id="cb5-25"><a href="#cb5-25" aria-hidden="true" tabindex="-1"></a><span class="do">## NA's :700 NA's :700 </span></span>
<span id="cb5-26"><a href="#cb5-26" aria-hidden="true" tabindex="-1"></a><span class="do">## VDB RPB MQB BQB </span></span>
<span id="cb5-27"><a href="#cb5-27" aria-hidden="true" tabindex="-1"></a><span class="do">## Min. :0.0005387 Min. :0.0000 Min. :0.0000 Min. :0.1153 </span></span>
<span id="cb5-28"><a href="#cb5-28" aria-hidden="true" tabindex="-1"></a><span class="do">## 1st Qu.:0.2180410 1st Qu.:0.3776 1st Qu.:0.1070 1st Qu.:0.6963 </span></span>
<span id="cb5-29"><a href="#cb5-29" aria-hidden="true" tabindex="-1"></a><span class="do">## Median :0.4827410 Median :0.8663 Median :0.2872 Median :0.8615 </span></span>
<span id="cb5-30"><a href="#cb5-30" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean :0.4926291 Mean :0.6970 Mean :0.5330 Mean :0.7784 </span></span>
<span id="cb5-31"><a href="#cb5-31" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:0.7598940 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 </span></span>
<span id="cb5-32"><a href="#cb5-32" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :0.9997130 Max. :1.0000 Max. :1.0000 Max. :1.0000 </span></span>
<span id="cb5-33"><a href="#cb5-33" aria-hidden="true" tabindex="-1"></a><span class="do">## NA's :773 NA's :773 NA's :773 </span></span>
<span id="cb5-34"><a href="#cb5-34" aria-hidden="true" tabindex="-1"></a><span class="do">## MQSB SGB MQ0F ICB </span></span>
<span id="cb5-35"><a href="#cb5-35" aria-hidden="true" tabindex="-1"></a><span class="do">## Min. :0.01348 Min. :-0.6931 Min. :0.00000 Mode:logical </span></span>
<span id="cb5-36"><a href="#cb5-36" aria-hidden="true" tabindex="-1"></a><span class="do">## 1st Qu.:0.95494 1st Qu.:-0.6762 1st Qu.:0.00000 NA's:801 </span></span>
<span id="cb5-37"><a href="#cb5-37" aria-hidden="true" tabindex="-1"></a><span class="do">## Median :1.00000 Median :-0.6620 Median :0.00000 </span></span>
<span id="cb5-38"><a href="#cb5-38" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean :0.96428 Mean :-0.6444 Mean :0.01127 </span></span>
<span id="cb5-39"><a href="#cb5-39" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:1.00000 3rd Qu.:-0.6364 3rd Qu.:0.00000 </span></span>
<span id="cb5-40"><a href="#cb5-40" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :1.01283 Max. :-0.4536 Max. :0.66667 </span></span>
<span id="cb5-41"><a href="#cb5-41" aria-hidden="true" tabindex="-1"></a><span class="do">## NA's :48 </span></span>
<span id="cb5-42"><a href="#cb5-42" aria-hidden="true" tabindex="-1"></a><span class="do">## HOB AC AN DP4 MQ </span></span>
<span id="cb5-43"><a href="#cb5-43" aria-hidden="true" tabindex="-1"></a><span class="do">## Mode:logical Min. :1 Min. :1 Length:801 Min. :10.00 </span></span>
<span id="cb5-44"><a href="#cb5-44" aria-hidden="true" tabindex="-1"></a><span class="do">## NA's:801 1st Qu.:1 1st Qu.:1 Class :character 1st Qu.:60.00 </span></span>
<span id="cb5-45"><a href="#cb5-45" aria-hidden="true" tabindex="-1"></a><span class="do">## Median :1 Median :1 Mode :character Median :60.00 </span></span>
<span id="cb5-46"><a href="#cb5-46" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean :1 Mean :1 Mean :58.19 </span></span>
<span id="cb5-47"><a href="#cb5-47" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:1 3rd Qu.:1 3rd Qu.:60.00 </span></span>
<span id="cb5-48"><a href="#cb5-48" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :1 Max. :1 Max. :60.00 </span></span>
<span id="cb5-49"><a href="#cb5-49" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb5-50"><a href="#cb5-50" aria-hidden="true" tabindex="-1"></a><span class="do">## Indiv gt_PL gt_GT gt_GT_alleles </span></span>
<span id="cb5-51"><a href="#cb5-51" aria-hidden="true" tabindex="-1"></a><span class="do">## Length:801 Length:801 Min. :1 Length:801 </span></span>
<span id="cb5-52"><a href="#cb5-52" aria-hidden="true" tabindex="-1"></a><span class="do">## Class :character Class :character 1st Qu.:1 Class :character </span></span>
<span id="cb5-53"><a href="#cb5-53" aria-hidden="true" tabindex="-1"></a><span class="do">## Mode :character Mode :character Median :1 Mode :character </span></span>
<span id="cb5-54"><a href="#cb5-54" aria-hidden="true" tabindex="-1"></a><span class="do">## Mean :1 </span></span>
<span id="cb5-55"><a href="#cb5-55" aria-hidden="true" tabindex="-1"></a><span class="do">## 3rd Qu.:1 </span></span>
<span id="cb5-56"><a href="#cb5-56" aria-hidden="true" tabindex="-1"></a><span class="do">## Max. :1 </span></span>
<span id="cb5-57"><a href="#cb5-57" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
<div class="cell">
<details>
<summary>Code</summary>
<div class="sourceCode cell-code" id="cb6"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a><span class="fu">str</span>(variants)</span>
<span id="cb6-2"><a href="#cb6-2" aria-hidden="true" tabindex="-1"></a><span class="do">## 'data.frame': 801 obs. of 29 variables:</span></span>
<span id="cb6-3"><a href="#cb6-3" aria-hidden="true" tabindex="-1"></a><span class="do">## $ sample_id : chr "SRR2584863" "SRR2584863" "SRR2584863" "SRR2584863" ...</span></span>
<span id="cb6-4"><a href="#cb6-4" aria-hidden="true" tabindex="-1"></a><span class="do">## $ CHROM : chr "CP000819.1" "CP000819.1" "CP000819.1" "CP000819.1" ...</span></span>
<span id="cb6-5"><a href="#cb6-5" aria-hidden="true" tabindex="-1"></a><span class="do">## $ POS : int 9972 263235 281923 433359 473901 648692 1331794 1733343 2103887 2333538 ...</span></span>
<span id="cb6-6"><a href="#cb6-6" aria-hidden="true" tabindex="-1"></a><span class="do">## $ ID : logi NA NA NA NA NA NA ...</span></span>
<span id="cb6-7"><a href="#cb6-7" aria-hidden="true" tabindex="-1"></a><span class="do">## $ REF : chr "T" "G" "G" "CTTTTTTT" ...</span></span>
<span id="cb6-8"><a href="#cb6-8" aria-hidden="true" tabindex="-1"></a><span class="do">## $ ALT : chr "G" "T" "T" "CTTTTTTTT" ...</span></span>
<span id="cb6-9"><a href="#cb6-9" aria-hidden="true" tabindex="-1"></a><span class="do">## $ QUAL : num 91 85 217 64 228 210 178 225 56 167 ...</span></span>
<span id="cb6-10"><a href="#cb6-10" aria-hidden="true" tabindex="-1"></a><span class="do">## $ FILTER : logi NA NA NA NA NA NA ...</span></span>
<span id="cb6-11"><a href="#cb6-11" aria-hidden="true" tabindex="-1"></a><span class="do">## $ INDEL : logi FALSE FALSE FALSE TRUE TRUE FALSE ...</span></span>
<span id="cb6-12"><a href="#cb6-12" aria-hidden="true" tabindex="-1"></a><span class="do">## $ IDV : int NA NA NA 12 9 NA NA NA 2 7 ...</span></span>
<span id="cb6-13"><a href="#cb6-13" aria-hidden="true" tabindex="-1"></a><span class="do">## $ IMF : num NA NA NA 1 0.9 ...</span></span>
<span id="cb6-14"><a href="#cb6-14" aria-hidden="true" tabindex="-1"></a><span class="do">## $ DP : int 4 6 10 12 10 10 8 11 3 7 ...</span></span>
<span id="cb6-15"><a href="#cb6-15" aria-hidden="true" tabindex="-1"></a><span class="do">## $ VDB : num 0.0257 0.0961 0.7741 0.4777 0.6595 ...</span></span>
<span id="cb6-16"><a href="#cb6-16" aria-hidden="true" tabindex="-1"></a><span class="do">## $ RPB : num NA 1 NA NA NA NA NA NA NA NA ...</span></span>
<span id="cb6-17"><a href="#cb6-17" aria-hidden="true" tabindex="-1"></a><span class="do">## $ MQB : num NA 1 NA NA NA NA NA NA NA NA ...</span></span>
<span id="cb6-18"><a href="#cb6-18" aria-hidden="true" tabindex="-1"></a><span class="do">## $ BQB : num NA 1 NA NA NA NA NA NA NA NA ...</span></span>
<span id="cb6-19"><a href="#cb6-19" aria-hidden="true" tabindex="-1"></a><span class="do">## $ MQSB : num NA NA 0.975 1 0.916 ...</span></span>
<span id="cb6-20"><a href="#cb6-20" aria-hidden="true" tabindex="-1"></a><span class="do">## $ SGB : num -0.556 -0.591 -0.662 -0.676 -0.662 ...</span></span>
<span id="cb6-21"><a href="#cb6-21" aria-hidden="true" tabindex="-1"></a><span class="do">## $ MQ0F : num 0 0.167 0 0 0 ...</span></span>
<span id="cb6-22"><a href="#cb6-22" aria-hidden="true" tabindex="-1"></a><span class="do">## $ ICB : logi NA NA NA NA NA NA ...</span></span>
<span id="cb6-23"><a href="#cb6-23" aria-hidden="true" tabindex="-1"></a><span class="do">## $ HOB : logi NA NA NA NA NA NA ...</span></span>
<span id="cb6-24"><a href="#cb6-24" aria-hidden="true" tabindex="-1"></a><span class="do">## $ AC : int 1 1 1 1 1 1 1 1 1 1 ...</span></span>
<span id="cb6-25"><a href="#cb6-25" aria-hidden="true" tabindex="-1"></a><span class="do">## $ AN : int 1 1 1 1 1 1 1 1 1 1 ...</span></span>
<span id="cb6-26"><a href="#cb6-26" aria-hidden="true" tabindex="-1"></a><span class="do">## $ DP4 : chr "0,0,0,4" "0,1,0,5" "0,0,4,5" "0,1,3,8" ...</span></span>
<span id="cb6-27"><a href="#cb6-27" aria-hidden="true" tabindex="-1"></a><span class="do">## $ MQ : int 60 33 60 60 60 60 60 60 60 60 ...</span></span>
<span id="cb6-28"><a href="#cb6-28" aria-hidden="true" tabindex="-1"></a><span class="do">## $ Indiv : chr "/home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam" "/home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam" "/home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam" "/home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam" ...</span></span>
<span id="cb6-29"><a href="#cb6-29" aria-hidden="true" tabindex="-1"></a><span class="do">## $ gt_PL : chr "121,0" "112,0" "247,0" "91,0" ...</span></span>
<span id="cb6-30"><a href="#cb6-30" aria-hidden="true" tabindex="-1"></a><span class="do">## $ gt_GT : int 1 1 1 1 1 1 1 1 1 1 ...</span></span>
<span id="cb6-31"><a href="#cb6-31" aria-hidden="true" tabindex="-1"></a><span class="do">## $ gt_GT_alleles: chr "G" "T" "T" "CTTTTTTTT" ...</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
<p>Let’s walk through some of the important columns:</p>
<section id="sample_id" class="level3">
<h3 class="anchored" data-anchor-id="sample_id"><strong><code>sample_id</code></strong></h3>
<p>The data should include only three samples:</p>
<table class="table">
<colgroup>
<col style="width: 17%">
<col style="width: 10%">
<col style="width: 13%">
<col style="width: 9%">
<col style="width: 15%">
<col style="width: 14%">
<col style="width: 19%">
</colgroup>
<thead>
<tr class="header">
<th style="text-align: left;">SRA Run Number</th>
<th style="text-align: left;">Clone</th>
<th style="text-align: left;">Generation</th>
<th style="text-align: left;">Cit</th>
<th style="text-align: left;">Hypermutable</th>
<th style="text-align: left;">Read Length</th>
<th style="text-align: left;">Sequencing Depth</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td style="text-align: left;">SRR2589044</td>
<td style="text-align: left;">REL2181A</td>
<td style="text-align: left;">5,000</td>
<td style="text-align: left;">Unknown</td>
<td style="text-align: left;">None</td>
<td style="text-align: left;">150</td>
<td style="text-align: left;">60.2</td>
</tr>
<tr class="even">
<td style="text-align: left;">SRR2584863</td>
<td style="text-align: left;">REL7179B</td>
<td style="text-align: left;">15,000</td>
<td style="text-align: left;">Unknown</td>
<td style="text-align: left;">None</td>
<td style="text-align: left;">150</td>
<td style="text-align: left;">88</td>
</tr>
<tr class="odd">
<td style="text-align: left;">SRR2584866</td>
<td style="text-align: left;">REL11365</td>
<td style="text-align: left;">50,000</td>
<td style="text-align: left;">Cit+</td>
<td style="text-align: left;">plus</td>
<td style="text-align: left;">150</td>
<td style="text-align: left;">138.3</td>
</tr>
</tbody>
</table>
<p>Let’s check if that is true.</p>
<p>Write a code to show the unique entries in <strong><code>sample_id</code></strong></p>
<div class="cell">
<details>
<summary>Code</summary>
<div class="sourceCode cell-code" id="cb7"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a><span class="fu">unique</span>(variants<span class="sc">$</span>sample_id)</span>
<span id="cb7-2"><a href="#cb7-2" aria-hidden="true" tabindex="-1"></a><span class="do">## [1] "SRR2584863" "SRR2584866" "SRR2589044"</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
</section>
<section id="chrom" class="level3">
<h3 class="anchored" data-anchor-id="chrom"><code>CHROM</code></h3>
<p>Remember that the files are generated from <em>E. coli</em>, which should have only one (circular) chromosome. Let’s confirm this in <code>CHROM</code>.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb8"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a><span class="co"># length display the lenght of a vector. </span></span>
<span id="cb8-2"><a href="#cb8-2" aria-hidden="true" tabindex="-1"></a><span class="co"># == 1 will give a logical result of whether it's T/F</span></span>
<span id="cb8-3"><a href="#cb8-3" aria-hidden="true" tabindex="-1"></a><span class="fu">length</span>(<span class="fu">unique</span>(variants<span class="sc">$</span>CHROM))<span class="sc">==</span><span class="dv">1</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code>[1] TRUE</code></pre>
</div>
</div>
</section>
<section id="pos" class="level3">
<h3 class="anchored" data-anchor-id="pos"><code>POS</code></h3>
<p><code>POS</code> is the position of where the variant occur. We can look at the count of these variant at different <code>POS</code> and separate it out by <code>sample_id.</code></p>
<div class="cell">
<div class="sourceCode cell-code" id="cb10"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(variants, <span class="fu">aes</span>(POS)) <span class="sc">+</span> </span>
<span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a> <span class="fu">geom_freqpoly</span>() <span class="sc">+</span> </span>
<span id="cb10-3"><a href="#cb10-3" aria-hidden="true" tabindex="-1"></a> <span class="fu">facet_grid</span>(<span class="fu">factor</span>(sample_id, <span class="at">levels =</span> <span class="fu">c</span>(<span class="st">"SRR2589044"</span>, <span class="st">"SRR2584863"</span>, <span class="st">"SRR2584866"</span>))<span class="sc">~</span>.) <span class="co"># The levels for sample_id are adjusted within ggplot. Previously we did it by modifying the dataframe.</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output-display">
<p><img src="R_vcf_files/figure-html/unnamed-chunk-8-1.png" class="img-fluid" width="672"></p>
</div>
</div>
</section>
<section id="refalt" class="level3">
<h3 class="anchored" data-anchor-id="refalt"><code>REF/ALT</code></h3>
<p><code>REF</code> is the nucleotide in the reference genome.</p>
<p><code>ALT</code> is the nucleotide in the sample genome.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb11"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a><span class="fu">unique</span>(variants<span class="sc">$</span>REF)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code> [1] "T" "G"
[3] "CTTTTTTT" "CCGC"
[5] "C" "ACAGCCAGCCAGCCAGCCAGCCAGCCAGCCAG"
[7] "AT" "A"
[9] "TGG" "AGGGG"
[11] "GTTTTTTTTT" "CAA"
[13] "GAA" "CTTTTTTTT"
[15] "CGGGGG" "CTTTTT"
[17] "TGGGGGGG" "AGGGGGGGG"
[19] "TGGGGGG" "CAAAAAAA"
[21] "AGGGGG" "GTTTT"
[23] "GTTTTTTT" "TGGGGG"
[25] "AGGGGGG" "ACCCCCCC"
[27] "CGGGGGG" "GAAAAAAAA"
[29] "AGGGGGGG" "GTTTTTTTT"
[31] "GCCCCCCC" "TCCCCCCC"
[33] "TCCCCCC" "TCCCCC"
[35] "GCCCCCC" "ATTTTTTTTT"
[37] "CTTTTTTTTT" "ATTTTTTTT"
[39] "ATT" "CATGATGATGATGAT"
[41] "TAAAAAA" "ACCCC"
[43] "CAAAAAAAA" "GCCCCC"
[45] "CAAAAAAAAA" "TAAAA"
[47] "GTT" "CAAAAAA"
[49] "GTTTTTT" "GCCCC"
[51] "CT" "GTTCTTCTTCTTC"
[53] "AGCGCGCGCGCG" "TAAAAAAAAA"
[55] "GTTTCGCTTTCGCT" "GA"
[57] "CCA" "CG"
[59] "AG" </code></pre>
</div>
</div>
<p>You’ll see that it includes more than just single nucleotide polymorphisms (SNPs).</p>
<p>For example, for the <code>REF</code> “CATGATGATGATGAT”, we can pull the row(s) that contains this, and display only the columns <code>REF</code> and <code>ALT</code>.</p>
<p>Turns out, the reference has four units of (GAT), but the variant has only 3.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb13"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a><span class="fu">filter</span>(variants, REF <span class="sc">==</span><span class="st">"CATGATGATGATGAT"</span>)[, <span class="fu">c</span>(<span class="st">"REF"</span>, <span class="st">"ALT"</span>)]</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code> REF ALT
1 CATGATGATGATGAT CATGATGATGAT</code></pre>
</div>
</div>
</section>
<section id="qual" class="level3">
<h3 class="anchored" data-anchor-id="qual"><code>QUAL</code></h3>
<p><code>QUAL</code> is a Phred-scaled probability that the observed variant exists at this site. Similar to that in fastq files, a higher score means a higher probability. We can look at the distribution of <code>QUAL</code> using <code>geom_histogram()</code>. Ideally we should filter the variants based on <code>QAUL</code> before analysis. But we are using them as is here.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb15"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb15-1"><a href="#cb15-1" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(variants, <span class="fu">aes</span>(QUAL)) <span class="sc">+</span> <span class="fu">geom_histogram</span>()</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output-display">
<p><img src="R_vcf_files/figure-html/unnamed-chunk-11-1.png" class="img-fluid" width="672"></p>
</div>
</div>
</section>
<section id="indel" class="level3">
<h3 class="anchored" data-anchor-id="indel"><code>INDEL</code></h3>
<p><code>INDEL</code> is a logical vector (T/F) that indicates whether the variant is a insertion/deletion.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb16"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb16-1"><a href="#cb16-1" aria-hidden="true" tabindex="-1"></a><span class="fu">class</span>(variants<span class="sc">$</span>INDEL)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code>[1] "logical"</code></pre>
</div>
<div class="sourceCode cell-code" id="cb18"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb18-1"><a href="#cb18-1" aria-hidden="true" tabindex="-1"></a><span class="fu">unique</span>(variants<span class="sc">$</span>INDEL)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code>[1] FALSE TRUE</code></pre>
</div>
</div>
<p>We can filter logical vector, just like other character vector:</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb20"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb20-1"><a href="#cb20-1" aria-hidden="true" tabindex="-1"></a><span class="fu">filter</span>(variants, INDEL<span class="sc">==</span>T) <span class="sc">%>%</span> <span class="fu">head</span>()</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code> sample_id CHROM POS ID REF
1 SRR2584863 CP000819.1 433359 NA CTTTTTTT
2 SRR2584863 CP000819.1 473901 NA CCGC
3 SRR2584863 CP000819.1 2103887 NA ACAGCCAGCCAGCCAGCCAGCCAGCCAGCCAG
4 SRR2584863 CP000819.1 2333538 NA AT
5 SRR2584863 CP000819.1 3901455 NA A
6 SRR2584863 CP000819.1 4431393 NA TGG
ALT QUAL FILTER
1 CTTTTTTTT 64.0000 NA
2 CCGCGC 228.0000 NA
3 ACAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAG 56.0000 NA
4 ATT 167.0000 NA
5 AC 43.4147 NA
6 T 228.0000 NA
INDEL IDV IMF DP VDB RPB MQB BQB MQSB SGB MQ0F ICB HOB AC
1 TRUE 12 1.000000 12 0.477704 NA NA NA 1.000000 -0.676189 0 NA NA 1
2 TRUE 9 0.900000 10 0.659505 NA NA NA 0.916482 -0.662043 0 NA NA 1
3 TRUE 2 0.666667 3 0.901652 NA NA NA 1.000000 -0.453602 0 NA NA 1
4 TRUE 7 1.000000 7 0.568173 NA NA NA 1.012830 -0.616816 0 NA NA 1
5 TRUE 2 1.000000 2 0.020000 NA NA NA NA -0.453602 0 NA NA 1
6 TRUE 10 1.000000 10 0.863210 NA NA NA 1.007750 -0.662043 0 NA NA 1
AN DP4 MQ
1 1 0,1,3,8 60
2 1 1,0,2,7 60
3 1 0,1,1,1 60
4 1 0,1,3,3 60
5 1 0,0,2,0 60
6 1 0,1,6,3 60
Indiv gt_PL
1 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 91,0
2 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 255,0
3 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 111,28
4 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 194,0
5 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 73,0
6 /home/dcuser/dc_workshop/results/bam/SRR2584863.aligned.sorted.bam 255,0
gt_GT gt_GT_alleles
1 1 CTTTTTTTT
2 1 CCGCGC
3 1 ACAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAGCCAG
4 1 ATT
5 1 AC
6 1 T</code></pre>
</div>
</div>
<p>We can modify the plot before to show the INDEL in different color. Note that the bars are stacked instead of overlapping.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb22"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb22-1"><a href="#cb22-1" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(variants, <span class="fu">aes</span>(POS, <span class="at">fill=</span>INDEL)) <span class="sc">+</span> </span>
<span id="cb22-2"><a href="#cb22-2" aria-hidden="true" tabindex="-1"></a> <span class="fu">geom_histogram</span>() <span class="sc">+</span> </span>
<span id="cb22-3"><a href="#cb22-3" aria-hidden="true" tabindex="-1"></a> <span class="fu">facet_grid</span>(<span class="fu">factor</span>(sample_id, <span class="at">levels =</span> <span class="fu">c</span>(<span class="st">"SRR2589044"</span>, <span class="st">"SRR2584863"</span>, <span class="st">"SRR2584866"</span>))<span class="sc">~</span>.) </span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output-display">
<p><img src="R_vcf_files/figure-html/unnamed-chunk-14-1.png" class="img-fluid" width="672"></p>
</div>
</div>
</section>
</section>
<section id="analysis" class="level1" data-number="4">
<h1 data-number="4"><span class="header-section-number">4</span> Analysis</h1>
<p>Remember that the sample_id represent <em>E. coli</em> colonies sampled from different generations, with different characters (e.g., Cit+ and hypermutable).</p>
<section id="lets-test-whether-the-number-of-genetic-variants-increase-with-generation." class="level3">
<h3 class="anchored" data-anchor-id="lets-test-whether-the-number-of-genetic-variants-increase-with-generation.">Let’s test whether the number of genetic variants increase with generation.</h3>
<p>We can make a dataframe to contain these information:</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb23"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb23-1"><a href="#cb23-1" aria-hidden="true" tabindex="-1"></a>exp_detail <span class="ot">=</span> <span class="fu">data.frame</span>(</span>
<span id="cb23-2"><a href="#cb23-2" aria-hidden="true" tabindex="-1"></a> <span class="at">sample_id =</span> <span class="fu">c</span>(<span class="st">"SRR2589044"</span>, <span class="st">"SRR2584863"</span>, <span class="st">"SRR2584866"</span>), </span>
<span id="cb23-3"><a href="#cb23-3" aria-hidden="true" tabindex="-1"></a> <span class="at">generation =</span> <span class="fu">c</span>(<span class="dv">5000</span>, <span class="dv">15000</span>, <span class="dv">50000</span>), </span>
<span id="cb23-4"><a href="#cb23-4" aria-hidden="true" tabindex="-1"></a> <span class="at">Cit =</span> <span class="fu">c</span>(<span class="st">"unknown"</span>, <span class="st">"unknown"</span>, <span class="st">"Cit+"</span>), </span>
<span id="cb23-5"><a href="#cb23-5" aria-hidden="true" tabindex="-1"></a> <span class="at">hypermutable =</span> <span class="fu">c</span>(<span class="st">"none"</span>, <span class="st">"none"</span>, <span class="st">"plus"</span>),</span>
<span id="cb23-6"><a href="#cb23-6" aria-hidden="true" tabindex="-1"></a> <span class="at">seq_depth =</span> <span class="fu">c</span>(<span class="fl">60.2</span>, <span class="dv">88</span>, <span class="fl">138.3</span>)</span>
<span id="cb23-7"><a href="#cb23-7" aria-hidden="true" tabindex="-1"></a> )</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Let’s count the number of variants in each <code>sample_id</code></p>
<p>We know that each row represent one variants. <code>table()</code> allows us to count the occurrence of each entry in a character vector.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb24"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb24-1"><a href="#cb24-1" aria-hidden="true" tabindex="-1"></a><span class="fu">table</span>(variants<span class="sc">$</span>sample_id)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code>
SRR2584863 SRR2584866 SRR2589044
25 766 10 </code></pre>
</div>
</div>
<p>We can convert this into a dataframe and combine it into exp_detail.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb26"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb26-1"><a href="#cb26-1" aria-hidden="true" tabindex="-1"></a>var_summary <span class="ot">=</span> <span class="fu">as.data.frame</span>(<span class="fu">table</span>(variants<span class="sc">$</span>sample_id))</span>
<span id="cb26-2"><a href="#cb26-2" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb26-3"><a href="#cb26-3" aria-hidden="true" tabindex="-1"></a><span class="co"># fix names</span></span>
<span id="cb26-4"><a href="#cb26-4" aria-hidden="true" tabindex="-1"></a><span class="fu">names</span>(var_summary) <span class="ot">=</span> <span class="fu">c</span>(<span class="st">"sample_id"</span>, <span class="st">"variant_count"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Now use <code>left_join()</code> to combine. (right_join will also work as we have the same sample_id in the two tables.)</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb27"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb27-1"><a href="#cb27-1" aria-hidden="true" tabindex="-1"></a><span class="co"># We'll save it as a new object</span></span>
<span id="cb27-2"><a href="#cb27-2" aria-hidden="true" tabindex="-1"></a>exp_result <span class="ot">=</span> <span class="fu">left_join</span>(exp_detail, var_summary, <span class="st">"sample_id"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>We can then plot and test whether the <code>variant_count</code> increases with <code>generation</code>.</p>
<p>Write a code to do the regression analysis with <code>log10(variant_count)</code>:</p>
<blockquote class="blockquote">
<p>Note that in the assignment before, we created a new column that has the log10 of a variable, then we run the test with that new column. An easier alternative is to directly use <code>log10(variant_count)</code> in the formula of <code>lm()</code> and in <code>ggplot</code>.</p>
</blockquote>
<div class="cell">
<details>
<summary>Code</summary>
<div class="sourceCode cell-code" id="cb28"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb28-1"><a href="#cb28-1" aria-hidden="true" tabindex="-1"></a><span class="fu">summary</span>(<span class="fu">lm</span>(<span class="fu">log10</span>(variant_count)<span class="sc">~</span>generation, exp_result))</span>
<span id="cb28-2"><a href="#cb28-2" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb28-3"><a href="#cb28-3" aria-hidden="true" tabindex="-1"></a><span class="do">## Call:</span></span>
<span id="cb28-4"><a href="#cb28-4" aria-hidden="true" tabindex="-1"></a><span class="do">## lm(formula = log10(variant_count) ~ generation, data = exp_result)</span></span>
<span id="cb28-5"><a href="#cb28-5" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb28-6"><a href="#cb28-6" aria-hidden="true" tabindex="-1"></a><span class="do">## Residuals:</span></span>
<span id="cb28-7"><a href="#cb28-7" aria-hidden="true" tabindex="-1"></a><span class="do">## 1 2 3 </span></span>
<span id="cb28-8"><a href="#cb28-8" aria-hidden="true" tabindex="-1"></a><span class="do">## 0.009769 -0.012560 0.002791 </span></span>
<span id="cb28-9"><a href="#cb28-9" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb28-10"><a href="#cb28-10" aria-hidden="true" tabindex="-1"></a><span class="do">## Coefficients:</span></span>
<span id="cb28-11"><a href="#cb28-11" aria-hidden="true" tabindex="-1"></a><span class="do">## Estimate Std. Error t value Pr(>|t|) </span></span>
<span id="cb28-12"><a href="#cb28-12" aria-hidden="true" tabindex="-1"></a><span class="do">## (Intercept) 7.801e-01 1.464e-02 53.30 0.01194 * </span></span>
<span id="cb28-13"><a href="#cb28-13" aria-hidden="true" tabindex="-1"></a><span class="do">## generation 4.203e-05 4.834e-07 86.94 0.00732 **</span></span>
<span id="cb28-14"><a href="#cb28-14" aria-hidden="true" tabindex="-1"></a><span class="do">## ---</span></span>
<span id="cb28-15"><a href="#cb28-15" aria-hidden="true" tabindex="-1"></a><span class="do">## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</span></span>
<span id="cb28-16"><a href="#cb28-16" aria-hidden="true" tabindex="-1"></a><span class="do">## </span></span>
<span id="cb28-17"><a href="#cb28-17" aria-hidden="true" tabindex="-1"></a><span class="do">## Residual standard error: 0.01615 on 1 degrees of freedom</span></span>
<span id="cb28-18"><a href="#cb28-18" aria-hidden="true" tabindex="-1"></a><span class="do">## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9997 </span></span>
<span id="cb28-19"><a href="#cb28-19" aria-hidden="true" tabindex="-1"></a><span class="do">## F-statistic: 7558 on 1 and 1 DF, p-value: 0.007322</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
<p>Write a code to plot <code>log10(variant_count)</code> against generations.</p>
<div class="cell">
<details>
<summary>Code</summary>
<div class="sourceCode cell-code" id="cb29"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb29-1"><a href="#cb29-1" aria-hidden="true" tabindex="-1"></a><span class="fu">ggplot</span>(exp_result, <span class="fu">aes</span>(<span class="at">x =</span> generation, <span class="at">y =</span> <span class="fu">log10</span>(variant_count))) <span class="sc">+</span> </span>
<span id="cb29-2"><a href="#cb29-2" aria-hidden="true" tabindex="-1"></a> <span class="fu">geom_point</span>() <span class="sc">+</span></span>
<span id="cb29-3"><a href="#cb29-3" aria-hidden="true" tabindex="-1"></a> <span class="fu">geom_smooth</span>(<span class="at">method =</span> lm) <span class="sc">+</span></span>
<span id="cb29-4"><a href="#cb29-4" aria-hidden="true" tabindex="-1"></a> <span class="fu">labs</span>(<span class="at">x =</span> <span class="st">"Generations"</span>, <span class="at">y =</span> <span class="st">"Log count of genetic variant"</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</details>
</div>
</section>
</section>
<section id="assignment" class="level1" data-number="5">
<h1 data-number="5"><span class="header-section-number">5</span> Assignment</h1>
<p>We know that there are many types of genetic variants in this VCF files. We have shown that the total number of variants increased with the number of generation. Here, the goal is to test:</p>
<blockquote class="blockquote">
<p><strong>Whether the number of SNPs increased with the number of generation.</strong></p>
</blockquote>
<p>Write codes to do the following and submit the R console commands & outputs in the assignment. Remember to annotate your codes.</p>
<ol type="1">
<li><p>Since we are only interested in SNPs, that means that <code>REF</code> should be either A, T, G, or C. Write some codes to <code>filter</code> <code>variant</code> to have only A/T/C/G in <code>REF</code>. Save it as a new object called <code>SNP</code></p></li>
<li><p>In the new object <code>SNP</code>, each row represent a SNP. Again there are still three unique <code>sample_id</code>, each representing a sample. Write some codes to generate a dataframe that shows the count of rows for each sample. Call this dataframe <code>SNP_summary</code>. It should have two columns, one named <code>sample_id</code>, one named <code>SNP_count</code>.</p></li>
<li><p>Now combine <code>SNP_summary</code> with <code>exp_result</code> based on <code>sample_id</code>. Save it as a new object called <code>exp_result2</code>.</p></li>
<li><p>Do a linear regression to test whether the <strong>log10</strong> of the number of SNP increases with generation. Write a statement to summarize the results and statistics.</p></li>
<li><p>Plot the data to visualize the regression result.</p></li>
</ol>
<div class="cell">
</div>
</section>
</main>
<!-- /main column -->
<script id="quarto-html-after-body" type="application/javascript">
window.document.addEventListener("DOMContentLoaded", function (event) {
const toggleBodyColorMode = (bsSheetEl) => {
const mode = bsSheetEl.getAttribute("data-mode");
const bodyEl = window.document.querySelector("body");
if (mode === "dark") {
bodyEl.classList.add("quarto-dark");
bodyEl.classList.remove("quarto-light");
} else {
bodyEl.classList.add("quarto-light");
bodyEl.classList.remove("quarto-dark");
}
}
const toggleBodyColorPrimary = () => {
const bsSheetEl = window.document.querySelector("link#quarto-bootstrap");
if (bsSheetEl) {
toggleBodyColorMode(bsSheetEl);
}
}
toggleBodyColorPrimary();
const icon = "";
const anchorJS = new window.AnchorJS();
anchorJS.options = {
placement: 'right',
icon: icon
};
anchorJS.add('.anchored');
const clipboard = new window.ClipboardJS('.code-copy-button', {
target: function(trigger) {
return trigger.previousElementSibling;
}
});
clipboard.on('success', function(e) {
// button target
const button = e.trigger;
// don't keep focus
button.blur();
// flash "checked"
button.classList.add('code-copy-button-checked');
var currentTitle = button.getAttribute("title");
button.setAttribute("title", "Copied!");
let tooltip;
if (window.bootstrap) {
button.setAttribute("data-bs-toggle", "tooltip");
button.setAttribute("data-bs-placement", "left");
button.setAttribute("data-bs-title", "Copied!");
tooltip = new bootstrap.Tooltip(button,
{ trigger: "manual",
customClass: "code-copy-button-tooltip",
offset: [0, -8]});
tooltip.show();
}
setTimeout(function() {
if (tooltip) {
tooltip.hide();
button.removeAttribute("data-bs-title");
button.removeAttribute("data-bs-toggle");
button.removeAttribute("data-bs-placement");
}
button.setAttribute("title", currentTitle);
button.classList.remove('code-copy-button-checked');
}, 1000);
// clear code selection
e.clearSelection();
});
function tippyHover(el, contentFn) {
const config = {
allowHTML: true,
content: contentFn,
maxWidth: 500,
delay: 100,
arrow: false,
appendTo: function(el) {
return el.parentElement;
},
interactive: true,
interactiveBorder: 10,
theme: 'quarto',
placement: 'bottom-start'
};
window.tippy(el, config);
}
const noterefs = window.document.querySelectorAll('a[role="doc-noteref"]');
for (var i=0; i<noterefs.length; i++) {
const ref = noterefs[i];
tippyHover(ref, function() {
// use id or data attribute instead here
let href = ref.getAttribute('data-footnote-href') || ref.getAttribute('href');
try { href = new URL(href).hash; } catch {}
const id = href.replace(/^#\/?/, "");
const note = window.document.getElementById(id);
return note.innerHTML;
});
}
const findCites = (el) => {
const parentEl = el.parentElement;
if (parentEl) {
const cites = parentEl.dataset.cites;
if (cites) {
return {
el,
cites: cites.split(' ')
};
} else {
return findCites(el.parentElement)
}
} else {
return undefined;
}
};
var bibliorefs = window.document.querySelectorAll('a[role="doc-biblioref"]');
for (var i=0; i<bibliorefs.length; i++) {
const ref = bibliorefs[i];
const citeInfo = findCites(ref);
if (citeInfo) {
tippyHover(citeInfo.el, function() {
var popup = window.document.createElement('div');
citeInfo.cites.forEach(function(cite) {
var citeDiv = window.document.createElement('div');
citeDiv.classList.add('hanging-indent');
citeDiv.classList.add('csl-entry');
var biblioDiv = window.document.getElementById('ref-' + cite);
if (biblioDiv) {
citeDiv.innerHTML = biblioDiv.innerHTML;
}
popup.appendChild(citeDiv);
});
return popup.innerHTML;
});
}
}
});
</script>
</div> <!-- /content -->
</body></html>