| Previous changeset 16:1710b0e874f1 (2017-10-21) Next changeset 18:8635a4cee6dd (2017-11-09) |
|
Commit message:
compare evaluation before and after trimming |
|
modified:
fastqc_report.Rmd |
| b |
| diff -r 1710b0e874f1 -r ac5c618e4d97 fastqc_report.Rmd --- a/fastqc_report.Rmd Sat Oct 21 09:25:49 2017 -0400 +++ b/fastqc_report.Rmd Mon Nov 06 16:53:14 2017 -0500 |
| [ |
| b'@@ -109,25 +109,27 @@\n \n ```{r}\n reads_1_summary = read.csv(\'REPORT_DIR/reads_1_summary.txt\', header = FALSE, sep = \'\\t\')[, 2:1]\n-reads_2_summary = read.csv(\'REPORT_DIR/reads_1_summary.txt\', header = FALSE, sep = \'\\t\')[, 1]\n+reads_2_summary = read.csv(\'REPORT_DIR/reads_2_summary.txt\', header = FALSE, sep = \'\\t\')[, 1]\n combined_summary = cbind(reads_1_summary, reads_2_summary)\n names(combined_summary) = c(\'MODULE\', paste0(opt$name_1, \'(before)\'), paste0(opt$name_2, \'(after)\'))\n+combined_summary[combined_summary == \'FAIL\'] = \'FAIL (X)\'\n+combined_summary[combined_summary == \'WARN\'] = \'WARN (!)\'\n knitr::kable(combined_summary)\n ```\n \n-## Visualization by module {.tabset}\n+## Visualization by data module {.tabset}\n \n * Define a function to extract outputs for each module from fastqc output\n \n ```{r \'function definition\'}\n-extract_data_module = function(fastqc_data, module_name) {\n+extract_data_module = function(fastqc_data, module_name, header = TRUE, comment.char = "") {\n f = readLines(fastqc_data)\n start_line = grep(module_name, f)\n end_module_lines = grep(\'END_MODULE\', f)\n end_line = end_module_lines[which(end_module_lines > start_line)[1]]\n module_data = f[(start_line+1):(end_line-1)]\n writeLines(module_data, \'temp.txt\')\n- read.csv(\'temp.txt\', sep = \'\\t\')\n+ read.csv(\'temp.txt\', sep = \'\\t\', header = header, comment.char = comment.char)\n }\n ```\n \n@@ -138,7 +140,7 @@\n pbsq_1 = extract_data_module(\'REPORT_DIR/reads_1_fastqc_data.txt\', \'Per base sequence quality\')\n pbsq_1$id = 1:length(pbsq_1$X.Base)\n \n-melt_pbsq_1 = filter(melt(pbsq_1, id=c(\'X.Base\', \'id\')), variable != \'X90th.Percentile\' & variable != \'X10th.Percentile\')\n+melt_pbsq_1 = filter(melt(pbsq_1, id=c(\'X.Base\', \'id\')), variable == \'Mean\')\n melt_pbsq_1$trim = \'before\'\n \n \n@@ -146,15 +148,17 @@\n pbsq_2 = extract_data_module(\'REPORT_DIR/reads_2_fastqc_data.txt\', \'Per base sequence quality\')\n pbsq_2$id = 1:length(pbsq_2$X.Base)\n \n-melt_pbsq_2 = filter(melt(pbsq_2, id=c(\'X.Base\', \'id\')), variable != \'X90th.Percentile\' & variable != \'X10th.Percentile\')\n+melt_pbsq_2 = filter(melt(pbsq_2, id=c(\'X.Base\', \'id\')), variable == \'Mean\')\n melt_pbsq_2$trim = \'after\'\n \n comb_pbsq = rbind(melt_pbsq_1, melt_pbsq_2)\n comb_pbsq$trim = factor(levels = c(\'before\', \'after\'), comb_pbsq$trim)\n+\n p = ggplot(data = comb_pbsq) +\n geom_line(mapping = aes(x = id, y = value, group = variable, color = variable)) +\n scale_x_continuous(breaks = pbsq_2$id, labels = pbsq_2$X.Base) + \n facet_grid(. ~ trim) + \n+ ylim(0, max(comb_pbsq$value) + 5) +\n theme(axis.text.x = element_text(angle=45))\n ggplotly(p)\n \n@@ -162,27 +166,249 @@\n \n ### Per tile sequence quality\n \n-```{r \'per tile sequence quality\'}\n+```{r \'per tile sequence quality\', fig.width=10}\n+## check if \'per tile sequence quality\' module exits or not\n+check_ptsq = grep(\'Per tile sequence quality\', readLines(\'REPORT_DIR/reads_1_fastqc_data.txt\'))\n+if (length(check_ptsq) > 0) {\n+ ## reads 1\n+ ptsq_1 = extract_data_module(\'REPORT_DIR/reads_1_fastqc_data.txt\', \'Per tile sequence quality\')\n+ ptsq_1$trim = \'before\'\n+ \n+ ## reads 2\n+ ptsq_2 = extract_data_module(\'REPORT_DIR/reads_2_fastqc_data.txt\', \'Per tile sequence quality\')\n+ ptsq_2$trim = \'after\'\n+ \n+ comb_ptsq = rbind(ptsq_1, ptsq_2)\n+ comb_ptsq$trim = factor(levels = c(\'before\', \'after\'), comb_ptsq$trim)\n+ comb_ptsq$Base = factor(levels = unique(comb_ptsq$Base), comb_ptsq$Base)\n+ \n+ # convert integers to charaters\n+ comb_ptsq$Tile = as.character(comb_ptsq$X.Tile)\n+ \n+ p = ggplot(data = comb_ptsq, aes(x = Base, y = Tile, fill = Mean)) +\n+ geom_raster() + \n+ facet_grid(. ~ trim) + \n+ xlab(\'Position in read (bp)\') + \n+ ylab(\'\') +\n+ theme(axis.text.x = element_text(angle=45))\n+ ggplotly(p)\n+} else {\n+ print(\'No "per tile sequence quality" data\')\n+}\n+\n+\n+```\n+\n+### Per sequence quality score\n+\n+```{r \'Per sequence quality score\', fig.width=10}\n ## reads 1\n-ptsq_1 = extract_data_module(\'REPORT_DIR/reads_1_fastqc_data.txt\', \'Per tile sequence quality\')'..b"xt', 'Per base N content')\n+pbNc_1$id = 1:length(pbNc_1$X.Base)\n+pbNc_1$trim = 'before'\n+\n+## reads 2\n+pbNc_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Per base N content')\n+pbNc_2$id = 1:length(pbNc_2$X.Base)\n+pbNc_2$trim = 'after'\n+\n+comb_pbNc = rbind(pbNc_1, pbNc_2)\n+comb_pbNc$trim = factor(levels = c('before', 'after'), comb_pbNc$trim)\n+\n+p = ggplot(data = comb_pbNc, aes(x = id, y = N.Count)) +\n+ geom_line(color = 'red') +\n+ scale_x_continuous(breaks = pbNc_2$id, labels = pbNc_2$X.Base) + \n+ facet_grid(. ~ trim) +\n+ ylim(0, 1) + \n+ xlab('N-Count') +\n+ ylab('') + \n theme(axis.text.x = element_text(angle=45))\n ggplotly(p)\n ```\n \n \n+### Sequence Length Distribution\n+\n+```{r 'Sequence Length Distribution', fig.width=10}\n+## reads 1\n+sld_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Sequence Length Distribution')\n+sld_1$id = 1:length(sld_1$X.Length)\n+sld_1$trim = 'before'\n+\n+## reads 2\n+sld_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Sequence Length Distribution')\n+sld_2$id = 1:length(sld_2$X.Length)\n+sld_2$trim = 'after'\n+\n+comb_sld = rbind(sld_1, sld_2)\n+comb_sld$trim = factor(levels = c('before', 'after'), comb_sld$trim)\n+\n+p = ggplot(data = comb_sld, aes(x = id, y = Count)) +\n+ geom_line(color = 'red') +\n+ scale_x_continuous(breaks = sld_2$id, labels = sld_2$X.Length) + \n+ facet_grid(. ~ trim) +\n+ xlab('Sequence Length (bp)') +\n+ ylab('') + \n+ theme(axis.text.x = element_text(angle=45))\n+ggplotly(p)\n+```\n+\n+### Sequence Duplication Levels\n+\n+```{r 'Sequence Duplication Levels', fig.width=10}\n+## reads 1\n+sdl_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Sequence Duplication Levels', header = FALSE, comment.char = '#')\n+names(sdl_1) = c('Duplication_Level', 'Percentage_of_deduplicated', 'Percentage_of_total')\n+sdl_1$id = 1:length(sdl_1$Duplication_Level)\n+\n+melt_sdl_1 = melt(sdl_1, id=c('Duplication_Level', 'id'))\n+melt_sdl_1$trim = 'before'\n+\n+\n+## reads 2\n+sdl_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Sequence Duplication Levels', header = FALSE, comment.char = '#')\n+names(sdl_2) = c('Duplication_Level', 'Percentage_of_deduplicated', 'Percentage_of_total')\n+sdl_2$id = 1:length(sdl_2$Duplication_Level)\n+\n+melt_sdl_2 = melt(sdl_2, id=c('Duplication_Level', 'id'))\n+melt_sdl_2$trim = 'after'\n+\n+comb_sdl = rbind(melt_sdl_1, melt_sdl_2)\n+comb_sdl$trim = factor(levels = c('before', 'after'), comb_sdl$trim)\n+\n+p = ggplot(data = comb_sdl, aes(x = id, y = value, color = variable)) +\n+ geom_line() +\n+ scale_x_continuous(breaks = sdl_2$id, labels = sdl_2$Duplication_Level) +\n+ facet_grid(. ~ trim) +\n+ xlab('Sequence Duplication Level') +\n+ ylab('') + \n+ theme(axis.text.x = element_text(angle=45))\n+ggplotly(p)\n+```\n+\n+### Adapter Content\n+\n+```{r 'Adapter Content', fig.width=10}\n+## reads 1\n+ac_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Adapter Content')\n+ac_1$id = 1:length(ac_1$X.Position)\n+\n+melt_ac_1 = melt(ac_1, id=c('X.Position', 'id'))\n+melt_ac_1$trim = 'before'\n+\n+## reads 2\n+ac_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Adapter Content')\n+ac_2$id = 1:length(ac_2$X.Position)\n+\n+melt_ac_2 = melt(ac_2, id=c('X.Position', 'id'))\n+melt_ac_2$trim = 'after'\n+\n+comb_ac = rbind(melt_ac_1, melt_ac_2)\n+comb_ac$trim = factor(levels = c('before', 'after'), comb_ac$trim)\n+\n+p = ggplot(data = comb_ac, aes(x = id, y = value, color = variable)) +\n+ geom_line() +\n+ facet_grid(. ~ trim) +\n+ xlim(min(comb_ac$id), max(comb_ac$id)) + \n+ ylim(0, 1) +\n+ xlab('Position in read (bp)') +\n+ ylab('')\n+ggplotly(p)\n+```\n+\n+### Kmer Content {.tabset}\n+\n+#### Before\n+\n+```{r 'Kmer Content (before)', fig.width=10}\n+kc_1 = extract_data_module('REPORT_DIR/reads_1_fastqc_data.txt', 'Kmer Content')\n+knitr::kable(kc_1)\n+```\n+\n+#### After\n+```{r 'Kmer Content (after)', fig.width=10}\n+kc_2 = extract_data_module('REPORT_DIR/reads_2_fastqc_data.txt', 'Kmer Content')\n+knitr::kable(kc_2)\n+```\n+\n \n # Session Info\n \n" |