2
|
1 ---
|
15
|
2 title: 'Short reads evaluation with [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)'
|
14
|
3 output:
|
|
4 html_document:
|
|
5 number_sections: true
|
|
6 toc: true
|
|
7 theme: cosmo
|
|
8 highlight: tango
|
2
|
9 ---
|
|
10
|
14
|
11 ```{r setup, include=FALSE, warning=FALSE, message=FALSE}
|
|
12 knitr::opts_chunk$set(
|
16
|
13 echo = ECHO,
|
|
14 error = TRUE
|
14
|
15 )
|
2
|
16 ```
|
|
17
|
|
18
|
16
|
19 # Fastqc Evaluation
|
14
|
20
|
16
|
21 ## Evaluation of reads before trimming
|
2
|
22
|
16
|
23 ```{r}
|
|
24 if ('READS_1' == 'None') {
|
|
25 stop("No pre-trimming reads provided!")
|
|
26 } else {
|
|
27 ## run fastqc evaluation
|
|
28 fastqc_command = paste0('fastqc ') %>%
|
|
29 (function(x) {
|
|
30 ifelse('CONTAMINANTS' != 'None', paste0(x, '-c CONTAMINANTS '), x)
|
|
31 }) %>%
|
|
32 (function(x) {
|
|
33 ifelse('LIMITS' != 'None', paste0(x, '-l LIMITS '), x)
|
|
34 }) %>%
|
|
35 (function(x) {
|
|
36 paste0(x, '-o REPORT_DIR ')
|
|
37 })
|
|
38 fastqc_command_reads_1 = paste0(fastqc_command, 'READS_1 > /dev/null 2>&1')
|
|
39 system(fastqc_command_reads_1, intern = TRUE)
|
|
40
|
|
41 # Original html report
|
|
42 reads_1_base = tail(strsplit('READS_1', '/')[[1]], 1)
|
|
43 original_html = tags$a(href=paste0(reads_1_base, '_fastqc.html'), paste0('HTML report: ', opt$name_1))
|
|
44
|
|
45 unzip(paste0('REPORT_DIR/', reads_1_base, '_fastqc.zip'), exdir = 'REPORT_DIR')
|
|
46 reads_1_unzip = paste0('REPORT_DIR/', reads_1_base, '_fastqc/')
|
|
47 # fastqc_data.txt
|
|
48 file.copy(paste0(reads_1_unzip, 'fastqc_data.txt'), 'REPORT_DIR/reads_1_fastqc_data.txt')
|
|
49 fastqc_data = tags$a(href='reads_1_fastqc_data.txt', paste0('fastqc_data.txt: ', opt$name_1))
|
|
50 # summary.txt
|
|
51 file.copy(paste0(reads_1_unzip, 'summary.txt'), 'REPORT_DIR/reads_1_summary.txt')
|
|
52 summary_data = tags$a(href='reads_1_summary.txt', paste0('summary.txt: ', opt$name_1))
|
|
53
|
|
54 tags$ul(
|
|
55 tags$li(original_html),
|
|
56 tags$li(fastqc_data),
|
|
57 tags$li(summary_data)
|
|
58 )
|
|
59 }
|
15
|
60 ```
|
|
61
|
|
62
|
16
|
63 ## Evaluation of reads after trimming
|
15
|
64
|
16
|
65 ```{r}
|
|
66 if ('READS_2' == 'None') {
|
|
67 stop("No pre-trimming reads provided!")
|
|
68 } else {
|
|
69 ## run fastqc evaluation
|
|
70 fastqc_command = paste0('fastqc ') %>%
|
|
71 (function(x) {
|
|
72 ifelse('CONTAMINANTS' != 'None', paste0(x, '-c CONTAMINANTS '), x)
|
|
73 }) %>%
|
|
74 (function(x) {
|
|
75 ifelse('LIMITS' != 'None', paste0(x, '-l LIMITS '), x)
|
|
76 }) %>%
|
|
77 (function(x) {
|
|
78 paste0(x, '-o REPORT_DIR ')
|
|
79 })
|
|
80 fastqc_command_reads_2 = paste0(fastqc_command, 'READS_2 > /dev/null 2>&1')
|
|
81 system(fastqc_command_reads_2, intern = TRUE)
|
|
82
|
|
83 # Original html report
|
|
84 reads_2_base = tail(strsplit('READS_2', '/')[[1]], 1)
|
|
85 original_html = tags$a(href=paste0(reads_2_base, '_fastqc.html'), paste0('HTML report: ', opt$name_2))
|
|
86
|
|
87 unzip(paste0('REPORT_DIR/', reads_2_base, '_fastqc.zip'), exdir = 'REPORT_DIR')
|
|
88 reads_2_unzip = paste0('REPORT_DIR/', reads_2_base, '_fastqc/')
|
|
89 # fastqc_data.txt
|
|
90 file.copy(paste0(reads_2_unzip, 'fastqc_data.txt'), 'REPORT_DIR/reads_2_fastqc_data.txt')
|
|
91 fastqc_data = tags$a(href='reads_2_fastqc_data.txt', paste0('fastqc_data.txt: ', opt$name_2))
|
|
92 # summary.txt
|
|
93 file.copy(paste0(reads_2_unzip, 'summary.txt'), 'REPORT_DIR/reads_2_summary.txt')
|
|
94 summary_data = tags$a(href='reads_2_summary.txt', paste0('summary.txt: ', opt$name_2))
|
|
95
|
|
96 tags$ul(
|
|
97 tags$li(original_html),
|
|
98 tags$li(fastqc_data),
|
|
99 tags$li(summary_data)
|
|
100 )
|
|
101 }
|
2
|
102 ```
|
|
103
|
15
|
104
|
|
105
|
|
106 # Fastqc output visualization
|
|
107
|
|
108 ## Overview
|
|
109
|
|
110 ```{r}
|
16
|
111 reads_1_summary = read.csv('REPORT_DIR/reads_1_summary.txt', header = FALSE, sep = '\t')[, 2:1]
|
17
|
112 reads_2_summary = read.csv('REPORT_DIR/reads_2_summary.txt', header = FALSE, sep = '\t')[, 1]
|
16
|
113 combined_summary = cbind(reads_1_summary, reads_2_summary)
|
|
114 names(combined_summary) = c('MODULE', paste0(opt$name_1, '(before)'), paste0(opt$name_2, '(after)'))
|
17
|
115 combined_summary[combined_summary == 'FAIL'] = 'FAIL (X)'
|
|
116 combined_summary[combined_summary == 'WARN'] = 'WARN (!)'
|
16
|
117 knitr::kable(combined_summary)
|
15
|
118 ```
|
|
119
|
17
|
120 ## Visualization by data module {.tabset}
|
2
|
121
|
14
|
122 * Define a function to extract outputs for each module from fastqc output
|
2
|
123
|
14
|
124 ```{r 'function definition'}
|
17
|
125 extract_data_module = function(fastqc_data, module_name, header = TRUE, comment.char = "") {
|
14
|
126 f = readLines(fastqc_data)
|
|
127 start_line = grep(module_name, f)
|
|
128 end_module_lines = grep('END_MODULE', f)
|
|
129 end_line = end_module_lines[which(end_module_lines > start_line)[1]]
|
|
130 module_data = f[(start_line+1):(end_line-1)]
|
|
131 writeLines(module_data, 'temp.txt')
|
17
|
132 read.csv('temp.txt', sep = '\t', header = header, comment.char = comment.char)
|
2
|
133 }
|
|
134 ```
|
|
135
|
15
|
136 ### Per base sequence quality
|
|
137
|
16
|
138 ```{r 'per base sequence quality', fig.width=10}
|
|
139 ## reads 1
|
|
140 pbsq_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Per base sequence quality')
|
|
141 pbsq_1$id = 1:length(pbsq_1$X.Base)
|
|
142
|
17
|
143 melt_pbsq_1 = filter(melt(pbsq_1, id=c('X.Base', 'id')), variable == 'Mean')
|
16
|
144 melt_pbsq_1$trim = 'before'
|
|
145
|
|
146
|
|
147 ## reads 2
|
|
148 pbsq_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per base sequence quality')
|
|
149 pbsq_2$id = 1:length(pbsq_2$X.Base)
|
|
150
|
17
|
151 melt_pbsq_2 = filter(melt(pbsq_2, id=c('X.Base', 'id')), variable == 'Mean')
|
16
|
152 melt_pbsq_2$trim = 'after'
|
|
153
|
|
154 comb_pbsq = rbind(melt_pbsq_1, melt_pbsq_2)
|
|
155 comb_pbsq$trim = factor(levels = c('before', 'after'), comb_pbsq$trim)
|
17
|
156
|
16
|
157 p = ggplot(data = comb_pbsq) +
|
|
158 geom_line(mapping = aes(x = id, y = value, group = variable, color = variable)) +
|
|
159 scale_x_continuous(breaks = pbsq_2$id, labels = pbsq_2$X.Base) +
|
|
160 facet_grid(. ~ trim) +
|
17
|
161 ylim(0, max(comb_pbsq$value) + 5) +
|
16
|
162 theme(axis.text.x = element_text(angle=45))
|
|
163 ggplotly(p)
|
|
164
|
15
|
165 ```
|
|
166
|
|
167 ### Per tile sequence quality
|
|
168
|
17
|
169 ```{r 'per tile sequence quality', fig.width=10}
|
|
170 ## check if 'per tile sequence quality' module exits or not
|
|
171 check_ptsq = grep('Per tile sequence quality', readLines('REPORT_DIR/reads_1_fastqc_data.txt'))
|
|
172 if (length(check_ptsq) > 0) {
|
|
173 ## reads 1
|
|
174 ptsq_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Per tile sequence quality')
|
|
175 ptsq_1$trim = 'before'
|
|
176
|
|
177 ## reads 2
|
|
178 ptsq_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per tile sequence quality')
|
|
179 ptsq_2$trim = 'after'
|
|
180
|
|
181 comb_ptsq = rbind(ptsq_1, ptsq_2)
|
|
182 comb_ptsq$trim = factor(levels = c('before', 'after'), comb_ptsq$trim)
|
|
183 comb_ptsq$Base = factor(levels = unique(comb_ptsq$Base), comb_ptsq$Base)
|
|
184
|
|
185 # convert integers to charaters
|
|
186 comb_ptsq$Tile = as.character(comb_ptsq$X.Tile)
|
|
187
|
|
188 p = ggplot(data = comb_ptsq, aes(x = Base, y = Tile, fill = Mean)) +
|
|
189 geom_raster() +
|
|
190 facet_grid(. ~ trim) +
|
|
191 xlab('Position in read (bp)') +
|
|
192 ylab('') +
|
|
193 theme(axis.text.x = element_text(angle=45))
|
|
194 ggplotly(p)
|
|
195 } else {
|
|
196 print('No "per tile sequence quality" data')
|
|
197 }
|
|
198
|
|
199
|
|
200 ```
|
|
201
|
|
202 ### Per sequence quality score
|
|
203
|
|
204 ```{r 'Per sequence quality score', fig.width=10}
|
16
|
205 ## reads 1
|
17
|
206 psqs_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Per sequence quality scores')
|
|
207 psqs_1$trim = 'before'
|
16
|
208
|
|
209 ## reads 2
|
17
|
210 psqs_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per sequence quality scores')
|
|
211 psqs_2$trim = 'after'
|
|
212
|
|
213 comb_psqs = rbind(psqs_1, psqs_2)
|
|
214 comb_psqs$trim = factor(levels = c('before', 'after'), comb_psqs$trim)
|
|
215
|
|
216 p = ggplot(data = comb_psqs, aes(x = X.Quality, y = Count)) +
|
|
217 geom_line(color = 'red') +
|
|
218 facet_grid(. ~ trim) +
|
|
219 xlim(min(comb_psqs$X.Quality), max(comb_psqs$X.Quality)) +
|
|
220 xlab('Mean Sequence Qaulity (Phred Score)') +
|
|
221 ylab('')
|
|
222 ggplotly(p)
|
|
223 ```
|
|
224
|
|
225
|
|
226 ### Per base sequence content
|
|
227
|
|
228 ```{r 'Per base sequence content', fig.width=10}
|
|
229 ## reads 1
|
|
230 pbsc_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Per base sequence content')
|
|
231 pbsc_1$id = 1:length(pbsc_1$X.Base)
|
|
232
|
|
233 melt_pbsc_1 = melt(pbsc_1, id=c('X.Base', 'id'))
|
|
234 melt_pbsc_1$trim = 'before'
|
|
235
|
|
236
|
|
237 ## reads 2
|
|
238 pbsc_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per base sequence content')
|
|
239 pbsc_2$id = 1:length(pbsc_2$X.Base)
|
|
240
|
|
241 melt_pbsc_2 = melt(pbsc_2, id=c('X.Base', 'id'))
|
|
242 melt_pbsc_2$trim = 'after'
|
|
243
|
|
244 comb_pbsc = rbind(melt_pbsc_1, melt_pbsc_2)
|
|
245 comb_pbsc$trim = factor(levels = c('before', 'after'), comb_pbsc$trim)
|
|
246
|
|
247 p = ggplot(data = comb_pbsc, aes(x = id, y = value, color = variable)) +
|
|
248 geom_line() +
|
|
249 facet_grid(. ~ trim) +
|
|
250 xlim(min(comb_pbsc$id), max(comb_pbsc$id)) +
|
|
251 ylim(0, 100) +
|
|
252 xlab('Position in read (bp)') +
|
|
253 ylab('')
|
|
254 ggplotly(p)
|
|
255 ```
|
16
|
256
|
17
|
257 ### Per sequence GC content
|
|
258
|
|
259 ```{r 'Per sequence GC content', fig.width=10}
|
|
260 ## reads 1
|
|
261 psGCc_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Per sequence GC content')
|
|
262 psGCc_1$trim = 'before'
|
|
263
|
|
264 ## reads 2
|
|
265 psGCc_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per sequence GC content')
|
|
266 psGCc_2$trim = 'after'
|
|
267
|
|
268 comb_psGCc = rbind(psGCc_1, psGCc_2)
|
|
269 comb_psGCc$trim = factor(levels = c('before', 'after'), comb_psGCc$trim)
|
|
270
|
|
271 p = ggplot(data = comb_psGCc, aes(x = X.GC.Content, y = Count)) +
|
|
272 geom_line(color = 'red') +
|
|
273 facet_grid(. ~ trim) +
|
|
274 xlab('Mean Sequence Qaulity (Phred Score)') +
|
|
275 ylab('')
|
|
276 ggplotly(p)
|
|
277 ```
|
|
278
|
16
|
279
|
17
|
280 ### Per base N content
|
|
281
|
|
282 ```{r 'Per base N content', fig.width=10}
|
|
283 ## reads 1
|
|
284 pbNc_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Per base N content')
|
|
285 pbNc_1$id = 1:length(pbNc_1$X.Base)
|
|
286 pbNc_1$trim = 'before'
|
|
287
|
|
288 ## reads 2
|
|
289 pbNc_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per base N content')
|
|
290 pbNc_2$id = 1:length(pbNc_2$X.Base)
|
|
291 pbNc_2$trim = 'after'
|
|
292
|
|
293 comb_pbNc = rbind(pbNc_1, pbNc_2)
|
|
294 comb_pbNc$trim = factor(levels = c('before', 'after'), comb_pbNc$trim)
|
|
295
|
|
296 p = ggplot(data = comb_pbNc, aes(x = id, y = N.Count)) +
|
|
297 geom_line(color = 'red') +
|
|
298 scale_x_continuous(breaks = pbNc_2$id, labels = pbNc_2$X.Base) +
|
|
299 facet_grid(. ~ trim) +
|
|
300 ylim(0, 1) +
|
|
301 xlab('N-Count') +
|
|
302 ylab('') +
|
16
|
303 theme(axis.text.x = element_text(angle=45))
|
|
304 ggplotly(p)
|
15
|
305 ```
|
|
306
|
|
307
|
17
|
308 ### Sequence Length Distribution
|
|
309
|
|
310 ```{r 'Sequence Length Distribution', fig.width=10}
|
|
311 ## reads 1
|
|
312 sld_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Sequence Length Distribution')
|
|
313 sld_1$id = 1:length(sld_1$X.Length)
|
|
314 sld_1$trim = 'before'
|
|
315
|
|
316 ## reads 2
|
|
317 sld_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Sequence Length Distribution')
|
|
318 sld_2$id = 1:length(sld_2$X.Length)
|
|
319 sld_2$trim = 'after'
|
|
320
|
|
321 comb_sld = rbind(sld_1, sld_2)
|
|
322 comb_sld$trim = factor(levels = c('before', 'after'), comb_sld$trim)
|
|
323
|
|
324 p = ggplot(data = comb_sld, aes(x = id, y = Count)) +
|
|
325 geom_line(color = 'red') +
|
|
326 scale_x_continuous(breaks = sld_2$id, labels = sld_2$X.Length) +
|
|
327 facet_grid(. ~ trim) +
|
|
328 xlab('Sequence Length (bp)') +
|
|
329 ylab('') +
|
|
330 theme(axis.text.x = element_text(angle=45))
|
|
331 ggplotly(p)
|
|
332 ```
|
|
333
|
|
334 ### Sequence Duplication Levels
|
|
335
|
|
336 ```{r 'Sequence Duplication Levels', fig.width=10}
|
|
337 ## reads 1
|
|
338 sdl_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Sequence Duplication Levels', header = FALSE, comment.char = '#')
|
|
339 names(sdl_1) = c('Duplication_Level', 'Percentage_of_deduplicated', 'Percentage_of_total')
|
|
340 sdl_1$id = 1:length(sdl_1$Duplication_Level)
|
|
341
|
|
342 melt_sdl_1 = melt(sdl_1, id=c('Duplication_Level', 'id'))
|
|
343 melt_sdl_1$trim = 'before'
|
|
344
|
|
345
|
|
346 ## reads 2
|
|
347 sdl_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Sequence Duplication Levels', header = FALSE, comment.char = '#')
|
|
348 names(sdl_2) = c('Duplication_Level', 'Percentage_of_deduplicated', 'Percentage_of_total')
|
|
349 sdl_2$id = 1:length(sdl_2$Duplication_Level)
|
|
350
|
|
351 melt_sdl_2 = melt(sdl_2, id=c('Duplication_Level', 'id'))
|
|
352 melt_sdl_2$trim = 'after'
|
|
353
|
|
354 comb_sdl = rbind(melt_sdl_1, melt_sdl_2)
|
|
355 comb_sdl$trim = factor(levels = c('before', 'after'), comb_sdl$trim)
|
|
356
|
|
357 p = ggplot(data = comb_sdl, aes(x = id, y = value, color = variable)) +
|
|
358 geom_line() +
|
|
359 scale_x_continuous(breaks = sdl_2$id, labels = sdl_2$Duplication_Level) +
|
|
360 facet_grid(. ~ trim) +
|
|
361 xlab('Sequence Duplication Level') +
|
|
362 ylab('') +
|
|
363 theme(axis.text.x = element_text(angle=45))
|
|
364 ggplotly(p)
|
|
365 ```
|
|
366
|
|
367 ### Adapter Content
|
|
368
|
|
369 ```{r 'Adapter Content', fig.width=10}
|
|
370 ## reads 1
|
|
371 ac_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Adapter Content')
|
|
372 ac_1$id = 1:length(ac_1$X.Position)
|
|
373
|
|
374 melt_ac_1 = melt(ac_1, id=c('X.Position', 'id'))
|
|
375 melt_ac_1$trim = 'before'
|
|
376
|
|
377 ## reads 2
|
|
378 ac_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Adapter Content')
|
|
379 ac_2$id = 1:length(ac_2$X.Position)
|
|
380
|
|
381 melt_ac_2 = melt(ac_2, id=c('X.Position', 'id'))
|
|
382 melt_ac_2$trim = 'after'
|
|
383
|
|
384 comb_ac = rbind(melt_ac_1, melt_ac_2)
|
|
385 comb_ac$trim = factor(levels = c('before', 'after'), comb_ac$trim)
|
|
386
|
|
387 p = ggplot(data = comb_ac, aes(x = id, y = value, color = variable)) +
|
|
388 geom_line() +
|
|
389 facet_grid(. ~ trim) +
|
|
390 xlim(min(comb_ac$id), max(comb_ac$id)) +
|
|
391 ylim(0, 1) +
|
|
392 xlab('Position in read (bp)') +
|
|
393 ylab('')
|
|
394 ggplotly(p)
|
|
395 ```
|
|
396
|
|
397 ### Kmer Content {.tabset}
|
|
398
|
|
399 #### Before
|
|
400
|
|
401 ```{r 'Kmer Content (before)', fig.width=10}
|
|
402 kc_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Kmer Content')
|
|
403 knitr::kable(kc_1)
|
|
404 ```
|
|
405
|
|
406 #### After
|
|
407 ```{r 'Kmer Content (after)', fig.width=10}
|
|
408 kc_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Kmer Content')
|
|
409 knitr::kable(kc_2)
|
|
410 ```
|
|
411
|
2
|
412
|
14
|
413 # Session Info
|
2
|
414
|
14
|
415 ```{r 'session info'}
|
|
416 sessionInfo()
|
2
|
417 ```
|
|
418
|